CN107315872A - 一种高效的结构频率响应拓扑优化方法 - Google Patents

Abstract

本发明属于结构拓扑优化设计相关技术领域，其公开了一种高效的结构频率响应拓扑优化方法，其包括以下步骤：(1)将待优化结构的动力学模型中的标准水平集函数中的时间和空间两个耦合变量解耦，同时将与时间相关的水平集函数表示为矩阵乘积形式；(2)将与时间相关的水平集函数的偏微分方程转化为常微分方程，进而得到新的线性系统，并求解获得与时间相关的水平集函数；(3)对宏观结构进行有限元分析，进而计算结构优化问题的目标函数与约束函数；(4)计算步骤(3)所得的目标函数与约束函数关于设计变量的敏度，进而更新设计变量后，判断所述目标函数是否收敛。该方法采用离散小波变换技术对插值矩阵进行再压缩，提高了效率，降低了成本。

Description

一种高效的结构频率响应拓扑优化方法

技术领域

本发明属于结构拓扑优化设计相关技术领域，更具体地，涉及一种高效的结构频率响应拓扑优化方法。

背景技术

基于水平集的边界描述技术具有独特的特点，例如边界形状光滑、清晰，能够方便地通过边界的融合与分裂灵活地描述其拓扑和形状的变化。由于水平集函数没有显式的解析解，整个设计域需要用矩形网格离散，并采用有限差分方法求解水平集方程，使得传统水平集方法应用于结构拓扑优化问题时尚存在多个缺陷。传统水平集方法拓扑形状优化存在多个问题，如求解Hamilton-Jacobi偏微分方程的逆风差分格式步长受Courant-Friedrichs-Lewy(CFL)条件限制，求解速度慢，优化效率不高；拓扑形状优化过程中需要不断对水平集函数进行耗时的周期性初始化，以保证求解的精度和数值稳定性；结构边界演化必须通过求解Hamilton-Jacobi偏微分方程来实现，不能同优化领域中成熟高效的优化算法相结合(如最优化准则法、数学规划法等)。因此，传统水平集函数离散计算的数值问题严重影响了水平集方法在结构优化应用中的优势。

此外，针对结构动力学优化设计问题，结构动力学拓扑优化主要包括特征值拓扑优化和频率响应拓扑优化，结构频率响应拓扑优化通常更为困难。首先，频率响应优化问题具有高度非线性和目标函数非凸的特点，降低了优化过程的稳定性与收敛性；其次，频率响应敏度求解困难，同时涉及频率和振型的敏度分析；再次，振动问题的有限元平衡方程求解复杂。这些困难导致结构频率响应拓扑优化发展缓慢，研究成果相对较少，目前在该领域研究多为单一频率激励下的动力响应优化，而在实际应用中，如火箭、导弹的雷达、跟踪仪等伺服系统、航空航天的精密设备、汽车减震降噪的结构与系统等通常会受到一定带宽范围内的频带激励，这类结构频率响应优化问题进一步增加了求解难度：一方面结构优化设计需要避开一个频带范围以避免共振，在优化过程中结构固有频率落入目标频带范围时可能产生多个相应峰值，影响迭代的稳定性；另一方面，频带激励下结构响应拓扑优化的目标函数为积分形式，这意味着采取高斯积分等方法进行数值处理时，必须多次求解复杂的振动方程，效率较低，极大地增加了计算成本，在问题规模较大时极易导致优化失败。相应地，本领域存在着发展一种效率较高的结构频率响应拓扑优化方法的技术需求。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种高效的结构频率响应拓扑优化方法，其利用高斯径向基函数插值技术使Hamilton-Jacobi偏微分方程转化为更易求解的常微分方程，同时采用离散小波变换技术对由高斯径向基函数构建的插值矩阵进行压缩，从而获得一种极其稀疏的插值系统，减少了优化过程的处理时间及其所需的计算机存储空间，解决了计算成本较高、效率较低的问题。

为实现上述目的，本发明提供了一种高效的结构频率响应拓扑优化方法，其包括以下步骤：

(1)采用高斯径向基函数将待优化结构的动力学模型中的标准水平集函数中的时间和空间两个耦合变量解耦，同时将与时间相关的水平集函数表示为高斯径向基函数构建的插值矩阵与扩展系数向量的矩阵乘积形式；

