CN107491585B - 以随机位移响应方差为目标的结构拓扑优化设计方法 - Google Patents

Abstract

本发明公开了一种以随机位移响应方差为目标的结构拓扑优化设计方法，用于解决现有结构拓扑优化设计方法实用性差的技术问题。技术方案是采用大质量法将多点加速度转化为力施加到结构上，采用虚拟激励法将随机响应分析转换为简谐响应分析，并使用模态叠加法求解位移响应。然后在考虑频段内以结构关心位置处位移响应功率谱的方差值最小为目标，以结构体积分数为约束进行设计。相比背景技术的设计方法，本发明方法考虑以随机位移响应方差来直观衡量频响曲线的平稳性，并以此为目标进行拓扑优化设计。最终能够设计得到清晰有效的结构构型，从而能够满足工程实际中对结构随机位移响应曲线尽可能平稳的设计需求，具有极强的工程实际应用价值。

Description

以随机位移响应方差为目标的结构拓扑优化设计方法

技术领域

本发明涉及一种结构拓扑优化设计方法，特别涉及一种以随机位移响应方差为目标的结构拓扑优化设计方法。

背景技术

实际工程应用中，许多结构经常承受着幅值、频率和方向等不确定的复杂振动激励，特别是随机振动激励，需要以统计的方法对其进行研究。如航空航天飞行器在服役时受到的气动、噪声激励，能源动力装置在动力转换过程中承受的冲击声振激励，以及自然界中存在的风振、地震激励等。长期工作在随机振动环境下的结构，疲劳损伤甚至破坏是其失效的主要形式，给生命财产带来不可估量的损失。因此在结构设计时考虑其在随机振动环境下的性能表现就显得十分重要。

文献“Zhang W.H.,Liu H,Gao T.Topology optimization of large-scalestructures subjected to stationary random excitation An efficientoptimization procedure integrating pseudo excitation method and modeacceleration method[J].Computers&Structures,2015,158:61-70.”公开了一种随机力激励下的结构拓扑优化方法。该方法采用虚拟激励法与模态加叠加法相结合的方法实现了随机力激励下大规模自由度结构的高精度、高效率随机响应分析，并进行了高效的响应灵敏度分析，解决了传统的完全二次结合法应用于实际工程结构动力响应拓扑优化设计时存在的效率低下的问题，向随机响应拓扑优化的工程应用迈出了一大步。

发明专利201610398316.X中公开了一种基于大质量法的随机加速度激励下的结构拓扑优化设计方法，解决了现有随机激励下的结构拓扑优化设计方法无法实现多点加速度载荷施加的技术问题。

文献与发明专利中公开的方法虽然能够实现工程结构的随机力激励下的结构拓扑优化设计，发明专利中还能施加类似于地震激励的多点加速度随机载荷。但以上两种方法中优化设计目标均是位移响应均方根最小化，无法直观衡量结构响应曲线的平稳性，无法满足工程结构的实际设计需求。

发明内容

为了克服现有结构拓扑优化设计方法实用性差的不足，本发明提供一种以随机位移响应方差为目标的结构拓扑优化设计方法。该方法采用大质量法将多点加速度转化为力施加到结构上，采用虚拟激励法将随机响应分析转换为简谐响应分析，并使用模态叠加法求解位移响应。然后在考虑频段内以结构关心位置处位移响应功率谱的方差值最小为目标，以结构体积分数为约束进行设计。相比背景技术的设计方法，本发明方法考虑以随机位移响应方差来直观衡量频响曲线的平稳性，并以此为目标进行拓扑优化设计。最终能够设计得到清晰有效的结构构型，从而能够满足工程实际中对结构随机位移响应曲线尽可能平稳的设计需求，具有极强的工程实际应用价值。

本发明解决其技术问题所采用的技术方案：一种以随机位移响应方差为目标的结构拓扑优化设计方法，其特点是包括以下步骤：

步骤一、对结构初始几何模型进行有限元网格划分，获得有限元模型。在拟施加激励位置外建一大质量点，大质量点与结构上承受加速度激励的节点之间通过刚性单元连接。

