本发明涉及结构拓扑优化设计技术领域，特别涉及一种力热耦合环境下层合板非概率可靠性拓扑优化方法，该方法考虑了力热耦合环境下材料弹性参数、热属性参数、外部结构力热载荷等参数的不确定性，建立了非概率可靠性指标，并基于simp插值模型对复合材料层合结构进行拓扑优化。
背景技术：
展望21世纪的航空航天蓝图，大幅减重仍是飞行器设计的迫切需求。碳纤维增强常刚度层合板是复合材料单层以特定的铺层方式层叠而成的重要结构复合材料，因其重量轻、制造技术成熟而被广泛应用。随着飞服役环境的日益严酷，飞行器结构精细化设计逐渐由单一的承载功能向承载、隔热、减振等多功能过渡，以高导热碳纤维为主要增强体的新型导热复合材料，可以有效地克服金属导热材料的缺点，实现更高的热导率，显示出广阔的应用前景。现有复合材料设计的研究主要集中于尺寸、几何层级，而复合材料层合板结构的拓扑优化设计还处于起步阶段，因此研究形成一套考虑承载防热多功能一体化层合板结构可靠性拓扑优化设计方法对改善复合材料层合结构的性能具有重要的理论意义。
技术实现要素：
本发明要解决的技术问题是：针对复合材料层合结构，克服现有技术的不足，提供一种针对复合材料结构在力热耦合环境下的非概率可靠性拓扑优化设计方法。本发明充分考虑实际工程中复杂力热环境的影响，构建考虑各类不确定性因素的非概率可靠性分析指标，所得到的结果更加符合真实情况，工程适用性更强。
本发明采用的技术方案为：一种力热耦合环境下层合板非概率可靠性拓扑优化方法，该方法实现步骤如下：
步骤一：根据给定的复合材料单向板材料参数和铺层方式，建立层合板力热耦合本构模型，利用一阶层合理论，计算层合板等效热导率，通过下式定义：
其中△q表示热量，△t表示时间，
mc,nn＝∑ωθ(mc,pcos2θ mc,tsin2θ)
式中，mc,p表示复合材料单向板沿着纤维方向的热导率，mc,t表示复合材料单向板垂直于纤维方向的热导率。根据复合材料经典层合理论，温度沿着层合板厚度方向线性变化，即：
t(x,y,z)＝t0(x,y) tz'(x,y)·z
其中t0(x,y)表示层合结构中面温度t′z(x,y)表示沿层合板厚度z方向的温度梯度。
以对称铺层的层合板为例，全局坐标系下第k层的应力通过下式得到：
其中σx,σy和τxy表示全局坐标系下的应力分量，εx，εy和γxy表示全局坐标系下的应变分量，αx，αy和αxy表示全局坐标系下的热膨胀系数。
式中e1，e2分别是单向板主方向1的杨氏模量、主方向2的杨氏模量，g12是主方向1和2的剪切模量，μ21和μ12分别是主方向1和2的泊松比；tk表示层合板第k层的转换矩阵，即：
式中，θk表示第k层的纤维角度。
复合材料层合板的单元弹性矩阵
其中，tk为层合板第k层的厚度，h为层合板的总厚度。在结构有限元法中，单元刚度矩阵ki可通过下式得到：
其中，b为单元应变-位移矩阵。复合材料层合板的总体刚度矩阵可通过单元刚度矩阵求和得到：
步骤二：根据实际工程中服役环境和层合板结构材料属性，选取描述层合板结构特性或外部力热载荷的不确定性参数b，根据有限元平衡方程：
ki(b)ui(b)＝fi(ti(b),b)
mi(b)ti(b)＝qi(b)
其中结构位移场ui包含温度场ti引起的位移。ki为层合板整体刚度区间矩阵，mi为整体热传导区间矩阵，fi为力载荷区间向量，包括温度场引起的外载荷以及机械载荷，qi为热载荷区间向量。使用区间参数顶点组合法，求出结构力热耦合响应在不确定性参数影响下的上下界，即：
其中，ui(b)为位移响应下界，
步骤三：基于结构非概率集合可靠性模型，将层合板位移约束处的实际位移区间(uj)i和许用位移区间
其中δuj和
φ(ρ,b)＝u*(b)-u(ρ,b)
来判断结构是否安全，当φ(ρ,b)<0时结构失效，当φ(ρ,b)≥0时结构安全。
步骤四：采用优化特征距离d作为衡量结构是否失效的非概率可靠性指标，其定义为：原失效平面到目标失效平面的移动距离，其中目标失效平面是与原失效平面平行的平面，并且其可靠度为一给定值rtarg，第j个约束的优化特征距离dj可以量化层合板满足该约束的非概率可靠度，表达式如下：
步骤五：采用基于simp方法的变密度拓扑优化模型，使用第四步得到的非概率可靠度指标即特征距离dj作为约束，建立复合材料层合结构非概率可靠性拓扑优化列式：
find:ρ＝(ρ1,ρ2,…,ρn)
s.t.k(ρ,b)u(ρ,b)＝f(t(ρ,b),ρ,b)
m(ρ,b)t(ρ,b)＝q(ρ,b)
dj(ρ,b)≤0,j＝1,2,…,m
0≤ρ≤ρi≤1
其中ρ＝(ρ1,ρ2,…,ρn)为单元伪密度，即拓扑优化问题的设计变量，n表示设计变量的个数，w表示结构重量，ρ为设计变量的下界，是为了防止刚度矩阵奇异而预设的一个小值。
步骤六：采用伴随向量法求解层合结构位移约束处响应上下界的灵敏度，再根据复合函数的链式求导法则得到非概率可靠性约束dj的灵敏度；
步骤七：采用移动渐近线优化算法(methodofmovingasymptotes)对设计变量进行迭代求解，在迭代过程中，如果当前设计不满足可靠度约束dj<0，或前后两个迭代步的设计变量变化绝对值之和大于容差ε时，则迭代步数加一，并返回第二步，否则，进行下一步；
步骤八：如果当前设计满足可靠性优化特征距离dj<0，而且前后两个迭代步的设计变量变化绝对值之和小于容差ε时，则迭代结束，得到层合结构非概率可靠性拓扑优化设计的最优构型。
本发明与现有技术相比的优点在于：本发明提供了一种针对复合材料层合结构的考虑力热耦合载荷影响下的非概率可靠性拓扑优化新方法，一方面，通过稳态热传导计算结构温度场，考虑温度场和位移场的耦合效应，弥补和完善了单一载荷形式下传统安全系数法的局限性。另一方面，本发明建立的非概率可靠性拓扑优化框架，为小样本试验数据条件下层合结构拓扑优化设计提供了有效的解决途径。
附图说明
图1是本发明针对力热耦合环境下复合材料层合结构的非概率可靠性拓扑优化方法流程图；
图2是本发明中所用到的承受力热耦合载荷的复合材料结构示意图
图3是本发明中所用到的单向板热导率计算示意图；
图4是本发明中所用到的经典层合板示意图；
图5是本发明中所用到的层合板位移约束处的实际位移响应区间和许用位移区间的一维干涉模型；
图6是本发明中所用到的参数标准化后的非概率可靠性模型的二维干涉示意图；
图7是本发明所用的非概率可靠性指标即优化特征距离dj的示意图。
具体实施方式
下面将结合本发明的附图，对本发明的技术方案进行清楚、完整的描述。
如图1所示，本发明提出了一种力热耦合环境下层合板非概率可靠性拓扑优化方法，对象为图2所示的承受力热耦合载荷的复合材料层合结构，其中，ω为设计域，q为热流率、热流密度等热载荷，fm为机械载荷，△t为热载荷引起的结构内部温度场与初始温度场的差。本发明包括以下步骤：
(1)根据给定的复合材料单向板材料参数和铺层方式，建立层合板力热耦合本构模型，利用一阶层合理论，计算层合板等效热导率，通过下式定义：
其中△q表示热量，△t表示时间，
mc,nn＝∑ωθ(mc,pcos2θ mc,tsin2θ)(2)
式中，mc,p表示复合材料单向板沿着纤维方向的热导率，mc,t表示复合材料单向板垂直于纤维方向的热导率。根据复合材料经典层合理论，温度沿着层合板厚度方向线性变化，即：
t(x,y,z)＝t0(x,y) tz'(x,y)·z(3)
其中t0(x,y)表示层合结构中面温度，t′z(x,y)表示沿层合板厚度z方向的温度梯度。
以对称铺层的层合板为例，全局坐标系下第k层的应力通过下式得到：
其中σx,σy和τxy表示全局坐标系下的应力分量，εx，εy和γxy表示全局坐标系下的应变分量，αx，αy和αxy表示全局坐标系下的热膨胀系数。
式中e1,，e2，g12，μ21和μ12分别是单向板主方向1的杨氏模量、主方向2的杨氏模量，主方向1和2的剪切模量以及主方向1和2的泊松比。tk表示层合板第k层的转换矩阵，即：
式中，θk表示第k层的纤维角度。
复合材料层合板的单元弹性矩阵
其中，tk为层合板第k层的厚度，h为层合板的总厚度。在结构有限元法中，单元刚度矩阵ki可通过下式得到：
其中，b为单元应变-位移矩阵。复合材料层合板的总体刚度矩阵可通过单元刚度矩阵求和得到：
(2)根据实际工程中复杂服役环境和层合板结构材料属性，选取描述层合板结构特性或外部力热载荷的不确定性参数b，根据有限元平衡方程
其中结构位移场ui包含温度场ti引起的位移。ki为层合板整体刚度区间矩阵，mi为整体热传导区间矩阵，fi为力载荷区间向量，包括温度场引起的外载荷以及机械载荷，qi为热载荷区间向量。计算结构力热耦合响应时，考虑层合结构温度场和位移场的单向耦合，即先通过稳态热传导计算结构的温度场分布，接着以此温度场计算温度引起的外载荷fth，最后将fth与温度无关的机械载荷fth叠加，根据层合板的热弹性理论计算结构的位移场。
使用区间参数顶点组合法，求出结构力热耦合响应在不确定性参数影响下的上下界，即：
其中，ui(b)为位移响应下界，
(3)基于结构非概率集合可靠性模型，将图5所示的层合板位移约束处的实际位移区间(uj)i和许用位移区间
其中δuj和
φ(ρ,b)＝u*(b)-u(ρ,b)(13)
来判断结构是否安全，如图6所示，当φ(ρ,b)<0时结构失效，当φ(ρ,b)≥0时结构安全。
(4)采用优化特征距离d作为衡量结构是否失效的非概率可靠性指标，其定义为：原失效平面到目标失效平面的移动距离，其中目标失效平面是与原失效平面平行的平面，并且其可靠度rtarg为一给定值，第j个约束的优化特征距离dj可以量化层合板满足该约束的非概率可靠度，表达式如下：
如图7所示，当dj>0时，失效平面在与目标可靠度rtarg对应的目标失效平面上方，此时由于安全区域的面积小于目标值，不满足要求。当dj≤0时，失效平面在与目标非概率可靠度rtarg对应的目标失效平面下方，此时由于安全区域的面积大于等于目标值，，满足设计要求；
(5)基于simp材料插值模型，将第i个单元的弹性模量、剪切模量、热导率表示为完全实心的层合板材料相对密度的函数：
其中p是惩罚因子，用于实现对中间密度单元的惩罚，一般取p＝3。
使用第(4)步得到的非概率可靠度指标即特征距离dj作为约束，建立基于simp方法的变密度拓扑优化模型复合材料层合结构非概率可靠性拓扑优化列式：
其中ρ＝(ρ1,ρ2,…,ρn)为单元伪密度，即拓扑优化问题的设计变量，n表示设计变量的个数，w表示结构重量，ρ为设计变量的下界，是为了防止刚度矩阵奇异而预设的一个小值。
(6)本发明采用移动渐进线方法(mma，methodofmovingasymptotes)这一梯度优化算法求解优化问题，因此需要获取目标函数和约束函数对设计变量的偏导数，即进行灵敏度分析。由于本发明的设计变量远多于约束函数个数，采用差分的方式求解灵敏度会带来巨额计算量。因此本发明先采用伴随向量法求解层合结构位移约束处响应上下界的灵敏度，再根据复合函数的链式求导法则得到非概率可靠性约束dj的灵敏度。
构建增广拉格朗日函数：
其中λj1为引入的拉格朗日乘子，(fm)ι和(fth)ι分别为机械载荷和热载荷引起的力区间向量，
上式对任意的λj1均成立，为了避免计算
即：
其中，
可以看到式在形式上和平衡方程kiui＝fi一致，因此为了求解伴随向量λj1，可以仅在第j个节点自由度处施加单位载荷通过一次有限元计算即可。
将λj1代入式得：
式中，
任一拓扑构型下区间总体刚度ki可写为：
因此有：
利用复合函数的链式求导法则，可靠性指标dj相对于设计变量ρi的灵敏度可以表示为：
结合式，有
和
其中，
此外，目标函数w对设计变量的偏导数为：
至此，目标函数和约束关于设计变量的灵敏度均已得到，可以调用梯度算法更新设计变量。
(7)采用移动渐近线优化算法(mma，methodofmovingasymptotes)对设计变量进行迭代求解，在迭代过程中，如果当前设计不满足可靠度约束dj<0，或前后两个迭代步的设计变量变化绝对值之和大于容差ε时，则迭代步数加一，并返回第(2)步，否则，执行第(8)步；
(8)同时考虑可靠度约束和相对变化量，如果当前设计满足可靠性优化特征距离dj<0，而且前后两个迭代步的设计变量变化绝对值之和小于容差ε时，则迭代结束，得到层合结构非概率可靠性拓扑优化设计的最优构型。
综上所述，本发明提出了一种力热耦合环境下层合板非概率可靠性拓扑优化方法。该方法首先根据层合结构一阶剪切变形理论，计算复合材料层合板等效弹性参数和热属性参数。接着考虑层合板材料参数和载荷等不确定性因素影响，利用simp材料插值模型，以结构重量最小为优化目标，以单元的相对密度作为设计变量，以力热耦合载荷作用下关注位置的许用位移响应可靠度为约束，构建层合结构可靠性拓扑优化框架，采用移动渐近线优化算法，通过迭代获得力热耦合环境下层合结构的最优构型。本发明提供了一种针对复合材料层合结构的考虑力热耦合载荷影响下的非概率可靠性拓扑优化新方法，一方面，通过稳态热传导计算结构温度场，考虑温度场和位移场的耦合效应，弥补和完善了单一载荷形式下传统安全系数法的局限性。另一方面，本发明建立的非概率可靠性拓扑优化框架，为小样本试验数据条件下层合结构拓扑优化设计提供了有效的解决途径。
以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制；其可扩展应用于力热耦合环境下的结构优化设计领域，凡采用等同变换或者等效替换而形成的技术方案，均落在本发明权利保护范围之内。
本发明未详细阐述部分属于本领域技术人员的公知技术。
1.一种力热耦合环境下层合板非概率可靠性拓扑优化方法，其特征在于，包括如下步骤：
第一步：根据给定的复合材料单向板材料参数和铺层方式，建立层合板力热耦合本构模型，利用一阶层合理论，计算层合板等效热导率、等效热传导矩阵和等效刚度矩阵；
第二步：根据实际工程中服役环境和层合板结构材料属性，选取描述层合板结构特性或外部力热载荷的不确定性参数b，根据有限元平衡方程：
ki(b)ui(b)＝fi(ti(b),b)
mi(b)ti(b)＝qi(b)
其中结构位移场ui包含温度场ti引起的位移；ki为层合板整体刚度区间矩阵，mi为整体热传导区间矩阵，fi为力载荷区间向量，qi为热载荷区间向量；使用区间参数顶点组合法，求出结构力热耦合响应在不确定性参数影响下的上、下界，即：
其中，ui(b)为位移响应下界，
第三步：基于结构非概率集合可靠性模型，将层合板位移约束处的实际位移区间(uj)i和许用位移区间
第四步：采用优化特征距离d作为衡量结构是否失效的非概率可靠性指标，其定义为：原失效平面到目标失效平面的移动距离，其中目标失效平面是与原失效平面平行的平面，并且其可靠度为一给定值，第j个约束的优化特征距离dj用于量化层合板满足该约束的非概率可靠度；
第五步：采用基于simp方法的变密度拓扑优化模型，使用第四步得到的非概率可靠度指标即特征距离dj作为约束，建立复合材料层合结构非概率可靠性拓扑优化列式；
第六步：采用伴随向量法求解层合结构位移约束处响应上下界的灵敏度，再根据复合函数的链式求导法则得到非概率可靠性约束dj的灵敏度；
第七步：采用移动渐近线优化算法对设计变量进行迭代求解，在迭代过程中，如果当前设计不满足可靠度约束dj<0，或前后两个迭代步的设计变量变化绝对值之和大于容差ε时，则迭代步数加一，并返回第二步，否则，进行下一步；
第八步：如果当前设计满足可靠性优化特征距离dj<0，而且前后两个迭代步的设计变量变化绝对值之和小于容差ε时，则迭代结束，得到层合结构非概率可靠性拓扑优化设计的最优构型。
2.根据权利要求1所述的一种力热耦合环境下层合板非概率可靠性拓扑优化方法，其特征在于：所述第一步中，热导率通过下式定义：
其中△q表示热量，△t表示时间，
mc,nn＝∑ωθ(mc,pcos2θ mc,tsin2θ)
式中，mc,p表示复合材料单向板沿着纤维方向的热导率，mc,t表示复合材料单向板垂直于纤维方向的热导率；根据复合材料经典层合理论，温度沿着层合板厚度方向线性变化，即：
t(x,y,z)＝t0(x,y) tz'(x,y)·z
其中t0(x,y)表示层合结构中面温度，tz'(x,y)表示沿层合板厚度z方向的温度梯度。
3.根据权利要求1所述的一种力热耦合环境下层合板非概率可靠性拓扑优化方法，其特征在于：所述第一步中，对于对称铺层的层合板，全局坐标系下第k层的应力通过下式得到：
其中σx,σy和τxy表示全局坐标系下的应力分量，εx，εy和γxy表示全局坐标系下的应变分量，αx，αy和αxy表示全局坐标系下的热膨胀系数；
式中e1，e2分别是单向板主方向1的杨氏模量、主方向2的杨氏模量，g12是主方向1和2的剪切模量，μ21和μ12分别是主方向1和2的泊松比；tk表示层合板第k层的转换矩阵，即：
式中，θk表示第k层的纤维角度；
复合材料层合板的单元弹性矩阵
其中，tk为层合板第k层的厚度，h为层合板的总厚度；在结构有限元法中，单元刚度矩阵ki通过下式得到：
其中，b为单元应变-位移矩阵，复合材料层合板的总体刚度矩阵通过单元刚度矩阵求和得到：
4.根据权利要求1所述的一种力热耦合环境下层合板非概率可靠性拓扑优化方法，其特征在于：所述第二步中考虑了层合结构温度场和位移场的单向耦合，即先通过稳态热传导计算结构的温度场分布，接着以此温度场计算温度引起的外载荷fth，最后将fth与温度无关的机械载荷fth叠加，根据层合板的热弹性理论计算结构的位移场。
5.根据权利要求1所述的一种力热耦合环境下层合板非概率可靠性拓扑优化方法，其特征在于：所述第三步中，标准化实际位移区间和许用位移区间的标准化方法如下：
其中δuj和
φ(ρ,b)＝u*(b)-u(ρ,b)
来判断结构是否安全，当φ(ρ,b)<0时结构失效，当φ(ρ,b)≥0时结构安全，ρ为伪密度。
6.根据权利要求1所述的一种力热耦合环境下层合板非概率可靠性拓扑优化方法，其特征在于：所述第四步中，采用了优化特征距离d作为衡量结构是否失效的非概率可靠性指标，对第j个约束对应的优化特征距离dj，其表达式为：
式中rtarg表示给定的可靠度指标。
7.根据权利要求1所述的一种力热耦合环境下层合板非概率可靠性拓扑优化方法，其特征在于：所述第五步中采用了simp(solidisotropicmaterialwithpenalization)模型来避免中间密度单元的产生，对于simp模型，第i个单元的弹性模量等参数是材料相对密度的函数：
其中p是惩罚因子，用于实现对中间密度单元的惩罚，取p＝3；
建立复合材料层合结构非概率可靠性拓扑优化列式：
find:ρ＝(ρ1,ρ2,…,ρn)
s.t.k(ρ,b)u(ρ,b)＝f(t(ρ,b),ρ,b)
m(ρ,b)t(ρ,b)＝q(ρ,b)
dj(ρ,b)≤0,j＝1,2,…,m
0≤ρ≤ρi≤1
其中ρ＝(ρ1,ρ2,…,ρn)为单元伪密度，即拓扑优化问题的设计变量，n表示设计变量的个数，w表示结构重量，ρ为设计变量的下界，是为了防止刚度矩阵奇异而预设的一个预定小值。
8.根据权利要求1所述的一种力热耦合环境下层合板非概率可靠性拓扑优化方法，其特征在于：所述第六步采用伴随向量法求解层合结构位移约束处响应上下界的灵敏度具体过程如下：构建增广拉格朗日函数：
其中
上式对任意的
即：
其中
将
式中，
任一拓扑构型下区间总体刚度ki写为：
因此有：
9.根据权利要求1所述的一种力热耦合环境下层合板非概率可靠性拓扑优化方法，其特征在于：所述第六步中利用复合函数的链式求导法则得到非概率可靠性约束dj的灵敏度具体过程如下：可靠性指标dj相对于设计变量ρi的灵敏度表示为：
结合可靠性指标dj的表达式，有：
和
其中，
