一种基于渐近变分法的复合材料层合板仿真及优化方法

摘要

本发明涉及一种基于渐近变分法的复合材料层合板仿真及优化方法，属于材料力学分析领域。具体包括以下步骤：基于旋转张量分解概念构建由一维广义应变和翘曲函数表示的三维板能量方程；基于渐近变分法将原三维问题分析严格拆分为非线性二维板分析和沿厚度方向的截面分析；利用板固有的小参数将降维模型的近似能量泛函渐近修正到第二阶，并通过平衡方程将近似能量泛函转换为Reissner模型形式以便于实际应用；由得到的全局响应的渐近修正翘曲函数精确重构三维应力/应变/位移场；采用极大化截面分析得到的弯曲和扭转刚度系数的优化策略对复合层合板进行优化。本发明实用性强，通用性高，可显著提高此类问题的解算速度和效率。

说明书

技术领域

本发明属于材料力学性能分析领域，尤其是一种能够有效模拟和准确重构复合材料层合板三维应力/应变/变形场和极大化刚度矩阵的弯曲和扭转项为优化策略、层合复合材料铺设倾角为设计变量的优化设计方法。

背景技术

近20年来，先进复合材料结构因其高强度、高模量、可设计性等优点已广泛应用于航天航空、机械、土木等领域。复合材料结构很多是厚度尺寸比其它两个方向尺寸小得多的平板铺层结构，复合材料铺层设计的灵活性为结构的应力计算、变形分析以及强度预测带来额外的复杂性。基于Kirchhoff假设(薄板)、由三维弹性理论推导出的古典层合理论是最简单的复合板分析理论，但若计入板厚，其精度较三维有限元分析下降很多。各国学者利用板厚相对基准面上的变形很小的特点，将板的厚度坐标从偏微分方程的独立变量中去掉，用精确的二维模型代替三维模型以对古典模型不足之处进行改进，但大部分研究都是基于特定的动力学假设(如沿厚度方向的位移分布假设等)，无法较好地反映层合板三维效应和铺层之间的相互作用。

随着近年来计算硬件和软件的发展，可使用三维有限元软件ANSYS或NASTRAN对这类复合层合板进行三维计算和分析，如王跃进等以ABAQUS有限元软件为平台，对三维复合层合板渐近损伤进行非线性分析；Apalak的分析模型中对复合材料层合板采用ANSYS软件中的八节点三维层状体单元进行应力求解。但仍有必要对这种结构进行简化分析，原因在于：首先，许多工程问题目前无法用三维有限元分析处理，如实际的转子叶片由200多层非常薄的复合层构成，若使用三维有限元模型，每一铺层至少需要一个单元进行模拟，该叶片模型很容易超过109个自由度。因此，由4个转子叶片构成的旋翼气弹性分析目前无法在任一台计算机上通过三维有限元分析解决；二是在对整个结构进行有限元分析的起始准备工作，包括材料种类的选择、层合和夹层结构的设计以及层合板铺设方式设计，有限元软件的实用性不大，当在层合板结构的不同层的级别上来对复合材料行为(应力、应变分布)进行细节研究时，有限元软件提供的后处理能力尤其有限。另外，尽管三维有限元在连续介质力学框架内提供很好的准确性，但计算太费时，且耗费大量的计算机资源，难以在规定的设计和分析时间内完成预定目标，因此迫切需要一种高效、快速的专业复合材料层合板分析方法和程序，以缩短设计时间，降低计算成本。本专利技术就是在这样的背景下展开的。

发明内容

针对现有技术存在分析效率低、精度差不足，尤其是沿厚度方向的应力分布无法精确预测之不足，本发明的目的是提供一种计算量小，占用计算机资源少，且效率高的基于渐近变分方法的复合材料层合板三维场仿真及优化方法。

本发明采用的技术方案如下：一种基于渐近变分法的复合材料层合板仿真及优化方法，包括如下步骤：

1)基于旋转张量分解方法将复合材料层合板三维连续方程用一维广义应变和翘曲函数表示；