步骤二、在大质量点处施加与加速度等效的力载荷。随机载荷采用限带白噪声激励bp(t)，即在整个激励频段上拥有完全相同的功率谱密度值。其中，p(t)为d维随机激励向量，其功率谱密度矩阵为S_p(ω)，d为载荷中力的个数，ω为激励圆频率。b为n×d的转换矩阵，用于将d维随机激励向量转换为n维激励分布向量，n为包含大质量节点自由度的结构总自由度数。载荷激励频段为[ω_a,ω_b]，ω_b、ω_a分别表示激励圆频率上下限。由于功率谱密度矩阵S_p(ω)为Hermitian矩阵，因此存在下式分解

其中，Q为功率谱密度矩阵S_p(ω)的秩，γ_q为d维列向量，表示第q个虚拟简谐激励，(γ_q)^*为其共轭矩阵。1≤q≤Q，上标T表示向量或矩阵的转置。

采用大质量法实现加速度载荷的施加，此时的动力学平衡方程按分块矩阵形式表示为

其中，M_ss、M_sb、M_bs、M_bb为结构整体质量矩阵M按结构基础节点和自由节点分块后得到的矩阵，其中下标s表示自由节点自由度，下标b表示基础节点自由度。同理C_ss、C_sb、C_bs、C_bb为结构整体阻尼矩阵C分块后得到的矩阵，K_ss、K_sb、K_bs、K_bb为结构整体刚度矩阵K分块后得到的矩阵，同时M_L为对应于基础节点质量矩阵M_bb的大质量矩阵。

为加速度幅值向量x的分块形式，同理

为速度幅值向量

的分块形式，x_s、x_b为位移幅值向量x的分块形式。

为拟施加的基础加速度载荷向量。

将上述公式中的第二个项展开

上式两侧左乘M_L的逆矩阵

当质量点的质量很大时，

中对角元素趋于零，则基础激励处实际获得的加速度为：

步骤三、设置拓扑设计变量η_h初始值，优化迭代时其值在0-1之间变化，h是正整数，表示设计域单元编号。给定实体材料杨氏模量E、密度ρ和泊松比μ。每次迭代后，根据当前设计变量值，更新结构有限元模型中的相应材料属性。为了降低局部模态对拓扑优化过程带来的负面影响，选取多项式插值对单元材料属性进行匹配，分别更新每一个有限元单元在当前迭代步下的杨氏模量E_h和密度ρ_h。

ρ_h＝η_hρ (6)

步骤四、从模态分析中提取每个单元的刚度矩阵K_h、质量矩阵M_h，以及结构的前l阶固有频率ω_i和模态振型矩阵

1≤i≤l。

由各阶模态振型向量

组成。采用瑞利阻尼时，结构的前l阶模态的阻尼比ζ_i按下式进行计算：

式中，α与β为瑞利阻尼系数。

采用虚拟激励法计算结构指定自由度r的随机位移响应功率谱密度

式中，u表示位移，(g_q(t))_r为结构指定自由度r在第q个虚拟简谐激励γ_q下的位移响应，||(g_q(t))_r||表示复数(g_q(t))_r的模，即位移响应幅值。(g_q(t))_r采用模态叠加的方法进行计算，其计算公式为

由于采用大质量法进行基础加速度激励时，大质量点激励方向的自由度是放开的，模态分析时会存在刚体模态。在进行模态固有频率和振型提取时需过滤掉刚体模态信息，模态叠加后所得位移为结构指定自由度r相对于基础大质量点的位移响应。式中a为n维列向量，除第r项为1外其余项均为0，T表示向量转置。

为各阶模态振型向量，1≤i≤l，l为进行模态叠加时提取的模态阶数。b为n×d转换矩阵，用于将d维随机激励向量转换为n维激励分布向量，n为包含大质量节点的结构总自由度数。假如p(t)中第k个力施加在第z个自由度上，则b的第k列第z行的元素值是1，k列中其它元素值均为0，1≤k≤d，1≤z≤n。e^jωt表示以自然常数e为底数的指数函数，ω为激励圆频率，j²＝-1。H_i为质量矩阵归一化进行解耦后的第i个单自由度系统的频响函数，计算公式为

