CN110532688A - 一种曲纤维铺层结构的力学性能分析方法 - Google Patents

Abstract

本发明提供一种曲纤维铺层结构的力学性能分析方法，包括以下步骤：定义曲纤维铺层平板结构中连续变角度纤维的铺设路径；构造曲纤维铺层平板结构的有限元模拟的四边形板壳单元；构造曲纤维铺层平板结构的非线性有限元全阶模型；构造曲纤维铺层平板结构的非线性有限元降阶模型；求解的曲纤维铺层平板结构的非线性有限元降阶模型，得到曲纤维铺层平板结构的非线性力学承载响应曲线。本发明提出的一种曲纤维铺层结构的力学性能分析方法，解决现有技术中曲纤维铺层结构离散后形成的有限元模型规模过于庞大，进而显著影响其非线性力学分析效率的问题；本发明能够实现高效准确的对曲纤维铺层结构的力学性能的非线性分析。

Description

一种曲纤维铺层结构的力学性能分析方法

技术领域

本发明属于结构力学建模与分析技术领域，具体涉及一种曲纤维铺层结构的力学性能分析方法。

背景技术

近几年逐渐兴起的曲纤维复合材料是带有连续变角度纤维铺层的复合材料。不断变化的纤维取向导致了铺层刚度在不同的空间位置各不相同，因此设计人员可实现对曲纤维铺层结构的三维刚度剪裁，据此调整其内部载荷分布，实现提高结构承载性能。

为充分摸清曲纤维铺层结构的承载机理，通过对纤维曲线角的科学设计来达成改进其承载性能的目标，设计人员需要基于大量精细的结构数值计算来细致研究其在载荷作用下的整个非线性承载过程。然而，由于在有限元数值建模过程中需要采用足够小的网格尺寸来保证变角度纤维铺放的连续性，当前普遍采用的非线性数值建模计算技术由于计算效率问题已难以支持精细的结构非线性有限元力学分析更多地应用于曲纤维铺层结构的工程设计领域。

发明内容

针对现有技术存在的缺陷，本发明提供一种曲纤维铺层结构的力学性能分析方法，可有效解决上述问题。

本发明采用的技术方案如下：

本发明提供一种

本发明提供的一种曲纤维铺层结构的力学性能分析方法具有以下优点：

本发明提出的一种曲纤维铺层结构的力学性能分析方法，解决现有技术中曲纤维铺层结构离散后形成的有限元模型规模过于庞大，进而显著影响其非线性力学分析效率的问题；本发明能够实现高效准确的对曲纤维铺层结构的力学性能的非线性分析。

附图说明

图1为本发明提供的一种曲纤维铺层结构的力学性能分析方法的流程示意图。

图2为连续变角度曲纤维铺设路径的示意图；

图3为曲纤维铺层平板结构示意图；

图4为结构承载响应曲线示意图。

具体实施方式

为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

当前结构非线性力学求解技术的庞大计算量归结于需要对曲纤维铺层结构离散后的大规模有限元模型进行迭代求解，因此，为了能够有效缩减非线性计算中有限元模型的自由度数目，本发明提出一种针对曲纤维铺层结构的非线性有限元降阶模型建立方法，属于结构力学建模与分析技术领域，从而实现高效准确的对曲纤维铺层结构的力学性能的非线性分析。

参考附图1-4，曲纤维铺层结构的力学性能分析方法，可以克服现有技术中针对曲纤维铺层结构离散后形成的有限元模型规模过于庞大，进而显著影响其非线性力学分析效率的问题，本发明包括以下步骤：

步骤1，定义曲纤维铺层平板结构中连续变角度纤维的铺设路径，包括：

