求最小特征频率问题的分片常数水平集方法_高祥静.pdf

第 卷第期杭 州 电 子 科 技 大 学 学 报（自然科学版）年月 （）：求最小特征频率问题的分片常数水平集方法高祥静，张郑芳，董燕强（杭州电子科技大学理学院，浙江 杭州 ）收稿日期：基金项目：浙江省自然科学基金资助项目（）作者简介：高祥静（），女，研究方向：偏微分方程数值解。：。通信作者：张郑芳，副教授，研究方向：偏微分方程数值解。：。摘要：通过引入虚拟材料，将线性弹性系统扩展到一个固定的背景域，研究结构振动最小特征频率的最大化问题。用分片常数水平集（，）方法表示不同材料区域，运用灵敏度分析得到 函数最小特征频率的导数。提出一种重新启动的共轭梯度惩罚算法对二维及三维悬臂梁进行数值实验，实验结果表明，方法可以在迭代过程中自动成核和消除空洞，验证了共轭梯度惩罚算法的有效性。关键词：特征频率；分片常数水平集方法；灵敏度分析；共轭梯度惩罚算法中图分类号：文献标志码：文章编号：（）引言结构形状及拓扑优化是解决结构设计问题效率最大化的重要工具。水平集方法是解决这类问题的传统方法，但在二维情况下，水平集方法很难生成核。分片常数水平集（，）方法是应用于形状和拓扑优化中的新方法。等引入 求解了一个以应变能为目标函数的结构形状和拓扑优化问题。等使用变分二元水平集方法处理应变能极小化问题。等将 方法推广到特征值优化问题的边界控制中，用于区分 边界和 边界，随后，提出一种基于 的线性弹性系统最小特征值最大化方法，并将 方法推广到离散形式上，用于处理谐波激励下阻尼层的拓扑优化问题。对于 方法，最速下降法在迭代开始时速度很快，但收敛速度很慢。共轭梯度法可以在有限的迭代次数内停止，但算法的收敛速度是线性的。本文研究线性弹性系统最小特征频率问题，提出一种重新启动的共轭梯度惩罚算法。特征频率问题假设线性弹性系统特征值如下：（）内上（）上（）式中，是（，）的光滑有界的连通子集，是外法线向量，是 边界，是 边界，它们的关系是，。表示胡克定律且有（），其中和是材料的拉梅模，是对称矩阵，表示单位矩阵。在平面应力问题中，（）（）且（），其中是杨氏模量，是泊松比，是材料的密度。（）.（）表示应变张量。（）表示特征振动频率，是相关的特征模。对式（）使用分部积分，得到：（）:（），（）其中，（）（）:，。最小特征频率表示为：（）（），（）:（）本文的优化目的是使最小特征频率最大化，因此考虑目标函数的极小化问题（）（）。设是中包括所有可容许形状的有界开集。定义固定体积的可容许形状为:。优化问题定义为：（）（）本文采用虚拟材料方法来求问题（）的最优解，将式（）的弹性方程中扩展到整个区域，并在 中填充一种虚拟材料来模拟空洞。拉梅模和密度函数分别为：内 内，内 内，内 内其中，和表示在 内虚拟材料的拉梅模，表示在 内虚拟材料的密度。因此，区域的弹性方程重新定义为：（）内 上（）上（）式中，胡克定律中和替换原来中的和。区域被分成部分：是 边界，是 边界。进一步限制可容许形状，满足，即边界是固定边界 的子集，边界可以在 之外的任何地方。对式（）使用分部积分，得到弱形式（）:（），（），其中，（）（）:，。研究如下优化问题：（）（）（）式中，（），（）:（）（）:（），（）是弹性方程（）的最小特征频率，是与（）相关的特征模。分片常数水平集方法和灵敏度分析 分片常数水平集方法采用 函数（）处理问题（），根据 函数，拉梅模和密度函数重新定义如下：（）（）（）（）（）（）（）第期高祥静，等：求最小特征频率问题的分片常数水平集方法运用 函数对线性弹性系统方程进行改进，得到：（）（）（）（）（）内（）上（）（）上（）由式（）可以得到如下弱形式：（）（）:（）（）（）（），（）（）根据 函数，的面积（体积）可以表示为（）。函数需要保存个常量中的个，因此引入新的约束（）（）（）。接下来，将依赖于的优化问题（）改为依赖于 函数的优化问题，即：（）（）（）使得（）（）（）:（）（）（），（），（），其中。灵敏度分析为了求解优化问题（），采用连续灵敏度分析推导最小特征频率（）相对于 函数的导数。首先求最小特征频率（）在方向上的一阶变分，（）:（）:（）（）:（）（）（）（）为了简化符号，这里及下面省略了变量。（）（）关于的一阶变分为（：（）（）（）（）（）此外，下式成立：（）:（）（）:（）（）在式（）中取，式（）中取 和，再结合式（），简化式（），可得：（）（）:（）（）:（）（）:（）（）（）（）（）:（）（）（）（）（）:（）对式（）中，和分别在方向上求导，且（），得到 关于在方向上的一阶变分，（）（）（）（）:（）（）杭州电子科技大学学报（自然科学版）年所以，关于 的导数为：（）（）（）（）:（）（）共轭梯度惩罚算法针对优化问题（），本文提出一种重新启动的共轭梯度惩罚算法，算法中的步骤借鉴 和 格式的共轭梯度算法。在目标函数中加入二次惩罚项和正则项，得到二次惩罚目标函数为：（；，）（）（）用投影限制（）中的，得到，.。因此，（；，）关于求导为：（；，）（）（）（）（）:（）（）（）（）（）在共轭过程中，下降方向（）；，），（）；，），其中迭代次数，（）；，）（）（）；，）（）。重新启动过程时，只需要将负梯度作为下降方向，即（）；，）。在第次迭代后更新为（）（），其中（）（）；，）。当特征频率对应多重时，需要选择合适的来计算（；，）的导数。受文献，的启发，本文选择接近上一次迭代中使用的，即（）（）（）其中，对于一个给定的，是与（）相关的特征模，（）:，:（）（）。重新启动的共轭梯度惩罚算法主要步骤如下：初始化：整个区域的（），已知参数，和，迭代次数，。（）采用有限元方法运用式（）计算得到最小特征频率（）以及对应的（）。如果（）（）是多重的，则从式（）中选择（）。（）由式（）以及（），（），（），求解（），）。（）由共轭梯度算法得到下降方向（），），其中（）；，）（）（）；，）（）。（）通过迭代来更新（）。本步骤中的迭代次数设为，且 。（）为了在数值实验中平衡 和惩罚项，假设惩罚项满足 （），（，设步长 （），），（，），计算（）；（）根据投影定义对作用投影，得到，；（）如果，（），设定（），然后进行下一步，否则，（）（），回到本步骤第一步，计算。（）如果满足停止条件，则输出最小特征频率。否则，取，回到步骤。第期高祥静，等：求最小特征频率问题的分片常数水平集方法数值实验结果及分析实验平台的操作系统为 ，测试工具为 。数值实验中，原始材料的杨氏模归一化为，泊松比为.，原始材料密度。假设虚拟材料有相同的泊松比，较小的杨氏模为 ，较小的密度为 ，固定面积。采用有限元方法，使用 中的 函数将线性弹性系统离散化为广义矩阵特征值问题。在 中运行本文提出的重新启动的共轭梯度惩罚算法，验证算法在二维矩形和三维长方体内解决最小特征频率问题的能力。二维矩形悬臂的二维矩形悬臂如图（）所示，左边是固定的零 边界，其余边界是自由的零 边界。在该区域右侧的中间填充较重的材料，其 ，假设该小方形区域在迭代过程中无需进行优化。将整个区域划分为 的网格并令 ，.。初始 函数在除了右边中间小正方形外的整个区域设定为（）（）。原始材料区域用黑色着色，虚拟材料用白色着色。，次迭代后，介质悬臂的演变情况如图（）图（）所示。图二维矩形悬臂的演化剖面由图可以看出，原始材料区域面积逐渐增加并趋于稳定，迭代开始时没有空洞，随后出现空洞，随着迭代次数的增加，二维矩形悬臂自动生成核并消除空洞。图显示了最小特征频率和原始材料区域面积的迭代历程。图最小特征频率 和原始材料区域面积的迭代历程由图可以看出，在前 次迭代中，最小特征频率和原始材料区域的面积都在增加；迭代到 次杭州电子科技大学学报（自然科学版）年后，条曲线趋于稳定，但仍存在一定的震荡，因为在一个单位元中，原始材料和虚拟材料之间的变化引起最小特征频率小幅变化，引起轻微震荡。三维长方体悬臂的三维长方体悬臂如图（），左边是固定的零 边界，其余边界是自由的零 边界。在该区域右侧的中间填充较重的材料，其 ，假设小立方体在迭代过程中无需进行优化。将整个曲域划分为 的小立方体并令 ，.。初始 函数在除了右边中间小立方外的整个区域设定为（）（）。黑色区域表示原材料，白色区域表示虚拟材料。，次迭代后，介质悬臂的演变情况如图（）图（）所示。图三维长方体悬臂的演化剖面由图可以看出，原始材料区域体积逐渐增加并趋于稳定。在迭代开始时没有空洞，随后出现空洞，随着迭代次数的增加，三维长方体悬臂自动生成核并消除空洞。图显示了最小特征频率 和原始材料区域体积的迭代历程。图最小特征频率 和原始材料区域体积的迭代历程由图可以看出，在前 次迭代过程中，特征频率 和原始材料区域的体积都在增加，之后趋于第期高祥静，等：求最小特征频率问题的分片常数水平集方法稳定。综上分析可以看出，用 方法表示不同材料区域，迭代过程中，本文提出的重新启动的共轭梯度惩罚算法能自动生成核并消除空洞，验证了算法的有效性。结束语本文研究了面积和体积约束下最小特征频率问题，提出一种基于 方法的重新启动的共轭梯度惩罚算法。迭代过程中，本文算法可自动生成核，并消除空洞，解决了特征频率多重的问题。在固定背景域内，采用虚拟材料求解线性弹性系统时，原始材料的分散并没有阻止迭代的进行。但是，本文没有考虑改变二维矩形悬臂和三维长方体悬臂的尺寸对最小特征频率的影响，后续将针对这个问题展开研究。参考文献 ，：，（）：，：，（）：，（）：，（）：，：，：，（）：，（）：，（，）：，（），：；杭州电子科技大学学报（自然科学版）年