其中，ω为激励圆频率，ω_i和ζ_i为第i阶固有频率和模态阻尼比，1≤i≤l，l为进行模态叠加时提取的模态阶数。

步骤五、按固有频率对考虑频段进行细分，得到若干个频率采样点。这些频率采样点的位移响应功率谱的方差值

按下式计算

其中，N为频段[ω_a，ω_b]内采样点总数，ω_b、ω_a分别表示激励圆频率上下限，ω_ξ为第ξ个采样点圆频率，1≤ξ≤N。

为指定自由度r处所有频率采样点位移响应功率谱的平均值，其计算式为：

显然，位移响应方差值

总是正值，其值越小表示指定自由度上的频响曲线越平稳。

步骤六、给定体积分数约束上限V_fU，采用下式计算当前迭代步的结构体积分数V_f。

其中，V_h表示第h个实体材料单元的体积，V₀为初始实体结构体积。

步骤七、以随机位移响应方差最小为目标的结构拓扑优化模型如下

式中，η代表设计变量的集合，η_h为第h个有限单元的设计变量值，其中1≤h≤enum，enum代表设计域单元总数。迭代过程中设计变量η_h值在0-1之间变化，表示实体材料的有无。为避免有限元分析时结构刚度矩阵的奇异，引入设计变量下限η_L。优化目标为结构指定自由度r处随机位移响应方差

最小，约束条件为体积约束。

步骤八、求得目标函数和约束条件关于设计变量的灵敏度。选取基于梯度的GCMMA优化算法，采用自编拓扑优化程序CommonOpt.exe进行优化迭代，得到以随机位移位移响应方差为目标的结构拓扑优化设计结果。

本发明的有益效果是：该方法采用大质量法将多点加速度转化为力施加到结构上，采用虚拟激励法将随机响应分析转换为简谐响应分析，并使用模态叠加法求解位移响应。然后在考虑频段内以结构关心位置处位移响应功率谱的方差值最小为目标，以结构体积分数为约束进行设计。相比背景技术的设计方法，本发明方法考虑以随机位移响应方差来直观衡量频响曲线的平稳性，并以此为目标进行拓扑优化设计。最终能够设计得到清晰有效的结构构型，从而能够满足工程实际中对结构随机位移响应曲线尽可能平稳的设计需求，具有极强的工程实际应用价值。

下面结合附图和具体实施方式对本发明作详细说明。

附图说明

图1是本发明以随机位移响应方差为目标的结构拓扑优化设计方法的流程图。

图2是本发明方法实施例中初始几何模型及载荷施加示意图。

图3是本发明方法实施例设计结果图。

具体实施方式

参照图1-3。本发明以随机位移响应方差为目标的结构拓扑优化设计方法具体步骤如下：

步骤一、建立有限元模型。采用ANSYS的参数化语言APDL编程，建立长宽厚分别为1m、0.48m、0.005m的矩形平板结构，将其划分为100×48的四边形网格，选择单元类型为plane42平面单元。节点编号从1-4949，共4949个节点9898个自由度。给定材料杨氏模量E＝206GPa，密度ρ＝7850Kg/m³,泊松比μ＝0.3。初始结构质量为18.84Kg，在拟施加激励位置外建立一大质量节点，编号为4950，在该节点上添加质量单元Mass21，施加质量大小为19×10⁸Kg。在关心响应的节点102处添加负载0.5Kg。大质量点与结构上承受加速度激励的节点之间通过MPC184多点约束单元连接。

步骤二、设置激励载荷。本方法中只在大质量点竖直方向施加加速度激励，大小为5000m/s²。因此通过大质量法转换的随机力激励向量p(t)的维数d＝1，对应的功率谱密度矩阵S_p(ω)为1维矩阵，其值为(5000×19×10⁸)²N²/(rad/s)。b为n×1的转换矩阵，将1维随机激励向量转换为整个结构的激励分布向量，n＝10000为包含大质量点的整个结构总自由度数。ω为激励圆频率，本方法中载荷激励频段为[0,3140]rad/s。

由于功率谱密度矩阵S_p(ω)为Hermitian矩阵，因此存在下式分解

其中Q＝1为功率谱密度矩阵S_p(ω)的秩，γ₁＝5000×19×10⁸N表示第1个虚拟简谐激励。