2)将非线性三维复合层合板模型分解为二维截面分析和沿法线方向的一维非线性分析；采用渐近变分算法实现简化模型的求解过程：通过二维截面分析得到本构关系，即刚度矩阵和截面特性和渐近修正翘曲函数；通过沿结线方向的一维非线性分析得到二维板的全局响应性能；

3)基于推导的重构关系，由计算得到的二维板变形、全局旋转张量分量和渐近修正翘曲函数重构三维变形；由一维广义应变和渐近修正翘曲函数重构三维应变场；使用材料本构关系得到三维应力；

4)以弯扭屈曲荷载最大为目标，以不发生屈曲为约束条件，层合复合材料铺设倾角为设计变量的复合梁优化设计。

相比现有技术，本发明具有如下有益效果：

本发明是一种有效模拟和准确重构复合材料层合板三维应力/应变/变形场和极大化刚度矩阵的弯曲和扭转项为优化策略、层合复合材料铺设倾角为设计变量的优化设计方法，借助渐近变分技术提供一种高效、快速、实用性强、通用性强，能够通过渐近修正翘曲函数建立复合材料层合板简化数值模型和三维重构关系的方法。

本发明首先建立用一维广义应变和翘曲函数表示的三维板能量方程数值模式；然后将渐近变分引入到复合材料层合板简化数值建模问题研究中，将原三维问题分析严格拆分为二维截面分析(等效单层板模型)和沿法线方向的一维非线性分析。本发明通过二维截面分析得到沿法线一维非线性分析所需的本构关系(刚度矩阵和截面特性)和渐近修正翘曲函数；利用一维非线性分析得到的板全局响应性能和翘曲函数精确重构三维场；采用极大化二维截面分析得到的弯曲和扭转刚度系数的优化策略对复合层合板进行优化，以达到通过改变复合材料铺设倾角的方法提高复合板所承受的最大弯扭屈曲荷载的目的。极大化刚度系数的优化策略能有效提高复合梁的屈曲临界荷载，并可节约优化搜索的时间和步骤，为工程技术人员提供了有效的设计手段。本发明实用性强，通用性高，可显著提高此类问题的解算速度和效率。

附图说明

图1为本发明三维场仿真及优化方法流程图；

图2为板变形前后坐标系图；

图3为复合材料层合板柱形弯曲问题模型图；

图4为σ₁₁随厚度坐标的变化关系图；

图5为σ₁₂随厚度坐标的变化关系图；

图6为σ₂₂随厚度坐标的变化关系图；

图7为σ₁₃随厚度坐标的变化关系图；

图8为σ₂₃随厚度坐标的变化关系图；

图9为σ₃₃随厚度坐标的变化关系图；

图10为ID3在极大化刚度系数优化策略下设计变量和目标函数的评估。

具体实施方式

下面结合具体实例和附图对本发明作进一步说明。

1、三维表达式

参见图2，是板变形前后坐标系图。复合层合板基准面上任意一点位置可由其笛卡尔坐标系x_i表示，其中，x_α是基准面上相互正交的坐标，x₃是法向坐标(字母下标i、j、k代表1、2、3；字母下标α、β代表1、2，下同)。引入一组沿x_i方向的正交参考坐标向量b_i，变形前板上任一点的位置可由固定点O到x_α确定点的位置向量描述：

$\hat{r}(x_1,x_2,x_3)=r(x_1,x_2)+x_3{b}_3---\left(1\right)$

当定义板的中间层为其基准面时

$<\hat{r}(x_1,x_2,x_3)>=r(x_1,x_2)---\left(2\right)$

其中尖括号表示沿板厚度的积分，下同。

板变形后，位置向量转换为变形后的位置向量，后者可由三维体的变形唯一确定。为推导方便，引入与变形板有关的另一组坐标系B_i。B_i和b_i之间的关系由方向余弦函数矩阵C(x₁，x₂)确定：

B_i＝C_ijb_j    C_ij＝B_i·b_j    (3)

变形板上任一点的位置向量可表示为