(2)将与时间相关的水平集函数的Hamilton-Jacobi偏微分方程转化为常微分方程，并对所述插值矩阵、所述扩展系数向量及所述与时间相关的水平集函数进行离散小波分解，以对高斯径向基函数构建的插值矩形进行压缩，进而得到新的线性系统，并求解获得与时间相关的水平集函数；

(3)根据求解获得的与时间相关的水平集函数对宏观结构进行有限元分析，进而计算结构优化问题的目标函数与约束函数；

(4)计算步骤(3)所得的目标函数与约束函数关于设计变量的敏度，进而更新设计变量后，判断所述目标函数是否收敛，若收敛，则输出待优化结构的最优宏观结构构型；否则，转至步骤(3)。

进一步地，步骤(1)之前还包括建立待优化结构的动力学优化模型的步骤，具体包括构造参数化水平集的结构局部频率响应优化模型和动柔度最小化问题优化模型。

进一步地，高斯径向基函数构建的插值矩阵A经小波分解后转换成矩阵步骤(2)还包括引入阈值以消除所述矩阵中的噪声元素，以得到压缩后的稀疏矩阵的步骤。

进一步地，所述阈值的计算公式如下：

式中，代表插值矩阵中的元素，代表插值矩阵中元素绝对值的平均值，κ为阀值参数。

进一步地，所述阀值参数κ的取值范围为1≤κ≤10。

进一步地，所述设计变量为扩展系数；所述目标函数与所述约束函数关于设计变量的敏度计算是通过形状导数与伴随变量法进行的。

进一步地，设计变量的更新是通过优化准则法进行的。

进一步地，与时间相关的水平集函数由高斯径向基函数构建的插值矩阵与扩展系数向量构成的矩阵乘积形式为：

Aα＝Φ

式中，Φ为与时间相关的水平集函数对应的矩阵，α为仅与时间相关的扩展系数向量，A为仅与空间变量相关的高斯径向基函数构建的插值矩阵。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，本发明提供的高效的结构频率响应拓扑优化方法主要具有以下有益效果：

(1)采用高斯径向基函数将动力学模型中的标准水平集函数中的时间和空间两个耦合变量解耦，同时将与时间相关的水平集函数表示为高斯径向基函数构建的插值矩阵与扩展系数向量的矩阵乘积形式，使得合适的高斯径向基函数保证了水平集函数的光滑性和连续性要求，无需反复对水平集函数进行重新初始化，减少了时间消耗；

(2)所述结构频率拓扑优化设计方法将与时间相关的水平集函数的Hamilton-Jacobi偏微分方程转化为常微分方程，避免了直接求解Hamilton-Jacobi偏微分方程所带来的数值问题，更易求解，简化了计算，减少了计算时间；

(3)采用小波分解技术对高斯径向基函数构建的插值矩阵进行压缩，从而获得一种极其稀疏的插值系统，减少优化过程的处理时间及其所需的计算机存储空间，极大地提高了效率。

附图说明

图1是本发明较佳实施方式提供的高效的结构频率响应拓扑优化方法的流程图；

图2是采用图1中的高效的结构频率响应拓扑优化设计方法来进行拓扑优化的二维结构的设计域示意图；

图3中的(a)、(b)、(c)、(d)、(e)、(f)分别是图2中的二维结构针对局部频率响应优化问题的拓扑优化的过程图；

图4中的(a)、(b)分别为图2中的二维结构的目标函数与约束条件收敛曲线、结构频率响应曲线。

图5是图2中的二维结构针对局部频率响应优化问题、采用不同插值方式及阈值参数时的最优结构拓扑形式；

图6中的(a)、(b)分别为图2中的二维结构针对动柔度的最优结构拓扑及3D水平集函数；

图7中的(a)、(b)分别是图2中的二维结构针对动柔度问题的目标函数与约束条件收敛曲线、结构频率响应曲线。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

请参阅图1，本发明较佳实施方式提供的高效的结构频率响应拓扑优化方法，所述结构频率响应拓扑优化方法主要包括以下步骤：

步骤一，建立待优化结构的动力学优化模型，具体包括构造参数化水平集的结构局部频率响应优化模型和动柔度最小化问题优化模型。

局部频率响应优化模型为：