采用大质量法实现加速度载荷的施加时，结构动力学平衡方程按分块矩阵形式可表示为

其中M_ss、M_sb、M_bs、M_bb为结构整体质量矩阵M按结构基础节点和自由节点分块后得到的矩阵，其中下标s表示自由节点自由度，下标b表示基础节点自由度。同理C_ss、C_sb、C_bs、C_bb为结构整体阻尼矩阵C分块后得到的矩阵，K_ss、K_sb、K_bs、K_bb为结构整体刚度矩阵K分块后得到的矩阵，同时M_L为对应于基础节点质量矩阵M_bb的大质量矩阵。

为加速度幅值向量

的分块形式，同理

为速度幅值向量

的分块形式，x_s、x_b为位移幅值向量x的分块形式。

为拟施加的基础加速度载荷向量。

将上述公式中的第二行展开

上式两侧左乘M_L的逆矩阵

当质量点的质量很大时，

中对角元素趋于零，则

基础激励处实际获得的加速度为：

这样就可以将加速度载荷等效为虚拟简谐力载荷施加到结构上。

步骤三、设置拓扑设计变量η_h初始值为0.5，优化迭代时其值在0-1之间变化，h是正整数，表示设计域单元编号。本实施例中所有平面单元均为设计单元，总数enum＝4800。给定实体材料杨氏模量E＝206GPa、密度ρ＝7850Kg/m³、泊松比μ＝0.3。每次迭代后，根据当前设计变量值，更新结构有限元模型中的相应材料属性。在本方法中，为了降低局部模态对拓扑优化过程的负面影响，选取多项式插值对单元材料属性进行匹配，分别更新每一个有限单元在当前迭代步下的杨氏模量E_h和ρ_h。

ρ_h＝η_hρ (6)

步骤四、响应分析计算。从模态分析中提取每个单元的刚度矩阵K_h和质量矩阵M_h，同时提取结构的前l＝30阶固有频率ω_i和模态振型矩阵

提取时不考虑大质量点及刚性单元的振型，

为9898行30列的矩阵，由各阶模态振型向量

组合而成，其中1≤i≤l。采用瑞利阻尼时，结构的前l＝30阶模态的阻尼比ζ_i可按下式进行计算：

式中α＝0.01与β＝0.00001，为瑞利阻尼系数。

采用虚拟激励法计算结构指定自由度r＝204的随机位移响应功率谱密度

式中u表示位移，(g_q(t))₂₀₄为结构指定自由度r＝204在第q个虚拟简谐激励γ_q下的位移响应，||(g_q(t))₂₀₄||表示复数(g_q(t))₂₀₄的模，即位移响应幅值。(g_q(t))₂₀₄采用模态叠加的方法进行计算，其计算公式为

由于采用大质量法进行基础加速度激励时，大质量点激励方向的自由度是放开的，模态分析时会存在刚体模态。在进行模态固有频率和振型提取时过滤掉刚体模态信息，模态叠加后所得位移为结构指定自由度r＝204相对基础大质量点的位移响应。式中a为n＝10000维列向量，除第r＝204项为1外其余项均为0，T表示向量转置。b为10000×1转换矩阵，用于将1维随机激励向量转换为10000维激励分布向量。p(t)中只有1个力施加在结构外大质量节点的第10000个自由度上，则b中只有第1列第10000行的元素值是1，其它元素值均为0。e^jωt表示以自然常数e为底数的指数函数，ω为激励圆频率，j²＝-1。H_i为质量矩阵归一化进行解耦后的第i个单自由度系统的频响函数，计算公式为

步骤五、考虑频段内位移响应方差计算。先按固有频率对考虑频段进行细分，再在两个相邻固有频率区间按固有频率细分点多的原则进行进一步细分，得到若干个频率采样点。这些频率采样点在结构指定自由度r＝204上的位移响应功率谱的方差值

可按下式计算

其中N为考虑频段[0,3140]rad/s内采样点总数，优化迭代过程中N是变化的，ω_ξ为第ξ个采样点圆频率，1≤ξ≤N。

为所有频率采样点在结构指定自由度r＝204上的位移响应功率谱的平均值，其计算式为：