曲纤维铺层平板结构由多个曲纤维铺层叠加粘合而成；对于曲纤维铺层平板结构，以曲纤维铺层平板结构的中心为坐标原点θ₀，建立全局坐标系(x,y)；另外，定义参考坐标系(g,f)，参考坐标系(g,f)和全局坐标系(x,y)共相同的坐标原点θ₀，参考坐标系(g,f)为相对于全局坐标系逆时针旋转角度η得到的坐标系；曲纤维铺层平板结构的长度为d，则：在参考坐标系(g,f)下，对于通过坐标原点θ₀的曲纤维铺设路径，在距坐标原点θ₀的d/2处的曲纤维角度为θ₁，曲纤维铺设路径上任意一点的曲纤维角度表示为θ(g)，则得到参考坐标系(g,f)下的曲纤维铺设路径的描述方式为：

因此，曲纤维铺层平板结构的任意一个曲纤维铺层，由与通过坐标原点θ₀的曲纤维铺设路径平行的多个曲纤维铺设路径组成，因此，曲纤维铺层平板结构的任意一个曲纤维铺层，采用符号η<θ₀|θ₁>来描述；

步骤2，构造曲纤维铺层平板结构的有限元模拟的四边形板壳单元，并将所述曲纤维铺层平板结构离散为多个单元网格，计算每个单元网格的单元应变向量ε以及单元本构矩阵C_m，具体包括以下步骤：

步骤2.1，单元应变向量ε为单元线性应变向量ε_l与单元非线性应变向量ε_nl之和，即：

ε＝ε_l+ε_nl (2)

其中：单元线性应变向量ε_l为：

单元非线性应变向量ε_nl为：

其中：

为单元面内线性应变向量；

ε^b为单元面外线性应变向量；

为单元面内非线性应变向量；

B_l为单元线性几何插值矩阵，是常数矩阵；

B_nl(q^e)为关于单元节点位移向量q^e的单元非线性几何插值矩阵；

q^e为单元节点位移向量；

T代表矩阵的转置；

Z为单元非线性几何插值矩阵中的单元节点位移无关项，进一步表示为：

其中：

K_xx为单元非线性几何插值矩阵中节点位移无关项的xx向分量，是常数矩阵；

K_yy为单元非线性几何插值矩阵中节点位移无关项的yy向分量，是常数矩阵；

K_xy为单元非线性几何插值矩阵中节点位移无关项的xy向分量，是常数矩阵；

步骤2.2，根据参考坐标系(g,f)下的曲纤维铺设路径，构造曲纤维铺层平板结构的有限元模拟的四边形板壳单元，其单元本构关系为：

其中：

N为单元面内力向量；

M为单元弯矩向量；

A(θ)是单元面内刚度；

G(θ)是单元面外刚度；

D(θ)是单元面内/面外耦合刚度；

A(θ)、G(θ)和D(θ)的值受曲纤维铺设方式影响而随铺层面内曲纤维角度θ(g)的变化而变化；

ε^m是单元面内应变向量；

ε^b是单元面外应变向量；

根据单元本构关系，求得单元面内力向量N的表示式；

步骤2.3，获得单元本构矩阵C_m关于θ的表示式为：

步骤3，构造曲纤维铺层平板结构的非线性有限元全阶模型，包括：

步骤3.1，每个单元网格的单元应变能U关于q^e的表达式为：

其中：

i＝1,2,...,6，代表6个方向，分别为：x方向，y方向，z方向，xy方向，xz方向，yz方向；

j＝1,2,...,6；代表6个方向，分别为：x方向，y方向，z方向，xy方向，xz方向，yz方向；

A_s是单元网格的面积；

C_mij是单元本构矩阵C_m曲纤维角度中的第i行第j列的元素，单元本构矩阵C_m为6行6列矩阵；

ε_li是单元线性应变向量ε_l中的i向分量，即：根据公式3计算得到的单元线性应变向量ε_l中的i向分量，ε_li是关于q^e的表达式；

ε_lj是单元线性应变向量ε_l中的j向分量，即：根据公式3计算得到的单元线性应变向量ε_l中的j向分量，ε_lj是关于q^e的表达式；

