一种可考虑频率约束的水平集拓扑优化方法

1.本发明属于拓扑领域，特别涉及一种可考虑频率约束的水平集拓扑优化方法。


背景技术：

2.水平集法(即level set法)是当前主要的集中拓扑优化方法之一，其主要优点是优化结果边界光滑，没有变密度法(即simp法)、或者双向演进法(beso)法等优化结果的边界锯齿化问题或者出现棋盘格现象。目前已有的水平集法一般以结构的柔顺度最小化为目标函数，以结构的体积分数(即优化后实体材料体积与原设计域体积之比)为约束条件。在实际工程中，部分情况下需要考虑优化后的结构间的共振问题，也即需要考虑频率约束问题。


技术实现要素：

3.有鉴于现有技术的上述缺陷，本发明的目的在于提供一种可考虑频率约束的水平集拓扑优化方法，以解决上述问题。
4.本发明提供如下的技术方案：
5.一种可考虑频率约束的水平集拓扑优化方法，该水平集方法是通过记录孔洞中心坐标的方法记录孔洞位置，通过计算孔洞中心和各网格节点横纵坐标的距离，建立水平集函数；将各个节点之间的横坐标和纵坐标距离与mq样条曲线基函数相结合，公式如下：
[0006][0007]
得到了各个中心点的水平集径向基函数：
[0008][0009]
其中αi(t)是广义展开系数，p(x,t)是一个随时间变化的一阶多项式系数：
[0010]
p(x,t)＝p0(t) p1(t)x p2(t)y
[0011]
为了保证水平集函数的rbf插值的唯一解，展开系数必须服从以下约束：
[0012][0013]
结合公式公式p(x,t)＝p0(t) p1(t)x p2(t)y和公式
[0014]
可重写为矩阵形式：
[0015]
gα(t)＝φ(t)；
[0016]
其中，
[0017]
[0018][0019][0020]
则广义展开系数就可转换成：
[0021]
α(t)＝g-1
φ(t)
[0022]
相对应公式的水平集函数可以简写为：
[0023]
φ(x,t)＝g(x)
t
α(t)
[0024]
这里的g(x)：
[0025]
g(x)＝{g1(x)
…gn
(x),1,x,y}
t
[0026]
将公式g(x)＝{g1(x)
…gn
(x),1,x,y}
t
带入hamilton-jacobi方程中，得到基于rbf展开系数的控制方程
[0027][0028]
其中
[0029][0030][0031][0032]
对于这个依赖时间的插值问题，必须结合公式相关约束，才能保证mq样条的条件正确性。在拓扑优化问题中，由于结点处初始水平集函数值是给定的，且矩阵g在理论上是可逆的，因此求解n 3个线性方程组，即可得到初始广义展开系数α(t0)：
[0033]
α(t0)＝g-1
φ(t0)
[0034]
采用一阶前向欧拉法求解，其近似值为：
[0035]
α(t
i 1
)＝α(ti) δtg-1
b(α(ti),ti)
[0036]
其中ti为第i个时间步长，δt为时间步长。在实际实现过程中，为了防止公式α(t
i 1
)＝α(ti) δtg-1
b(α(ti),ti)可能会造成不必要的数值不稳定性，将其速度重申为：
[0037]
b(ti)＝{vn(x1,ti)
…vn
(xn,ti),0,0,0}
t
[0038]
为了防止边界附近的|
▽
φ|值过大或者过小，将整个设计域中的水平集函数φ更新为：
[0039][0040]
由于φ和α关系是线性的，所以将α更新成αu：
[0041][0042]
为了避免水平集函数的无界增长，在公式α(t
i 1
)＝α(ti) δtg-1
b(α(ti),ti)的演化方案中引入一个近似的δ(φ)函数：
[0043][0044]
因此，更新方程修改为：
[0045][0046]
其中可以表示为ti时的速度矢量，它可以写成：
[0047][0048]
为了简单而又不失一般性，接下来讨论在体积约束下静载线弹性结构的柔度最小化问题。下面是基于数学模型：
[0049][0050][0051]
上述公式的符号定义为：
[0052][0053][0054]
式中j(φ)为目标函数，u为位移场，ε为线性化应变张量，c为胡克弹性张量，u0，τ和b分别表示给定位移、牵引力和体力常数值。线性弹性平衡方程以能量形式a(u,v,φ)和载荷形式l(v,φ)分别表示，其中v表示运动允许位移场空间u中的虚位移场。g(φ)是限制材料使用的约束，v
max
是设计域d的最大允许体积。h(φ)表示物质在水平集函数值φ上某一
点存在，其定义如下：
[0055][0056]
对于柔度最小问题可以导出法向加速度：
[0057]vn
＝ε(u):c:ε(u)-λ
[0058]
式中λ为拉格朗日乘子，用于处理体积分数约束。或者采用下面增广拉格朗日更新方案：
[0059][0060]
其中μ和γk为第k次迭代的参数，γk的更新方案为：
[0061]
γ
k 1
＝min(γk δγ,γ
max
)k＞nr[0062]
式中δγ为增量，γ
max
为γ的上限。而基于水平集的优化方法中初始设计的体积分数通常不满足规定的体积分数，因此在第一次迭代中，体积约束被放宽如下：
[0063][0064]
以上过程仅考虑了体积分数约束，在工程中为了保证结构变形不超限，通常需要设置频率约束，也即结构整体频率维持在频率限值之上。因此本发明将原有的水平集法的拉格朗日函数进行拓展，增加频率约束项，具体如下：
[0065]
本发明在原有的水平集目标函数基础上添加了相应的频率拉格朗日乘子和相应的频率约束，如下公式：
[0066][0067]
公式中：c(x)为目标函数，λv为体积拉格朗日乘子，为体积约束，μ2为频率拉格朗日乘子，μ2(ω
*-ω)为频率约束，ω为结构第一阶的自然频率，ω
*
为其下界。
[0068]
当结构第一阶的自然频率小于其下界值时，作者采用如下公式对频率进行相应惩罚：
[0069][0070]
为了进一步控制增量速度，作者做出了进一步的控制，如下:
[0071][0072]
频率拉格朗日乘子更新公式如下：
[0073][0074]
本发明的有益效果是：在原水平集法的基础上，提出了一种通过构造拉格朗日方程以及更新拉格朗日乘子的方法，使水平集法能够在拓扑优化过程中考虑频率约束，并通过算例验证了该方法的有效性。
附图说明
[0075]
图1是本发明实施例中的平面悬臂梁示意图；
[0076]
图2是本发明实施例中约束情况1的示意图；
[0077]
图3是本发明实施例中约束情况2的示意图；
[0078]
图4是本发明实施例中约束情况3的示意图；
[0079]
图5是本发明实施例中约束情况4的示意图。
具体实施方式
[0080]
下面结合附图对本发明的实施例作详细说明，下述的实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
[0081]
实施例
[0082]
参见附图1-5，本发明实施例提供的可考虑频率约束的水平集拓扑优化方法，该水平集方法是通过记录孔洞中心坐标的方法记录孔洞位置，通过计算孔洞中心和各网格节点横纵坐标的距离，建立水平集函数；将各个节点之间的横坐标和纵坐标距离与mq样条曲线基函数相结合，公式如下：
[0083][0084]
得到了各个中心点的水平集径向基函数：
[0085][0086]
其中αi(t)是广义展开系数，p(x,t)是一个随时间变化的一阶多项式系数：
[0087]
p(x,t)＝p0(t) p1(t)x p2(t)y
[0088]
为了保证水平集函数的rbf插值的唯一解，展开系数必须服从以下约束：
[0089][0090]
结合公式公式p(x,t)＝p0(t) p1(t)x p2(t)y和公式可重写为矩阵形式：
[0091]
gα(t)＝φ(t)；
[0092]
其中，
[0093][0094]
[0095][0096]
则广义展开系数就可转换成：
[0097]
α(t)＝g-1
φ(t)
[0098]
相对应公式的水平集函数可以简写为：
[0099]
φ(x,t)＝g(x)
t
α(t)
[0100]
这里的g(x)：
[0101]
g(x)＝{g1(x)
…gn
(x),1,x,y}
t
[0102]
将公式g(x)＝{g1(x)
…gn
(x),1,x,y}
t
带入hamilton-jacobi方程中，得到基于rbf展开系数的控制方程
[0103][0104]
其中
[0105][0106][0107][0108]
对于这个依赖时间的插值问题，必须结合公式相关约束，才能保证mq样条的条件正确性。在拓扑优化问题中，由于结点处初始水平集函数值是给定的，且矩阵g在理论上是可逆的，因此求解n 3个线性方程组，即可得到初始广义展开系数α(t0)：
[0109]
α(t0)＝g-1
φ(t0)
[0110]
采用一阶前向欧拉法求解，其近似值为：
[0111]
α(t
i 1
)＝α(ti) δtg-1
b(α(ti),ti)
[0112]
其中ti为第i个时间步长，δt为时间步长。在实际实现过程中，为了防止公式α(t
i 1
)＝α(ti) δtg-1
b(α(ti),ti)可能会造成不必要的数值不稳定性，将其速度重申为：
[0113]
b(ti)＝{vn(x1,ti)
…vn
(xn,ti),0,0,0}
t
[0114]
为了防止边界附近的|
▽
φ|值过大或者过小，将整个设计域中的水平集函数φ更新为：
[0115][0116]
由于φ和α关系是线性的，所以将α更新成αu：
[0117][0118]
为了避免水平集函数的无界增长，在公式α(t
i 1
)＝α(ti) δtg-1
b(α(ti),ti)的演化方案中引入一个近似的δ(φ)函数：
[0119][0120]
因此，更新方程修改为：
[0121][0122]
其中可以表示为ti时的速度矢量，它可以写成：
[0123][0124]
为了简单而又不失一般性，接下来讨论在体积约束下静载线弹性结构的柔度最小化问题。下面是基于数学模型：
[0125][0126][0127]
上述公式的符号定义为：
[0128][0129][0130]
式中j(φ)为目标函数，u为位移场，ε为线性化应变张量，c为胡克弹性张量，u0，τ和b分别表示给定位移、牵引力和体力常数值。线性弹性平衡方程以能量形式a(u,v,φ)和载荷形式l(v,φ)分别表示，其中v表示运动允许位移场空间u中的虚位移场。g(φ)是限制材料使用的约束，v
max
是设计域d的最大允许体积。h(φ)表示物质在水平集函数值φ上某一点存在，其定义如下：
[0131]
[0132]
对于柔度最小问题可以导出法向加速度：
[0133]vn
＝ε(u):c:ε(u)-λ
[0134]
式中λ为拉格朗日乘子，用于处理体积分数约束。或者采用下面增广拉格朗日更新方案：
[0135][0136]
其中μ和γk为第k次迭代的参数，γk的更新方案为：
[0137]
γ
k 1
＝min(γk δγ,γ
max
)k＞nr[0138]
式中δγ为增量，γ
max
为γ的上限。而基于水平集的优化方法中初始设计的体积分数通常不满足规定的体积分数，因此在第一次迭代中，体积约束被放宽如下：
[0139][0140]
以上过程仅考虑了体积分数约束，在工程中为了保证结构变形不超限，通常需要设置频率约束，也即结构整体频率维持在频率限值之上。因此本发明将原有的水平集法的拉格朗日函数进行拓展，增加频率约束项，具体如下：
[0141]
本发明在原有的水平集目标函数基础上添加了相应的频率拉格朗日乘子和相应的频率约束，如下公式：
[0142][0143]
公式中：c(x)为目标函数，λv为体积拉格朗日乘子，为体积约束，μ2为频率拉格朗日乘子，μ2(ω
*-ω)为频率约束，ω为结构第一阶的自然频率，ω
*
为其下界。
[0144]
当结构第一阶的自然频率小于其下界值时，作者采用如下公式对频率进行相应惩罚：
[0145][0146]
为了进一步控制增量速度，作者做出了进一步的控制，如下:
[0147][0148]
频率拉格朗日乘子更新公式如下：
[0149][0150]
下面结合算例，说明该方法的有效性：
[0151]
如图1所示，以一个80
×
60的平面悬臂梁1为例，其左侧有一个固定约束，右侧中点处有一个f＝-100的集中荷载，如图1所示。材料性能：固体材料杨氏弹性模量e＝1，空隙材料的弹性模量e＝10，泊松比υ＝0.3，质量密度7.8
×
10-3
，规定的体积分数设为50％。
[0152]
表1给出了四种不同频率限制约束情况和其分别对应的优化结果，其中表1中的表示结构第一阶的自然频率。
[0153]
表1不同频率限制约束及相应的优化设计结果
[0154][0155]
如图2-5所示的约束情况示意图可以看出，不同频率限制约束情况都能被很好的满足，验证了所添加的频率约束的有效性。
[0156]
以上详细描述了本发声明的较佳具体实施例。应当理解，本领域的普通技术无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的试验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。