一种材料不确定性结构稳健性拓扑优化设计方法与流程

1.本发明涉及一种基于随机场模型的材料不确定性结构稳健性拓扑优化设计方法，属于结构设计领域。


背景技术：

2.近年来，随着计算机的高速发展以及有限元方法和优化理论的不断完善，结构拓扑优化设计迅速发展。但传统拓扑优化方法大多都建立在参数确定性假设上，忽略了实际工程中存在的材料性能、荷载以及几何不确定等因素对结构构型的影响，影响结构的安全性。因此考虑不确定性因素对结构拓扑优化结果影响受到了国内外学者的广泛关注，目前已有学者研究了结构材料性能、荷载及几何不确定性对拓扑优化结果的影响，主要是把不确定性参数当做随机变量来处理，且大多进行的是可靠性拓扑优化研究。
3.但实际工程中不确定性参数具有分散性，这些不确定性不仅随着时间的变化而波动，而且还依赖于空间位置的变化，因此引入随机场模型处理具有空间分散特征的不确定性参数(如材料性能)变得尤为重要，将不确定性变量简单地描述为随机变量可能会导致拓扑优化所得构型不一定是最优的，从而影响结构的安全性、经济性和适用性。而采用稳健性拓扑优化设计，能够让结构性能指标期望值达到预定目标值的同时降低性能指标对不确定性因素的敏感性，是结构设计领域的重要分支。


技术实现要素：