ε_nli是单元非线性应变向量ε_nl中的i向分量，即：根据公式4计算得到的单元非线性应变向量ε_nl中的i向分量，ε_nli是关于q^e的表达式；

ε_nlj是单元非线性应变向量ε_nl中的j向分量，即：根据公式4计算得到的单元非线性应变向量ε_nl中的j向分量，ε_nlj是关于q^e的表达式；

步骤3.2，计算单元应变能U关于单元节点位移向量q^e的一阶导数，得到单元节点内力向量f^e；其中，单元节点是指每个单元网格的顶点；

其中：

为常规的单元线性刚度矩阵；

N_nl是单元面内力向量N的非线性部分；

T代表矩阵的转置；

步骤3.3，对曲纤维铺层平板结构的每个单元节点的单元节点内力向量f^e进行组装，获得曲纤维铺层平板结构的节点内力向量f(q)；

由此得到曲纤维铺层平板结构的非线性有限元全阶模型：

f(q)＝f_ext (10)

其中：f_ext为曲纤维铺层平板结构的节点外载荷向量；

曲纤维铺层平板结构的非线性有限元全阶模型，也被称作曲纤维铺层平板结构的非线性平衡方程，f_ext为曲纤维铺层平板结构的节点外载荷向量。上述结构非线性有限元全阶模型的自由度数为N_s。对于曲纤维铺层结构，其非线性有限元全阶模型的规模通常较大，即N_s>2000。采用常规的路径跟踪求解技术，如Newton-Raphson法或弧长法，求解结构非线性有限元全阶模型所构成的非线性平衡方程(10)即可获得结构的非线性承载响应曲线，进而得到结构的承载性能。

步骤4，构造曲纤维铺层平板结构的非线性有限元降阶模型，包括：

步骤4.1，计算单元应变能U关于单元节点位移向量q^e的二阶导数，得到第一种单元高阶物理量S^e：

其中：

N_x为单元面内力向量N的x向分量；

N_y为单元面内力向量N的y向分量；

N_xy为单元面内力向量N的xy向分量；

步骤4.2，计算单元应变能U关于单元节点位移向量q^e的三阶导数，将三阶层数与任意的以及相乘，得到第二种单元高阶物理量

其中：

下标α＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

下标β＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

代表与第α阶密集屈曲模态对应的单元节点位移向量；

代表与第β阶密集屈曲模态对应的单元节点位移向量；其中，β≠α；

为关于的单元非线性几何插值矩阵；

为关于的单元非线性几何插值矩阵；

为关于的单元面内力向量；

为关于的单元面内力向量；

为关于和的单元面内力向量；

步骤4.3，计算单元应变能U关于单元节点位移向量q^e的四阶导数，将四阶导数与任意的以及相乘，得到第三种单元高阶物理量

其中：

下标α＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

下标β＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

下标γ＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

下标δ＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

其中：α≠β≠γ≠δ；

代表与第γ阶密集屈曲模态对应的单元节点位移向量；

代表与第δ阶密集屈曲模态对应的单元节点位移向量；

同时定义符号：

B_nl( )为关于括号内参数的单元非线性几何插值矩阵；

N( )为关于括号内参数的单元面内力向量；

步骤4.4，在计算得到曲纤维铺层平板结构的每个单元对应的第一种单元高阶物理量S^e、第二种单元高阶物理量和第三种单元高阶物理量后，对曲纤维铺层平板结构的各个单元进行组装，获得曲纤维铺层平板结构的整体的第一种结构高阶物理量S、第二种结构高阶物理量Q(q_α,q_β)和第三种结构高阶物理量P(q_α,q_β,q_γ,q_δ)；

步骤4.5，将第一种结构高阶物理量S、第二种结构高阶物理量Q(q_α,q_β)和第三种结构高阶物理量P(q_α,q_β,q_γ,q_δ)带入到以下三组线性方程组，