式中，J_l为优化模型的目标函数，表达为一种双重积分形式；[ω_s,ω_e]表示激励载荷的振动频带范围；Ω_r∈Ω为所考虑的响应输出端，即待优化的局部区域；u和v为结构的频率响应，是关于频率ω的函数；u₀为响应许用值；H为Heaviside函数；Φ为水平集函数，G为材料体积分数约束；V_max为最大的材料使用量；i＝1,2,…,N为设计变量的编号；α_max和α_min分别表示设计变量的上下限；U定义了运动学允许的位移空间；j为虚数；k(u,v,Φ)为结构刚度半双线性形式；c(u,v,Φ)为结构阻尼半双线性形式；m(u,v,Φ)为结构质量半双线性形式；l(v,Φ)为载荷半线性形式。在频带范围内对目标函数的积分可采用高斯积分法进行计算。

动柔度最小化问题优化模型：

假设p代表结构外部激励的载荷大小，基于参数化水平集的材料体积分数约束下，结构全局频率响应最小化的拓扑优化模型为：

可以看到，全局频率响应拓扑优化模型采用了与局部频率响应拓扑优化模型不同的目标函数J_g，该目标函数J_g即为结构的动态柔度。

步骤二，采用高斯径向基函数将所述动力学模型中的标准水平集函数中的时间和空间两个耦合变量解耦，同时将与时间相关的水平集函数表示为高斯径向基函数构建的插值矩阵与扩展系数向量的矩阵乘积形式。具体地，所述的矩阵乘积形式为：

其中，x为关于位置的变量，t为关于时间的变量。高斯径向基函数构建的插值矩阵为：

扩散系数构成的向量为：

α(t)＝[α₁(t),α₂(t),...,α_N(t)]^T∈R^N (5)

为了阐述方便，可将公式(3)表述为以下矩阵形式：

Aα＝Φ (6)

式中，Φ＝[Φ₁,Φ₂,…,Φ_N]^T代表N个插值点处的水平集函数值；A为一个由高斯径向基函数构成的N×N可逆矩阵：

此时，与时间相关的水平集函数由两部分构成：(1)仅与空间变量相关的高斯径向基函数；(2)仅与时间变量相关的扩展系数α。原本时间变量和空间变量耦合的水平集函数Φ(x,t)通过高斯径向基函数插值实现了时间变量与空间变量的解耦，可以看到插值矩阵A仅与所选取的高斯径向基函数及控制点的位置坐标相关，也就是说在插值控制点固定的前提下，矩阵A由一系列常量构成，故其在迭代过程中仅需计算一次。

步骤三，将与时间相关的水平集函数的Hamilton-Jacobi(雅可比)偏微分方程转化为常微分方程，并对所述插值矩阵、所述扩展系数向量及所述与时间相关的水平集函数进行离散小波分解以对高斯径向基函数构建的插值矩形进行压缩，进而得到新的线性系统，求解获得与时间相关的水平集函数。步骤三具体包括以下子步骤：

(31)将与时间相关的水平集函数的Hamilton-Jacobi偏微分方程转化为常微分方程形式(如下式所示)，避免直接求解Hamilton-Jacobi偏微分方程所带来的数值问题。

式中，▽为微分算子。

由公式(8)可建立速度场v_n与扩展系数α间的数学关系：

可以看到，公式(9)右边的所有子项均需要在全部插值控制点上进行计算，因此，相较于传统水平集方法仅计算边界上的速度场，初始边界的速度场v_n被扩展到了整个设计域中。

至此，已将传统的水平集方法转化为一种参数化的形式，在保留传统水平集方法优点的同时，能够有效克服其应用缺陷。具体来讲，高斯径向基函数插值使得原本由偏微分方程驱动的水平集函数改变为相对容易的常微分方程驱动，从而避免了直接求解Hamilton-Jacobi偏微分方程所带来的数值问题，具体地，(1)合适的高斯径向基函数保证了水平集函数的光滑性和连续性要求，无需反复对水平集函数进行重新初始化；(2)步长不受迎风差分格式中的CFL条件限制，不会被最小网格尺寸所影响，因此设计域离散时可以使用更加精细的网格以提高精度，求解大规模问题时效率更高，尤其是针对动力学复杂的计算，可大量减少计算时间；(3)如前所述，速度场v_n被自动地扩展到整个设计域中；(4)基于梯度的优化算法能够得以直接应用。

(32)采用离散小波分解技术对所述插值矩阵、所述扩展系数向量及与时间相关的水平集函数进行小波分解，以对高斯径向基函数构建的插值矩阵A进行压缩，减少优化过程的处理时间及所需的计算机存储空间。