$\hat{R}(x_1,x_2,x_3)=R(x_1,x_2)+x_3{B}_3(x_1,x_2)+w_i(x_1,x_2,x_3)B_i(x_1,x_2)---\left(4\right)$

其中w_i是翘曲分量，在这里视为未知的三维函数求解，以考虑包括古典层合理论未考虑的局部翘曲变形在内的所有变形。翘曲的引入使式(4)有6次冗余，需6个约束来求解方程，可通过适当定义和B_i消除冗余，与式(2)类似，可定义在板的中间层，这样，翘曲函数必须满足3个约束：

<w_i(x₁，x₂，x₃)>＝0                         (5)

将B₃设为与板变形面垂直可确定其它两个约束。

基于旋转张量分解理论，由局部小旋转的条件，Jauman-Biot-Cauchy应变分量可表示为

Γ_ij＝1/2(F_ij+F_ji)-δ_ij                     (6)

其中，F_ij是变形梯度张量的混合基分量

假设应变很小，忽略一维广义应变和翘曲的乘积项，可得到一维广义应变和翘曲函数表示的三维应变场：

$\Gamma ={\Gamma }_{h}w+{\Gamma }_{ϵ}ϵ+{\Gamma }_{l_1}{w}_{,1}+{\Gamma }_{l_2}{w}_{,2}---\left(7\right)$

其中：w_，1表示w对x_i的偏导，同理，下同；

Γ＝[Γ₁₁ 2Γ₁₂ 2Γ₂₂ 2Γ₁₃ 2Γ₂₃ Γ₃₃]^T    (8)

w＝[w₁ w₂ w₃]^T                              (9)

ε＝[ε₁₁ 2ε₁₂ ε₂₂ K₁₁ K₁₂+K₂₁ K₂₂]^T      (10)

式中Γ_h，Γ_ε，为积分算子。

板每单位面积上的应变能可表示为

$U=\frac{1}{2}<{\Gamma }^T\mathrm{D\Gamma }>---\left(11\right)$

其中D是6×6阶对称材料参数矩阵。

由于应变很小，可安全地忽略虚拟旋转中翘曲与荷载的乘积项，得到式(5)约束下荷载τ_iB_i，β_iB_i，φ_iB_i(分别作用在板的顶面、底面和厚度方向)虚拟位移产生的虚功：

W＝-τ^Tw|_x3＝h/2-β^Tw|_x3＝-h/2-<φ^Tw>       (12)

根据虚功原理，层合板的三维能量问题的完整表述为

δU-δW＝0                                  (13)

总能量的泛函为

∏＝U+W                                     (14)

δ∏＝0                                     (15)

在(5)式约束下的总能量中，仅翘曲是变化的。未知的翘曲函数可由对总能量最小化求得。

(14)式仅是原三维板弹性问题的另一表达形式。若直接进行求解，将会遇到与原三维板弹性问题相同的困难。若板由多层复合层构成，该计算过程将变得十分繁琐。本文借助渐近变分法渐近计算三维翘曲函数，可使计算得到简化。为考虑任何横截面形状和各向异性材料，首先应用有限元方法将三维翘曲场离散为一维有限元单元形式

w(x₁，x₂，x₃)＝S(x₃)V(x₁，x₂)              (16)

其中：S(x₃)表示单元形函数，V为横向法线方向翘曲场的节点位移。

将(16)代入(13)中，可得到离散形式的总能量为

$2\Pi =V^T\mathrm{EV}+2V^T(D_{\mathrm{hϵ}}ϵ+D_{h{l}_1}{V}_{,1}+D_{{\mathrm{hl}}_2}{V}_{,2})$

$+{ϵ}^T{D}_{\mathrm{ϵϵ}}ϵ+V_{,1}^T{D}_{l_1{l}_1}{V}_{,1}+V_{,2}^T{D}_{l_2{l}_2}{V}_{,2}---\left(17\right)$$+2(V_{,1}^T{D}_{l_1ϵ}ϵ+V_{,2}^T{D}_{l_2ϵ}ϵ+V_{,1}^T{D}_{l_1{l}_2}{V}_2)+2V^{T}L$