其中：

矩阵Y为摄动矩阵，其各列向量分别由曲纤维铺层平板结构的前u阶密集屈曲模态与结构几何刚度矩阵的乘积计算获得；

单位基向量E_α，为与第α阶密集屈曲模态对应的单位基向量，α＝1,2,...,u，其第α个分量为1，而其余各分量均为零；

u_α,u_β,u_γ,u_δ分别为：与第α阶密集屈曲模态对应的结构一阶位移场；与第β阶密集屈曲模态对应的结构一阶位移场；与第γ阶密集屈曲模态对应的结构一阶位移场；与第δ阶密集屈曲模态对应的结构一阶位移场；

u_αβ,u_δγ,u_βγ,u_δα,u_γα,u_δβ分别为：与第α和β阶密集屈曲模态对应的结构二阶位移场；与第δ和γ阶密集屈曲模态对应的结构二阶位移场；与第β和γ阶密集屈曲模态对应的结构二阶位移场；与第δ和α阶密集屈曲模态对应的结构二阶位移场；与第γ和α阶密集屈曲模态对应的结构二阶位移场；与第δ和β阶密集屈曲模态对应的结构二阶位移场；

S( )为第一种结构高阶物理量S和括号内的结构二阶位移场的乘积形式；

为第一降阶模型物理量；

为第二降阶模型物理量；

为第三降阶模型物理量；

通过求解方程组(14)～(16)，得到结构一阶位移场u_α,u_β,u_γ,u_δ、结构二阶位移场u_αβ,u_δγ,u_βγ,u_δα,u_γα,u_δβ，以及第一降阶模型物理量分量第二降阶模型物理量分量和第三降阶模型物理量分量

步骤4.6，将第一降阶模型物理量分量第二降阶模型物理量分量和第三降阶模型物理量分量代入下式，得到曲纤维铺层平板结构的非线性有限元降阶模型：

其中：

ξ为降阶模型所对应的结构节点广义位移向量；

为结构载荷系数向量；

曲纤维铺层平板结构的非线性有限元降阶模型实质上是关于结构节点广义位移向量ξ的非线性方程组，其非线性有限元降阶模型的规模通常很小，且非线性有限元降阶模型的自由度数等于结构密集的屈曲模态个数u，通常有限元降阶模型的自由度数N_y<10。

步骤5，采用弧长求解技术求解步骤4得到的曲纤维铺层平板结构的非线性有限元降阶模型，得到曲纤维铺层平板结构的非线性力学承载响应曲线，所述曲纤维铺层平板结构的非线性力学承载响应曲线，为曲纤维铺层平板结构的非线性有限元降阶模型所对应的结构节点广义位移向量ξ随结构载荷系数向量的变化规律曲线；

具体的，采用与求解非线性有限元全阶模型类似的常规路径跟踪求解技术，如Newton-Raphson法或弧长法，即可获得非线性有限元降阶模型所对应的结构节点广义位移向量ξ随结构载荷系数向量的变化规律曲线。结构节点位移向量q可由结构节点广义位移向量ξ以及结构一阶位移场和结构二阶位移场组合得到，具体见步骤6。

步骤6，曲纤维铺层平板结构的结构节点位移向量q由结构节点广义位移向量ξ以及结构一阶位移场和结构二阶位移场组合得到，即：

q＝u_αξ_α+u_αβξ_αξ_β (18)

其中：ξ_α为与第α阶密集屈曲模态对应的结构节点广义位移向量ξ，通过结构节点广义位移向量ξ计算得到；

ξ_β为与第β阶密集屈曲模态对应的结构节点广义位移向量ξ，通过结构节点广义位移向量ξ计算得到；

步骤7，根据公式18，可将公式17换算为结构节点位移向量q随结构载荷系数向量的变化规律曲线；通过分析结构节点位移向量q随结构载荷系数向量的变化规律曲线，分析曲纤维铺层结构的力学性能。