很显然，位移响应方差值

总是正值，其值越小表示指定自由度上的频响曲线越平稳。

步骤六、给定体积分数约束上限V_f＝50％，采用以下计算式计算当前迭代步的结构体积分数V_f。

其中V_h表示第h个实体材料单元的体积，V₀为初始实体结构体积。

步骤七、定义拓扑优化模型。

本方法中采用的以随机位移响应方差最小为目标的结构拓扑优化模型如下

式中η代表设计变量的集合，η_h为第h个有限单元的设计变量值，其中1≤h≤enum，enum代表设计域单元总数4800。初始迭代步设计变量值η_h均为0.5，迭代过程中其值在0-1之间变化，表示实体材料的有无。为避免有限元分析时结构刚度矩阵的奇异，引入设计变量下限η_L＝0.001。优化目标为结构指定自由度r＝204处随机位移响应方差

最小，约束条件为体积约束，体积分数约束上限为50％。另外优化时对设计变量η_h做了对称处理，保证结构的对称性。

步骤八、进行灵敏度分析与优化迭代。

求得目标函数和约束条件关于设计变量的灵敏度。选取基于梯度的GCMMA优化算法，采用自编拓扑优化程序CommonOpt.exe进行优化迭代，得到以随机响应方差为目标的结构拓扑优化设计结果。特别说明，其他的基于梯度的优化算法如ConLin，GCM，MDPA，SLP，QP等均能实现优化迭代。此外，一些其他的优化方法如优化准则法、数学规划法、渐进结构优化法等也均能进行本发明方法的优化设计。

本实施例经过235步迭代后得到优化设计结果。考虑激励频段[0,3140]rad/s，初始结构指定自由度r＝204的随机位移响应方差值为0.9475m⁴，优化设计得到的结构满足体积约束的前提下指定自由度r＝204的随机位移响应方差值为5.72e-8m⁴，随机位移响应方差值降幅达到99.99％。位移响应功率谱幅值从4.19e^-2m²降至4.71e^-6m²，响应曲线峰值大幅下降，响应曲线也更为平稳。最终设计得到的结构构型清晰有效，考虑了实际工程中对结构响应平稳性的需求，给工程设计提供了思路。

Claims (1)

1.一种以随机位移响应方差为目标的结构拓扑优化设计方法，其特征在于包括以下步骤：

步骤一、对结构初始几何模型进行有限元网格划分，获得有限元模型；在拟施加激励位置外建一大质量点，大质量点与结构上承受加速度激励的节点之间通过刚性单元连接；

步骤二、在大质量点处施加与加速度等效的力载荷；随机载荷采用限带白噪声激励bp(t)，即在整个激励频段上拥有完全相同的功率谱密度值；其中，p(t)为d维随机激励向量，其功率谱密度矩阵为S_p(ω)，d为载荷中力的个数，ω为激励圆频率；b为n×d的转换矩阵，用于将d维随机激励向量转换为n维激励分布向量，n为包含大质量节点自由度的结构总自由度数；载荷激励频段为[ω_a,ω_b]，ω_b、ω_a分别表示激励圆频率上下限；由于功率谱密度矩阵S_p(ω)为Hermitian矩阵，因此存在下式分解

其中，Q为功率谱密度矩阵S_p(ω)的秩，γ_q为d维列向量，表示第q个虚拟简谐激励，(γ_q)^*为其共轭矩阵；1≤q≤Q，上标T表示向量或矩阵的转置；

采用大质量法实现加速度载荷的施加，此时的动力学平衡方程按分块矩阵形式表示为

其中，M_ss、M_sb、M_bs、M_bb为结构整体质量矩阵M按结构基础节点和自由节点分块后得到的矩阵，其中下标s表示自由节点自由度，下标b表示基础节点自由度；同理C_ss、C_sb、C_bs、C_bb为结构整体阻尼矩阵C分块后得到的矩阵，K_ss、K_sb、K_bs、K_bb为结构整体刚度矩阵K分块后得到的矩阵，同时M_L为对应于基础节点质量矩阵M_bb的大质量矩阵；

为加速度幅值向量

的分块形式，同理

为速度幅值向量

的分块形式，x_s、x_b为位移幅值向量x的分块形式；