其中L＝-S^+Tτ-S^-Tβ-<S^Tφ>为荷载相关项。

新引入的几何形状和材料属性有关的变量包括：

E＝<[Γ_hS]^T D[Γ_hS]>，D_hε＝<[Γ_hS]^T DΓ_ε>

$D_{h{l}_1}=<{\left[{\Gamma }_{h}S\right]}^{T}D\left[{\Gamma }_{l_1}S\right]>,$$D_{{\mathrm{hl}}_2}=<{\left[{\Gamma }_{h}S\right]}^{T}D\left[{\Gamma }_{l_2}S\right]>$

D_εε＝<Γ_ε^T DΓ_ε>，$D_{l_1{l}_1}=<{\left[{\Gamma }_{l_1}S\right]}^{T}D\left[{\Gamma }_{l_1}S\right]>---\left(18\right)$

$D_{{l_{1}l}_2}=<{\left[{\Gamma }_{l_1}S\right]}^{T}D\left[{\Gamma }_{l_2}S\right]>,$$D_{{l_{2}l}_2}=<{\left[{\Gamma }_{l_2}S\right]}^{T}D\left[{\Gamma }_{l_2}S\right]>$

$D_{l_1ϵ}=<{\left[{\Gamma }_{l_1}S\right]}^{T}D{\Gamma }_{ϵ}>,$$D_{l_2ϵ}=<{\left[{\Gamma }_{l_2}S\right]}^{T}D{\Gamma }_{ϵ}>$

式(5)的离散形式为

V^THψ＝0              (19)

其中H＝<S^TS>和ψ是E的正交化核心矩阵，ψ^THψ＝1。三维翘曲函数的求解问题转化为在式(19)约束下(17)式最小化问题。

2、降维后的近似能量公式推导

从所周知，弹性体的形态完全由其能量所确定。渐近变分法是一种有效的数学工具，可用二维能量尽可能准确地再现三维能量。对复合层合板结构，可选择h/l(h为板厚，l为板翘曲变形值)和广义二维应变的阶数ε作为渐近变分计算所需的小参数。

应用渐近变分法，首先需根据泛函的不同阶数找到主导项。由于总能量中只有翘曲是变化的，只需找到仅含有翘曲的主导项或翘曲与其它量乘积项(如广义应变和荷载)的主导项。

对式(17)零阶近似后泛函的主导项为

$2{\Pi }_0^{*}=V^T\mathrm{EV}+2V^T{D}_{\mathrm{hϵ}}ϵ---\left(20\right)$

在式(19)约束下(20)式的欧拉-拉格朗日方程可借助拉格朗日乘子Λ通过对变量的常规计算得到

EV+D_hεε＝HψΛ      (21)

考虑内核矩阵ψ的属性，拉格朗日乘子Λ为

A＝ψ^TD_hεε          (22)

将式(22)代回式(21)，可得

EV＝(Hψψ^T-I)D_hεε  (23)

由于式(23)等式右边与V的零空间E正交，存在唯一解与零空间E线性独立，可选择任何约束以方便求解：

V＝V^*+ψλ            (24)

其中：V^*为线性体系解，λ可由式(19)确定

λ＝-ψ^THV^*           (25)

式(25)代入式(24)，则式(19)约束下式(17)的最小化解为

$V=(I-\psi {\psi }^{T}H)V^{*}={\hat{V}}_0ϵ=V_0---\left(26\right)$

将式(26)代入式(20)得到渐近修正到零阶近似的总能量泛函为

$2{\Pi }_0={ϵ}^T({\hat{V}}_0^T{D}_{\mathrm{hϵ}}+D_{\mathrm{ϵϵ}})ϵ---\left(27\right)$