下面列举一个具体实施例：

图3所示的数值试验为受面内轴压载荷作用的曲纤维铺层平板结构。曲纤维铺层平板结构的几何尺寸、加载和约束边界条件如图3所示。正交各项异性材料的属性为：E1＝47.0GPa、E2＝9.0GPa、v＝0.32、G12＝G13＝4.6GPa、G23＝3.12GPa。

本实施例根据用户要求的结构承载响应曲线长度，对该曲纤维铺层平板结构进行非线性有限元降阶建模，并开展非线性力学承载响应分析，所得到的加载点沿载荷方向位移随载荷变化曲线如图4所示。

本实施例具体步骤如下：

步骤1，依据曲纤维铺层平板结构中连续变角度曲纤维铺设路径的定义方式，即公式1，给出该曲纤维铺层平板结构的每个曲纤维铺层的信息。具体的，本例中，曲纤维铺层平板结构共含有4个曲纤维铺层，在建立坐标系时，η为0，因此，得到各个曲纤维铺层信息分别为：<0°|45°>/<0°|-45°>/<0°|-45°>/<0°|-45°>。

步骤2，构造曲纤维铺层平板结构的有限元模拟的四边形板壳单元，根据公式2计算得到单元应变向量ε，根据公式6得到单元本构关系。

本例中，该曲纤维铺层平板结构离散为1500个单元网格，总自由度数为3620。

步骤3，构造曲纤维铺层平板结构的非线性有限元全阶模型。

计算每个单元网格的单元应变能U，单元应变能U关于单元节点位移向量q^e计算至一阶导数后，得到单元节点内力向量f^e，然后对结构内部的所有单元进行组装，即可得到公式10所示的曲纤维铺层结构的非线性有限元全阶模型，该模型的总自由度数为3620。

步骤4，构造曲纤维铺层平板结构的非线性有限元降阶模型。

在步骤3的基础上，进一步对单元应变能U关于单元节点位移向量q^e计算至各高阶(二阶、三阶、四阶)导数，得到三种单元高阶物理量，即：S^e、和并在曲纤维铺层平板结构内部对所有单元的单元高阶物理量进行组装，得到曲纤维铺层结构的三种结构高阶物理量S、Q(q_α,q_β)和P(q_α,q_β,q_γ,q_δ)。

将计算获得的曲纤维铺层结构的三种结构高阶物理量S、Q(q_α,q_β)和P(q_α,q_β,q_γ,q_δ)带入到线性方程组(14)～(16)，求解线性方程组后即可得到非线性有限元降阶模型各物理量进而获得非线性有限元降阶模型的具体形式，即公式17。

步骤5，采用弧长求解技术求解曲纤维铺层结构的非线性有限元降阶模型，获得曲纤维铺层平板结构的非线性力学承载响应曲线，如图4。

由图4可知，本实施例通过求解曲纤维铺层结构的非线性有限元降阶模型与求解曲纤维铺层结构的非线性有限元全阶模型所得到的非线性力学承载响应曲线完全吻合，其跟踪计算精度达到要求。

为获得图4所示曲线，求解曲纤维铺层结构的非线性有限元全阶模型需要的CPU计算时间为56s，而本发明方法的计算时间仅为9s，计算时间仅约为传统方法的1/6，大幅提高了曲纤维铺层结构非线性力学分析的计算效率。

本发明与现有技术相比的优点在于：

曲纤维铺层结构在有限元离散后的模型规模较大(自由度数通常大于2000)，常规的非线性有限元分析技术要对结构非线性有限元全阶模型进行反复迭代求解，计算规模巨大，不利于曲纤维铺层结构的力学参数设计。为此，本发明旨在建立曲纤维铺层结构的线性有限元降阶模型，其自由度数通常小于10，采用对该小规模的降阶模型进行求解，可以快速精确地获得曲纤维铺层结构非线性承载的力学响应。