4.本发明的目的在于解决上述问题，提出一种基于随机场模型的材料不确定性结构稳健性拓扑优化设计方法，该方法将材料随机场模型与确定性拓扑优化模型结合，建立考虑弹性模量随机场的拓扑优化模型，再采用simp(solid isotropic microstructures with penalization)法进行稳健性拓扑优化。该方法结果准确、计算效率高，采用自编的fortran程序予以实现。
5.为实现上述目的，本发明采用技术方案如下：
6.步骤1，输入设计域尺寸、边界条件、荷载条件和体积分数，划分有限元网格并设置初始拓扑密度，建立确定性拓扑优化模型；
7.步骤2，对随机场进行离散，根据给定的弹性模量的概率分布特性，采用蒙特卡洛法生成m组弹模随机场数据；
8.步骤3，用步骤2中生成的m组弹性模量随机场数据一一替换确定性拓扑优化模型中对应单元的弹性模量，完成随机场在确定性拓扑优化模型中的实现，再进行确定性有限元分析得到柔顺度的均值与标准差，并建立稳健性拓扑优化目标；
9.步骤4，根据步骤3中建立的拓扑优化目标，采用simp法进行稳健性拓扑优化迭代计算，当拓扑优化计算收敛时输出优化结果。
10.作为本发明的进一步技术方案，步骤1中建立的确定性拓扑优化模型具体为：
[0011][0012]
式中，ne表示有限单元数量；ρ表示结构的拓扑密度；c
w
表示结构的最不利柔顺度；d(ρ)表示结构的整体位移矩阵；f表示结构的整体荷载矩阵；ρ
i
表示第i个单元的单元密度值；ρ
min
表示单元密度值变化下限；ρ
max
表示单元密度值变化上限；v0表示设计域的总体积；v*表示可用材料的总体积；表示体积分数；v
i
表示第i个单元的体积。
[0013]
作为本发明的进一步技术方案，步骤2中随机场的离散采用的方法为谱表示法，生成弹性模量随机场数据由matlab程序予以实现。
[0014]
作为本发明的进一步技术方案，步骤3具体为：
[0015]
(3
‑
1)读取步骤2中生成的m组弹性模量随机场数据，一一替换确定性拓扑优化模型中对应单元的弹性模量；
[0016]
(3
‑
2)进行有限元分析，分别计算m组弹性模量随机场下结构的柔顺度，再求出柔顺度的均值和标准差；
[0017]
(3
‑
3)建立稳健性拓扑优化目标函数，求出目标函数对拓扑密度的灵敏度表达式。
[0018]
作为本发明的进一步技术方案，步骤(3
‑
3)中，稳健性拓扑优化的目标函数为μ
f
(ρ) ωσ
f
(ρ)，其中ρ表示结构的拓扑密度，均值μ
f
(ρ)和标准差σ
f
(ρ)对拓扑密度的灵敏度表达式分别为：
[0019][0020][0021]
式中：ρ
i
表示第i个有限单元的单元密度值，与分别表示均值与标准差对ρ
i
的偏导；表示第h步迭代时第n组弹模随机场下的柔顺度对ρ
i
的偏导，n＝1,2,
…
,m。
[0022]
作为本发明的进一步技术方案，步骤4具体为：
[0023]
(4
‑
1)利用simp变密度法，构建拓扑密度与结构弹性模量之间的插值函数；
[0024]
(4
‑
2)根据步骤3中构建的稳健性拓扑优化目标，建立考虑弹性模量随机场的simp拓扑优化模型；
[0025]
(4
‑
3)采用灵敏度过滤技术对单元密度梯度进行修正，以解决拓扑优化过程中的数值不稳定问题；
[0026]
(4
‑
4)采用库恩塔克k
‑
t条件优化准则法对拓扑密度进行更新求解；
[0027]
(4
‑
5)迭代计算更新拓扑密度，当拓扑优化收敛后输出最终优化结果。
[0028]
作为本发明的进一步技术方案，步骤(4
‑
1)的插值函数为：
[0029]
e
i
＝e0(ρ
i
)
p
,i＝1,2,
…
,ne
[0030]
式中：ne表示有限单元数量；e
i
表示第i个单元的弹性模量；e0表示单元中材料完全保留时的单元弹性模量；ρ
i
表示第i个单元的单元密度值，p表示惩罚因子。
[0031]
作为本发明的进一步技术方案，步骤(4
‑
2)中构建的考虑弹性模量随机场的simp拓扑优化模型为：
[0032][0033]
式中ω为权系数，ω∈[0,1]；ξ为随机变量向量；ρ表示结构的拓扑密度；ρ
i
表示第i个单元的单元密度值；ρ
min
表示单元密度值变化下限；ρ
max
表示单元密度值变化上限；v0表示设计域的总体积；v*表示可用材料的总体积；v
i
表示第i个单元的体积；表示体积分数；k(ρ,ξ)与u(ρ,ξ)分别表示刚度矩阵与位移矩阵；μ
f
(ρ)与σ
f
(ρ)分别为柔顺度的均值和标准差。
[0034]
作为本发明的进一步技术方案，步骤(4
‑
3)中采用灵敏度过滤技术对单元密度梯度进行修正，第i个单元目标函数对单元密度梯度的修正公式为：
[0035][0036]
式中：表示修正后的第h步迭代时第i个单元密度梯度；c表示稳健性拓扑优化目标函数μ
f
(ρ) ωσ
f
(ρ)，其中μ
f
(ρ)与σ
f
(ρ)分别为柔顺度的均值和标准差；分别表示第h步迭代第i个单元和第j个单元的单元密度值；n
e
表示单元中心与第i个单元中心距离小于过滤半径r
min
的单元集合；h
ej
表示n
e
中第j个单元的权重，h
ej
＝max(0,r
min
‑
d(i,j))；r
min
表示过滤半径；d(i,j)表示第j个单元中心与第i个单元中心点间的距离。
[0037]
作为本发明的进一步技术方案，步骤(4
‑
4)中拓扑密度的计算表达式为：
[0038][0039]
式中：p为惩罚系数；k
i
表示第i个单元的单元刚度矩阵，d
i
为第i个单元的单元位移矩阵，k
0i
表示第i个单元为实体单元时的单元刚度矩阵。
[0040]
作为本发明的进一步技术方案，步骤(4
‑
5)中的收敛条件为：
[0041][0042]
式中：c
(h)
表示第h步迭代时的目标函数值；c
(h 1)
表示第h 1步迭代时的目标函数值；ε
t
为设定阈值。
[0043]
本发明采用以上技术方案与现有技术相比，具有以下技术效果：
[0044]
1.采用随机场来描述不确定性参数，能够处理不确定性参数具备空间分散特征的问题，比如材料空间分布不均、几何边界在空间中的扰动等，更符合实际工程情况；
[0045]
2.采用非侵入式随机有限元方法将不确定性参数随机场的模拟与确定性有限元分析分别独立进行，便于修改模型数据以及利用成熟的商业软件进行随机场模拟和拓扑优化计算，提高了计算效率；
[0046]
3.采用蒙特卡洛模拟来描述不确定性参数条件下结构的随机响应，方法清晰明了，计算精度较高。
附图说明
[0047]
图1为一种基于随机场模型的材料不确定性结构稳健性拓扑优化设计方法的流程图；
[0048]
图2为michell结构的设计域尺寸、边界条件和荷载条件；
[0049]
图3为michell结构网格划分图；
[0050]
图4为灵敏度过滤示意图；
[0051]
图5为michell结构的稳健性拓扑优化结果图。
具体实施方式
[0052]
下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。
[0053]
本技术领域技术人员可以理解的是，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。
[0054]
下面结合附图对本发明的技术方案做进一步的详细说明：
[0055]
一种基于随机场模型的材料不确定性结构稳健性拓扑优化设计方法，如图1所示，实现步骤如下：
[0056]
步骤1，在fortran程序中输入设计域尺寸、边界条件、荷载条件和体积分数，划分有限元网格，并设置初始拓扑密度等于体积分数，建立确定性拓扑优化模型；
[0057]
采用simp法建立的确定性拓扑优化数学模型具体为：
[0058][0059]
式中：ne表示有限单元数量；ρ表示结构的拓扑密度；c
w
表示结构的最不利柔顺度；d(ρ)表示结构的整体位移矩阵；f表示结构的整体荷载矩阵；ρ
i
表示第i个单元的单元密度值；ρ
min
表示单元密度值变化下限；ρ
max
表示单元密度值变化上限；v0表示设计域的总体积；v*表示可用材料的总体积；表示体积分数；v
i
表示第i个单元的体积。
[0060]
本发明实例中，michell结构设计域如图2所示，设计域尺寸为12m
×
6m，左下角点和右下角点为简支支承，结构上部承受均布荷载q＝10kn/m。有限元分析采用平面三角形单元，网格划分如图3所示，单元数为3600，节点数为1891，过滤半径取0.3m，体积分数取
[0061]
步骤2，采用谱表示法对随机场进行离散，根据给定的弹性模量的概率分布特性，采用蒙特卡洛法通过matlab程序模拟生成m组弹模随机场数据。
[0062]
本发明实例中，设计域内所有单元的弹性模量都是不确定的，服从对数正态分布，弹模均值为180gpa，变异系数为0.2，相关距离d
x
＝d
y
＝1m，泊松比v＝0.3，采用蒙特卡洛法生成400组随机场数据。
[0063]
步骤3，采用fortran程序，用步骤2中生成的m组弹性模量随机场中的参数一一替换确定性拓扑优化模型中对应单元的弹性模量，完成随机场在拓扑优化模型中的实现，再进行确定性有限元分析得到柔顺度均值和标准差，并建立稳健性拓扑优化目标。所述步骤3的具体步骤为：
[0064]
(3
‑
1)在fortran程序中读取步骤2中生成的m组弹性模量随机场数据，一一替换确定性拓扑优化模型中对应单元的弹性模量。
[0065]
(3
‑
2)进行有限元分析，分别计算m组弹性模量随机场下结构的柔顺度，再进行统计分析，求出柔顺度的均值μ和标准差σ：
[0066][0067][0068]
式中：h表示第h个迭代步；表示第h步时第m组弹模随机场下的柔顺度。
[0069]
本发明实例中，由抽样得到的400组随机场数据分别计算出其对应的柔顺度值，再计算400组样本对应柔顺度的均值与标准差。
[0070]
(3
‑
3)根据稳健性拓扑优化的目标函数μ
f
(ρ) ωσ
f
(ρ)，求出目标函数对拓扑密度的灵敏度表达式：
[0071]
[0072][0073]
式中：ρ
i
表示第i个有限单元的单元密度值，与分别为均值与标准差对ρ
i
的偏导；表示第h步迭代时第n组弹模随机场下的柔顺度对ρ
i
的偏导，n＝1,2,
…
,m。
[0074]
本发明实例中，权系数ω＝1，即目标函数为μ
f
(ρ) σ
f
(ρ)，目标函数相对拓扑密度的灵敏度表达式为
[0075]
步骤4，根据步骤3中建立的拓扑优化目标，采用simp法进行稳健性拓扑优化迭代计算，当拓扑优化计算收敛时输出优化结果。所述步骤4的具体步骤为：
[0076]
(4
‑
1)利用simp变密度法，构建拓扑密度与结构弹性模量之间的插值函数：
[0077]
e
i
＝e0(ρ
i
)
p
,i＝1,2,
…
,ne
[0078]
式中：ne表示有限单元数量；e
i
表示第i个单元的弹性模量；e0表示单元中材料完全保留时的单元弹性模量；ρ
i
表示第i个单元的单元密度值，p表示惩罚因子，一般取p＝3。
[0079]
(4
‑
2)根据步骤3中构建的稳健性拓扑优化目标，建立考虑弹性模量随机场的simp拓扑优化模型：
[0080][0081]
式中：ω为权系数，ω∈[0,1]。ξ为随机变量向量；ρ表示结构的拓扑密度；ρ
i
表示第i个单元的单元密度值；ρ
min
表示单元密度值变化下限；ρ
max
表示单元密度值变化上限；v0表示设计域的总体积；v*表示可用材料的总体积；v
i
表示第i个单元的体积；表示体积分数；k(ρ,ξ)与u(ρ,ξ)分别表示刚度矩阵与位移矩阵；μ
f
(ρ)与σ
f
(ρ)分别为柔顺度的均值和标准差。
[0082]
本发明实例中，取ρ
min
＝0.001，ρ
max
＝1；权系数ω＝1，即稳健性拓扑优化目标函数为μ
f
(ρ) σ
f
(ρ)。
[0083]
(4
‑
3)采用灵敏度过滤技术解决灰度单元、棋盘格和网格依赖性等数值不稳定问题，针对第i个单元目标函数对单元密度梯度的修正公式为：
[0084][0085]
式中：表示修正后的第h步迭代时第i个单元密度梯度；c表示稳健性拓扑优化目标函数μ
f
(ρ) ωσ
f
(ρ)，其中μ
f
(ρ)与σ
f
(ρ)分别为柔顺度的均值和标准差；分别表示第h步迭代第i个单元和第j个单元的单元密度值；n
e
表示单元中心与第i个单元中心距离小于过滤半径的单元集合，如图4所示，圆点表示各个单元的中心点，灰色单元为单元集合n
e
内的单元，其他单元为不属于单元集合n
e
的单元；h
ej
表示第j个单元的权重，其大小与i、j两单元中心点间的距离有关，可表示为：
[0086]
h
ej
＝max(0,r
min
‑
d(i,j))
[0087]
式中：r
min
表示过滤半径；d(i,j)表示第j个单元中心与第i个单元中心点间的距离。
[0088]
本发明实例中，过滤半径r
min
等于单元尺寸的1.5倍。
[0089]
(4
‑
4)利用simp法建立拓扑优化模型后，采用库恩塔克k
‑
t条件优化准则法对拓扑密度进行更新求解：
[0090][0091]
式中：p为惩罚系数；k
i
表示第i个单元的单元刚度矩阵，d
i
为第i个单元的单元位移矩阵，k
0i
表示第i个单元为实体单元时的单元刚度矩阵。
[0092]
(4
‑
5)迭代计算更新拓扑密度，当满足下式时表示整个拓扑优化达到了收敛，优化进程结束，输出最终优化结果。
[0093][0094]
式中：c
(h)
表示第h步迭代时的目标函数值；c
(h 1)
表示第h 1步的目标函数值；ε
t
为一足够小的实数。
[0095]
本发明实例中，取ε
t
＝0.001，得到优化结果如图5所示。
[0096]
本发明方案所公开的技术手段不仅限于上述技术手段所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。
[0097]
以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。