这一近似能量与古典分层板理论得出的结果相吻合，但未建立在Kirchhoff动力学假设基础上，且尽管能量相同，但推导出的横向法线方向应变并不为零。注意到零阶近似得到的翘曲是ε阶，根据渐近变分法，为接受零阶近似结果，必须检查下一阶近似是否比这还要高。为得到一阶近似翘曲，对零阶近似翘曲进行如下摄动

V＝V₀+V₁              (28)

将式(28)代入式(7)及(17)，可得到总能量泛函一阶近似的主导项

$2{\Pi }_1^{*}=V_1^{T}E{V}_1+2V_{\alpha }^T{D}_{\alpha }{ϵ}_{,\alpha }+2V_1^{T}L---\left(29\right)$

其中L为荷载相关项。

同零阶近似一样，可求解一阶翘曲场

V₁＝V_1αε_，α+V_1L    (30)

代入(29)式，可得到渐近修正到二阶的总能量泛函：

$2{\Pi }_1^{*}={ϵ}^T\mathrm{Aϵ}+{ϵ}_{,1}^{T}B{ϵ}_{,1}+2{ϵ}_{,1}^T{\mathrm{Cϵ}}_{,2}+{ϵ}_{,2}^T{\mathrm{Dϵ}}_{,2}+2{ϵ}^{T}F+P---\left(31\right)$

其中：

$A={\hat{V}}_0^T{D}_{\mathrm{hϵ}}+D_{\mathrm{ϵϵ}},$$B={\hat{V}}_0^T{D}_{l_1{l}_1}{\hat{V}}_0+V_{11}^T{D}_1,$$C={\hat{V}}_0^T{D}_{l_1{l}_2}{\hat{V}}_0+1/2(V_{11}^T{D}_2+D_1^T{V}_{12}),$$C={\hat{V}}_0^T{D}_{l_1{l}_2}{\hat{V}}_0+V_{12}^T{D}_2$

                                                   (32)

$F={\hat{V}}_0^{T}L-1/2(D_1^T{V}_{1L,1}+V_{11}^T{L}_{,1}+D_2^T{V}_{1L,2}+V_{12}^T{L}_{,2})$$P=V_{1L}^{T}L$

3、将近似能量转换为Reissner模型形式

尽管总能量经二次渐近修正，但因包含广义应变量的导数，涉及到更加复杂的边界条件。为得到实用的能量方程，可将这种近似能量转换为Reissner模型形式。

在Reissner模型中有两个增加的横向剪切应变自由度γ＝[2γ₁₃ 2γ₂₃]^T，将其纳入横向法线方向的旋转变量中。古典应变量的Reissner形式可表示为

ε＝R-D_αγ_，α       (33)

其中：$D_1={\left(\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\end{array}\right)}^T,$$D_2={\left(\begin{array}{cccccc}0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)}^T$

$R={\left(\begin{array}{cccccc}{ϵ}_{11}^{*}& 2{ϵ}_{12}^{*}& {ϵ}_{22}^{*}& K_{11}^{*}& K_{12}^{*}+K_{21}^{*}& K_{22}^{*}\end{array}\right)}^T---\left(34\right)$

将(33)式代入(31)式，可得到Reissner应变形式表示的修正到第二阶的总能量泛函

$2{\Pi }_R=R^T\mathrm{AR}-2R^{T}A{D}_{\alpha }{\gamma }_{,\alpha }+R_{,1}^{T}B{R}_{,1}+2R_{,1}^T{\mathrm{CR}}_{,2}+R_{,2}^T{\mathrm{DR}}_{,2}+2R^{T}F+P---\left(35\right)$

4、三维场重构关系推导

经渐近修正的Reissner模型形式的近似能量尽可能接近总能。但这还远不足以准确重现复合层合板的三维位移场，还需提供重构关系以完善降维模型。本文通过得到二维变形和翘曲函数准确重构复合层合板三维场。

对经二阶渐近修正的应变能，可基于渐近修正的严格定义，由式(1)，(3)，(4)得到一阶渐近修正的三维变形

$U_{3d}=u_{2d}+x_3\left(\begin{array}{c}{C}_{31}\\ C_{32}\\ C_{33}-1\end{array}\right)+{\mathrm{SV}}_0+S{\bar{V}}_1---\left(36\right)$