本发明具有对曲纤维铺层结构非线性有限元全阶模型进行精确降阶的特点，可准确经济地评估曲纤维铺层结构的真实承载能力。

综上所述，本发明提出的一种曲纤维铺层结构的力学性能分析方法，解决现有技术中曲纤维铺层结构离散后形成的有限元模型规模过于庞大，进而显著影响其非线性力学分析效率的问题；本发明能够实现高效准确的对曲纤维铺层结构的力学性能的非线性分析。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (2)

1.一种曲纤维铺层结构的力学性能分析方法，其特征在于，包括以下步骤：

步骤1，定义曲纤维铺层平板结构中连续变角度纤维的铺设路径，包括：

曲纤维铺层平板结构由多个曲纤维铺层叠加粘合而成；对于曲纤维铺层平板结构，以曲纤维铺层平板结构的中心为坐标原点θ₀，建立全局坐标系(x,y)；另外，定义参考坐标系(g,f)，参考坐标系(g,f)和全局坐标系(x,y)共相同的坐标原点θ₀，参考坐标系(g,f)为相对于全局坐标系逆时针旋转角度η得到的坐标系；曲纤维铺层平板结构的长度为d，则：在参考坐标系(g,f)下，对于通过坐标原点θ₀的曲纤维铺设路径，在距坐标原点θ₀的d/2处的曲纤维角度为θ₁，曲纤维铺设路径上任意一点的曲纤维角度表示为θ(g)，则得到参考坐标系(g,f)下的曲纤维铺设路径的描述方式为：

因此，曲纤维铺层平板结构的任意一个曲纤维铺层，由与通过坐标原点θ₀的曲纤维铺设路径平行的多个曲纤维铺设路径组成，因此，曲纤维铺层平板结构的任意一个曲纤维铺层，采用符号η<θ₀|θ₁>来描述；

步骤2，构造曲纤维铺层平板结构的有限元模拟的四边形板壳单元，并将所述曲纤维铺层平板结构离散为多个单元网格，计算每个单元网格的单元应变向量ε以及单元本构矩阵C_m，具体包括以下步骤：

步骤2.1，单元应变向量ε为单元线性应变向量ε_l与单元非线性应变向量ε_nl之和，即：

ε＝ε_l+ε_nl (2)

其中：单元线性应变向量ε_l为：

单元非线性应变向量ε_nl为：

其中：

为单元面内线性应变向量；

ε^b为单元面外线性应变向量；

为单元面内非线性应变向量；

B_l为单元线性几何插值矩阵，是常数矩阵；

B_nl(q^e)为关于单元节点位移向量q^e的单元非线性几何插值矩阵；

q^e为单元节点位移向量；

T代表矩阵的转置；

Z为单元非线性几何插值矩阵中的单元节点位移无关项，进一步表示为：

其中：

K_xx为单元非线性几何插值矩阵中节点位移无关项的xx向分量，是常数矩阵；

K_yy为单元非线性几何插值矩阵中节点位移无关项的yy向分量，是常数矩阵；

K_xy为单元非线性几何插值矩阵中节点位移无关项的xy向分量，是常数矩阵；

步骤2.2，根据参考坐标系(g,f)下的曲纤维铺设路径，构造曲纤维铺层平板结构的有限元模拟的四边形板壳单元，其单元本构关系为：

其中：

N为单元面内力向量；

M为单元弯矩向量；

A(θ)是单元面内刚度；

G(θ)是单元面外刚度；

D(θ)是单元面内/面外耦合刚度；

A(θ)、G(θ)和D(θ)的值受曲纤维铺设方式影响而随铺层面内曲纤维角度θ(g)的变化而变化；

ε^m是单元面内应变向量；

ε^b是单元面外应变向量；

根据单元本构关系，求得单元面内力向量N的表示式；