具体地，引入正交矩阵W，针对扩展系数向量α和与时间相关的水平集函数矩阵Φ的小波变换过程为和矩阵A的变换过程则可类似地表述为小波分解后的矩阵突出了矩阵A中的数据特征，能够利用中的少量但重要的元素来反映A中元素的内在联系。

(33)引入合适的阀值以消除矩阵中的噪声元素，得到压缩后的稀疏矩阵通过矩阵中存在的少量非零数据元素即可实现对原始矩阵A中全部数据的操作。阀值计算公式如下所示：

式中，代表插值矩阵中的元素，代表元素绝对值的平均值，κ为阀值参数。

(34)阀值可通过改变阈值参数κ≥0的值进行调整，为了得到较高的插值系统稀疏度和较高的计算效率，阀值参数κ的取值范围为1≤κ≤10，所得到的插值矩阵A中零元素的数量占比达到98％以上。

(35)利用小波变换后的向量和矩阵以及压缩后的矩阵来得到新的线性系统：由于矩阵的稀疏性，参数化水平集方法中的扩展系数α可利用重构操作得到，水平集函数Φ则可利用得到。

步骤四，根据求解获得的与时间相关的水平集函数对宏观结构进行有限元分析，进而构建结构优化问题的目标函数与约束函数。具体地，根据求解获得的与时间相关的水平集函数对宏观结构进行有限元分析，以获取结构整体位移场U，进而计算结构优化的目标函数与约束函数。

步骤五，计算步骤四所得到的目标函数与约束函数关于设计变量的敏度，进而更新设计变量后，判断所述目标函数是否收敛，如果收敛，则输出最优宏观结构构型；否则，转至步骤四。

具体地，首先，采用形状导数与伴随变量法计算所述目标函数与约束函数关于设计变量的敏度，其中所述设计变量为扩展系数。目标函数关于设计变量的敏度计算公式如下：

局部频率响应：

式中，β(u,w,Φ)＝pw-(1+jωA)E_ijklε_kl(u)ε_ij(w)-(jωB-ω²)ρu_iw_i+div(τwn)。u,w为结构的频率响应；p为动态载荷；A、B为Rayleigh阻尼参数；E_ijkl为材料弹性张量；ρ表示材料的密度ε_kl(u),ε_ij(w)为应变；u_i,w_i为单元的频率响应；div为散度算子；τ为结构的边界牵引力分量；n为法向量；为高斯径向基函数；δ(Φ)为Heaviside函数的导数。

动柔度：

式中，γ(u,u,Φ)＝2pu-(1+jωA)E_ijklε_kl(u)ε_ij(u)-(jωB-ω²)ρu_iu_i+2div(τun)。

之后，基于目标函数与约束函数及其对应的灵敏度，采用优化准则法更新设计变量。具体地，通过引入Lagrange乘子将带约束优化问题转化为无约束优化问题：

式中，Λ、λ、分别为体积约束、平衡方程、设计变量下限、设计变量上限的Lagrange乘子。

基于Kuhn-Tucker条件，建立Lagrange函数L的驻值条件。进而建立设计变量更新条件。

最后，判断所述目标函数是否收敛，如果收敛，则以水平集函数的形式输出待优化结构的最优宏观结构构型；否则，转至步骤四。

请参阅图2至图7，以下以二维结构的频率响应优化设计来进一步说明本发明。如图2所示，所述的二维结构的设计域的尺度为140cm×20cm×1cm，材料自身弹性模量E＝1000GPa，泊松比等于0.3，结构的左右端面固定，结构中心点处有F＝1000e^iωt竖直向下的简谐激励；给定频率带为Ω_freq＝[0Hz,100Hz]。优化目标为负载点的响应最小以及全局动柔度最小，限定材料使用量为50％。

如图3所示，在给定阈值参数κ＝1时，对应的结构最优拓扑形式迭代变化情况，从整个过程来看，结构的形式会不断的变化至最优形式。图4给出了局部频率响应下、目标函数以及约束条件的收敛曲线，可以看出迭代非常稳定，在100步以内基本上可以达到最优解，这说明了本发明提供的结构频率响应拓扑优化方法的高效率性，优化后的结构特征频率要大于初始的特征频率，这也说明了该结构频率响应拓扑优化方法的合理性，使得结构的稳定性增强。图5给出的是阈值参数κ在不同值的情况下，结构的最优形式变化曲线，从图中可以看出，结构的最优形式以及3D水平集函数图基本上没有变化，这说明了该结构频率响应拓扑优化方法的稳定性，在阈值取值的一定范围内对结构的优化影响不大。图6与图7同样给出了针对动柔度问题下的结构优化设计结果，可以看出，基本上与频率响应问题得到的结论保持一致，即优化后的设计较之初始设计，在激励点处的频率响应有显著的降低，意味着优化设计能够有效地起到减振效果。此外，优化后的设计具有更高的结构基频，且处于所考虑的激励频带范围之外，从而避免了频带范围内产生共振而导致的结构失稳。