其中：U_3d是三维变形的列阵；u_2d是板的二维变形；C_ij由式(3)得到的全局旋转张量分量。

由式(7)可得到一阶渐近修正的三维应变场

$\Gamma ={\Gamma }_{h}S(V_0+{\bar{V}}_1)+{\Gamma }_{ϵ}ϵ+{\Gamma }_{l_1}S{V}_{\mathrm{0,1}}+{\Gamma }_{l_2}{\mathrm{SV}}_{\mathrm{0,2}}---\left(37\right)$

然后，可使用材料本构关系得到三维应力σ_ij。

5、优化问题和算法

优化的主要目的是寻找复合层合梁最优铺层倾角，以提高弯扭屈曲载荷。在这里考虑在屈曲临界荷载约束下极大化刚度矩阵的弯曲和扭转项。由刚度矩阵的扭转项s₄₄和弯曲项s₅₅作为目标函数

minmize    ψ₀(θ)＝1/(D₄₄+D₅₅)    (38)

优化算法选用以梯度为基础的序列二次规划算法，在当前设计点将目标函数做二阶近似，求解后得到新的设计变量，重复这一过程直至迭代收敛。该方法在连续的、凸设计空间被证明是有效的，在非凸空间可能会收敛到局部最优解，可从不同出发点进行优化设计以得到全局收敛解。

为了验证上述方法的可靠性，对一具有20层对称复合层合板([30°/-30°/30°/-30°]₅，从底部到顶部编号)柱形弯曲问题进行分析。板长L＝10cm(沿x₁轴方向)，板厚h＝2.5cm，长厚比L/h＝4。板两端简支，承受正弦变化荷载：

${\tau }_3={\beta }_3=\frac{p_0}{2}\sin \left(\frac{\pi x_1}{L}\right);$τ_α＝β_α＝0                          (39)

参见图4～9，为重构的沿厚度方向6个应力分量与精确解的对比结果。可看出，所有模型的平面应力分量σ₁₁，σ₁₂，σ₂₂均与精确解吻合较好；对于面外横向应力分量σ₁₃，σ₂₃，本方法解可获得比古典层合理论和一阶剪切变形理论解更精确的结果；对于面外横向正应力σ₃₃，古典层合理论和一阶剪切变形理论解皆无法获得满意的数值计算结果，而本方法解与精确解一致，且可在一般台式机上1s完成所有计算，而三维有限元计算需30s(台式机配置为P41.7G，512M内存，320G硬盘)。

优化过程所需要的初始设计变量和屈曲荷载列于表2，共考虑三种不同的初始铺层倾角：ID1、ID2、ID3。设计变量的优化范围定义为[-90°，90°]。采用极大化刚度系数的优化策略、铺层倾角为设计变量、屈曲荷载最大为优化目标进行优化，优化结果见表3。对比分析表明：ID1、D2、ID3的屈曲荷载经优化后都有明显增加，ID3的屈曲荷载由初始的400N增加到586N，增幅达46.5％；优化过程中使用不同的初始设计会产生不同的优化设计结果，这是因为在非凸问题中使用确定性算法的结果。在这种情况下，优化往往只能找到局部最优解，而不是所希望的整体优化设计。可通过从几个不同的初始点进行优化设计以获取最佳解。

参见图10，为极大化刚度系数的优化策略下ID3的目标函数和设计变量随迭代次数的收敛情况。从图10(a)可看出，经第一次函数调用后目标函数就基本收敛到最优解收敛速度很快。由图10(b)可看到，设计变量θ₂由60°变化到20°，θ₄由-45°变化到-5°，且基本上都是在第一次函数调用后就收敛到最佳倾角值，大大节约了优化搜索的时间和步骤，同时也表明该优化策略的高效性。

表1为复合材料工程弹性常数表；

表2为复合层合板的初始设计变量和屈曲荷载；

表3为使用极大化刚度系数优化策略进行优化设计结果。

表1

表2

表3