步骤2.3，获得单元本构矩阵C_m关于θ的表示式为：

步骤3，构造曲纤维铺层平板结构的非线性有限元全阶模型，包括：

步骤3.1，每个单元网格的单元应变能U关于q^e的表达式为：

其中：

i＝1,2,...,6，代表6个方向，分别为：x方向，y方向，z方向，xy方向，xz方向，yz方向；

j＝1,2,...,6；代表6个方向，分别为：x方向，y方向，z方向，xy方向，xz方向，yz方向；

A_s是单元网格的面积；

C_mij是单元本构矩阵C_m曲纤维角度中的第i行第j列的元素，单元本构矩阵C_m为6行6列矩阵；

ε_li是单元线性应变向量ε_l中的i向分量，即：根据公式3计算得到的单元线性应变向量ε_l中的i向分量，ε_li是关于q^e的表达式；

ε_lj是单元线性应变向量ε_l中的j向分量，即：根据公式3计算得到的单元线性应变向量ε_l中的j向分量，ε_lj是关于q^e的表达式；

ε_nli是单元非线性应变向量ε_nl中的i向分量，即：根据公式4计算得到的单元非线性应变向量ε_nl中的i向分量，ε_nli是关于q^e的表达式；

ε_nlj是单元非线性应变向量ε_nl中的j向分量，即：根据公式4计算得到的单元非线性应变向量ε_nl中的j向分量，ε_nlj是关于q^e的表达式；

步骤3.2，计算单元应变能U关于单元节点位移向量q^e的一阶导数，得到单元节点内力向量f^e；其中，单元节点是指每个单元网格的顶点；

其中：

为常规的单元线性刚度矩阵；

N_nl是单元面内力向量N的非线性部分；

T代表矩阵的转置；

步骤3.3，对曲纤维铺层平板结构的每个单元节点的单元节点内力向量f^e进行组装，获得曲纤维铺层平板结构的节点内力向量f(q)；

由此得到曲纤维铺层平板结构的非线性有限元全阶模型：

f(q)＝f_ext (10)

其中：f_ext为曲纤维铺层平板结构的节点外载荷向量；

步骤4，构造曲纤维铺层平板结构的非线性有限元降阶模型，包括：

步骤4.1，计算单元应变能U关于单元节点位移向量q^e的二阶导数，得到第一种单元高阶物理量S^e：

其中：

N_x为单元面内力向量N的x向分量；

N_y为单元面内力向量N的y向分量；

N_xy为单元面内力向量N的xy向分量；

步骤4.2，计算单元应变能U关于单元节点位移向量q^e的三阶导数，将三阶层数与任意的以及相乘，得到第二种单元高阶物理量

其中：

下标α＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

下标β＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

代表与第α阶密集屈曲模态对应的单元节点位移向量；

代表与第β阶密集屈曲模态对应的单元节点位移向量；其中，β≠α；

为关于的单元非线性几何插值矩阵；

为关于的单元非线性几何插值矩阵；

为关于的单元面内力向量；

为关于的单元面内力向量；

为关于和的单元面内力向量；

步骤4.3，计算单元应变能U关于单元节点位移向量q^e的四阶导数，将四阶导数与任意的以及相乘，得到第三种单元高阶物理量

其中：

下标α＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

下标β＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

下标γ＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

下标δ＝1,2,…,u，u为曲纤维铺层平板结构的前u阶密集屈曲模态；

其中：α≠β≠γ≠δ；

代表与第γ阶密集屈曲模态对应的单元节点位移向量；

代表与第δ阶密集屈曲模态对应的单元节点位移向量；

同时定义符号：

B_nl()为关于括号内参数的单元非线性几何插值矩阵；

N()为关于括号内参数的单元面内力向量；