本发明提供的高效的结构频率响应拓扑优化方法，其采用高斯径向基函数插值水平集函数来实现水平集函数的参数化表达，同时采用离散小波变换技术来实现对原结构优化设计系统的有效数据提取，大大减少了数据冗余度，实现结构优化效率的显著提高，能够很好地解决结构频率响应优化问题，优化设计结果重量轻、边界清晰、动态性能优异。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种高效的结构频率响应拓扑优化方法，其特征在于，该方法包括以下步骤：

(1)采用高斯径向基函数将待优化结构的动力学模型中的标准水平集函数中的时间和空间两个耦合变量解耦，同时将与时间相关的水平集函数表示为高斯径向基函数构建的插值矩阵与扩展系数向量的矩阵乘积形式；

(2)将与时间相关的水平集函数的Hamilton-Jacobi偏微分方程转化为常微分方程，并对所述插值矩阵、所述扩展系数向量及所述与时间相关的水平集函数进行离散小波分解，以对高斯径向基函数构建的插值矩形进行压缩，进而得到新的线性系统，并求解获得与时间相关的水平集函数；

(3)根据求解获得的与时间相关的水平集函数对宏观结构进行有限元分析，进而计算结构优化问题的目标函数与约束函数；

(4)计算步骤(3)所得的目标函数与约束函数关于设计变量的敏度，进而更新设计变量后，判断所述目标函数是否收敛，若收敛，则输出待优化结构的最优宏观结构构型；否则，转至步骤(3)。

2.如权利要求1所述的高效的结构频率响应拓扑优化方法，其特征在于：步骤(1)之前还包括建立待优化结构的动力学优化模型的步骤，具体包括构造参数化水平集的结构局部频率响应优化模型和动柔度最小化问题优化模型。

3.如权利要求1所述的高效的结构频率响应拓扑优化方法，其特征在于：高斯径向基函数构建的插值矩阵A经小波分解后转换成矩阵步骤(2)还包括引入阈值以消除所述矩阵中的噪声元素，以得到压缩后的稀疏矩阵的步骤。

4.如权利要求3所述的高效的结构频率响应拓扑优化方法，其特征在于：所述阈值的计算公式如下：

<mrow> <mi>t</mi> <mi>h</mi> <mi>r</mi> <mi>e</mi> <mi>s</mi> <mi>h</mi> <mi>o</mi> <mi>l</mi> <mi>d</mi> <mrow> <mo>(</mo> <msub> <mover> <mi>a</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mover> <mi>a</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mo>|</mo> <msub> <mover> <mi>a</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>|</mo> <mo>&amp;GreaterEqual;</mo> <mi>&amp;kappa;</mi> <mover> <mi>m</mi> <mo>&amp;OverBar;</mo> </mover> <mo>;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>0</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>o</mi> <mi>t</mi> <mi>h</mi> <mi>e</mi> <mi>r</mi> <mi>w</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> <mo>.</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

式中，代表插值矩阵中的元素，代表插值矩阵中元素绝对值的平均值，κ为阀值参数。

5.如权利要求4所述的高效的结构频率响应拓扑优化方法，其特征在于：所述阀值参数κ的取值范围为1≤κ≤10。

6.如权利要求1-5任一项所述的高效的结构频率响应拓扑优化方法，其特征在于：所述设计变量为扩展系数；所述目标函数与所述约束函数关于设计变量的敏度计算是通过形状导数与伴随变量法进行的。

7.如权利要求1-5任一项所述的高效的结构频率响应拓扑优化方法，其特征在于：设计变量的更新是通过优化准则法进行的。

8.如权利要求1-5任一项所述的高效的结构频率响应拓扑优化方法，其特征在于：与时间相关的水平集函数由高斯径向基函数构建的插值矩阵与扩展系数向量构成的矩阵乘积形式为：

Aα＝Φ

式中，Φ为与时间相关的水平集函数对应的矩阵，α为仅与时间相关的扩展系数向量，A为仅与空间变量相关的高斯径向基函数构建的插值矩阵。