为拟施加的基础加速度载荷向量；

将上述公式中的第二行展开

上式两侧左乘M_L的逆矩阵

当质量点的质量很大时，

中对角元素趋于零，则基础激励处实际获得的加速度为：

步骤三、设置拓扑设计变量η_h初始值，优化迭代时其值在0-1之间变化，h是正整数，表示设计域单元编号；给定实体材料杨氏模量E和密度ρ；每次迭代后，根据当前设计变量值，更新结构有限元模型中的相应材料属性；为了降低局部模态对拓扑优化过程带来的负面影响，选取多项式插值对单元材料属性进行匹配，分别更新每一个有限元单元在当前迭代步下的杨氏模量E_h和密度ρ_h；

ρ_h＝η_hρ (6)

步骤四、从模态分析中提取每个单元的刚度矩阵K_h、质量矩阵M_h，以及结构的前l阶固有频率ω_i和模态振型矩阵

1≤i≤l；

由各阶模态振型向量

组成；采用瑞利阻尼时，结构的前l阶模态的阻尼比ζ_i按下式进行计算：

式中，α与β为瑞利阻尼系数；

采用虚拟激励法计算结构指定自由度r的随机位移响应功率谱密度

式中，u表示位移，(g_q(t))_r为结构指定自由度r在第q个虚拟简谐激励γ_q下的位移响应，||(g_q(t))_r||表示复数(g_q(t))_r的模，即位移响应幅值；(g_q(t))_r采用模态叠加的方法进行计算，其计算公式为

由于采用大质量法进行基础加速度激励时，大质量点激励方向的自由度是放开的，模态分析时会存在刚体模态；在进行模态固有频率和振型提取时需过滤掉刚体模态信息，模态叠加后所得位移为结构指定自由度r相对于基础大质量点的位移响应；式中a为n维列向量，除第r项为1外其余项均为0，T表示向量转置；

为各阶模态振型向量，1≤i≤l，l为进行模态叠加时提取的模态阶数；b为n×d转换矩阵，用于将d维随机激励向量转换为n维激励分布向量，n为包含大质量节点的结构总自由度数；假如p(t)中第k个力施加在第z个自由度上，则b的第k列第z行的元素值是1，k列中其它元素值均为0，1≤k≤d，1≤z≤n；e^jωt表示以自然常数e为底数的指数函数，ω为激励圆频率，j²＝-1；H_i为质量矩阵归一化进行解耦后的第i个单自由度系统的频响函数，计算公式为

其中，ω为激励圆频率，ω_i和ζ_i为第i阶固有频率和模态阻尼比，1≤i≤l，l为进行模态叠加时提取的模态阶数；

步骤五、按固有频率对考虑频段进行细分，得到若干个频率采样点；这些频率采样点的位移响应功率谱的方差值

按下式计算

其中，N为频段[ω_a，ω_b]内采样点总数，ω_b、ω_a分别表示激励圆频率上下限，ω_ξ为第ξ个采样点圆频率，1≤ξ≤N；

为指定自由度r处所有频率采样点位移响应功率谱的平均值，其计算式为：

显然，位移响应方差值

总是正值，其值越小表示指定自由度上的频响曲线越平稳；

步骤六、给定体积分数约束上限V_fU，采用下式计算当前迭代步的结构体积分数V_f；

其中，V_h表示第h个实体材料单元的体积，V₀为初始实体结构体积；

步骤七、以随机位移响应方差最小为目标的结构拓扑优化模型如下

式中，η代表设计变量的集合，η_h为第h个有限单元的设计变量值，其中1≤h≤enum，enum代表设计域单元总数；迭代过程中设计变量η_h值在0-1之间变化，表示实体材料的有无；为避免有限元分析时结构刚度矩阵的奇异，引入设计变量下限η_L；优化目标为结构指定自由度r处随机位移响应方差

最小，约束条件为体积约束；

步骤八、求得目标函数和约束条件关于设计变量的灵敏度；选取基于梯度的GCMMA优化算法，采用自编拓扑优化程序CommonOpt.exe进行优化迭代，得到以随机位移位移响应方差为目标的结构拓扑优化设计结果。

Images