步骤4.4，在计算得到曲纤维铺层平板结构的每个单元对应的第一种单元高阶物理量S^e、第二种单元高阶物理量和第三种单元高阶物理量后，对曲纤维铺层平板结构的各个单元进行组装，获得曲纤维铺层平板结构的整体的第一种结构高阶物理量S、第二种结构高阶物理量Q(q_α,q_β)和第三种结构高阶物理量P(q_α,q_β,q_γ,q_δ)；

步骤4.5，将第一种结构高阶物理量S、第二种结构高阶物理量Q(q_α,q_β)和第三种结构高阶物理量P(q_α,q_β,q_γ,q_δ)带入到以下三组线性方程组，

其中：

矩阵Y为摄动矩阵，其各列向量分别由曲纤维铺层平板结构的前u阶密集屈曲模态与结构几何刚度矩阵的乘积计算获得；

单位基向量E_α，为与第α阶密集屈曲模态对应的单位基向量，α＝1,2,...,u，其第α个分量为1，而其余各分量均为零；

u_α,u_β,u_γ,u_δ分别为：与第α阶密集屈曲模态对应的结构一阶位移场；与第β阶密集屈曲模态对应的结构一阶位移场；与第γ阶密集屈曲模态对应的结构一阶位移场；与第δ阶密集屈曲模态对应的结构一阶位移场；

u_αβ,u_δγ,u_βγ,u_δα,u_γα,u_δβ分别为：与第α和β阶密集屈曲模态对应的结构二阶位移场；与第δ和γ阶密集屈曲模态对应的结构二阶位移场；与第β和γ阶密集屈曲模态对应的结构二阶位移场；与第δ和α阶密集屈曲模态对应的结构二阶位移场；与第γ和α阶密集屈曲模态对应的结构二阶位移场；与第δ和β阶密集屈曲模态对应的结构二阶位移场；

S()为第一种结构高阶物理量S和括号内的结构二阶位移场的乘积形式；

为第一降阶模型物理量；

为第二降阶模型物理量；

为第三降阶模型物理量；

通过求解方程组(14)～(16)，得到结构一阶位移场u_α,u_β,u_γ,u_δ、结构二阶位移场u_αβ,u_δγ,u_βγ,u_δα,u_γα,u_δβ，以及第一降阶模型物理量分量第二降阶模型物理量分量和第三降阶模型物理量分量

步骤4.6，将第一降阶模型物理量分量第二降阶模型物理量分量和第三降阶模型物理量分量代入下式，得到曲纤维铺层平板结构的非线性有限元降阶模型：

其中：

ξ为降阶模型所对应的结构节点广义位移向量；

为结构载荷系数向量；

步骤5，求解步骤4得到的曲纤维铺层平板结构的非线性有限元降阶模型，得到曲纤维铺层平板结构的非线性力学承载响应曲线，所述曲纤维铺层平板结构的非线性力学承载响应曲线，为曲纤维铺层平板结构的非线性有限元降阶模型所对应的结构节点广义位移向量ξ随结构载荷系数向量的变化规律曲线；

步骤6，曲纤维铺层平板结构的结构节点位移向量q由结构节点广义位移向量ξ以及结构一阶位移场和结构二阶位移场组合得到，即：

q＝u_αξ_α+u_αβξ_αξ_β (18)

其中：ξ_α为与第α阶密集屈曲模态对应的结构节点广义位移向量ξ，通过结构节点广义位移向量ξ计算得到；

ξ_β为与第β阶密集屈曲模态对应的结构节点广义位移向量ξ，通过结构节点广义位移向量ξ计算得到；

步骤7，根据公式18，可将公式17换算为结构节点位移向量q随结构载荷系数向量的变化规律曲线；通过分析结构节点位移向量q随结构载荷系数向量的变化规律曲线，分析曲纤维铺层结构的力学性能。

2.根据权利要求1所述的一种曲纤维铺层结构的力学性能分析方法，其特征在于，步骤5中，采用弧长求解技术求解步骤4得到的曲纤维铺层平板结构的非线性有限元降阶模型。