CN115879346A - 基于改进型四节点逆有限元理论的结构应变场反演方法 - Google Patents

Abstract

本发明公开一种基于改进型四节点逆有限元理论的结构应变场反演方法，属于结构健康监测技术领域；包括以下步骤，步骤一：确定结构单元网格划分与传感器网络布局；步骤二：根据结构单元网格划分计算各个单元形函数矩阵；步骤三：求解结构表面应变；步骤四：构造误差函数；步骤五：重构结构应变场；本发明所述方法改进了四节点单元插值函数，通过对拟合得到的高阶位移场函数进行求导，得到板壳类结构的应变场，用较少的传感器同时得到结构的变形与应变响应且反演精度较高；本发明方法简单方便，实时性强。

Description

基于改进型四节点逆有限元理论的结构应变场反演方法

技术领域

本发明属于结构健康监测的应变场反演技术领域，具体提出了一种基于改进型四节点逆有限元理论的壁板结构应变场反演方法。

背景技术

近年来随着中国科技的飞速发展，航空航天事业技术得以不断发展，板壳等薄壁结构由于其重量轻、性能好等优点，在航空航天事业中得以广泛应用，例如在航空航天领域中，机体的承力构件大多为板壳结构，板壳结构的加工质量和装配精度决定了飞行器的质量，因此板壳结构受力变形正逐渐成为研究者们较为关注的话题。

飞行器在飞行服役过程中，常常遇到恶劣的工作环境，并承受着各种形式的载荷共同作用，因此很容易产生不同程度的、看不见的损伤，服役时间长了也会造成结构强度失效，比如符合材料板会产生内部纤维断裂等问题，金属板壳结构会产生局部変形等问题，使得飞机典型承力部件(机体、机翼)等结构的力学性能大幅下降，直接威胁飞行器的飞行安全。因此需要及时获取板壳结构表面应变、变形信息，为结构健康状态评估以及维修提供可靠依据。

目前，结构应变场的重构方法主要有两类：1)根据结构表面有限个离散点的应变，通过插值拟合得到整个结构表面的应变值；2)根据结构受力特点进行理论分析，根据边界条件和力学公式推导出结构的应变场。模态叠加原理首先要通过试验或有限元分析得到结构各阶应变模态和位移模态，结合应变模态和位移模态之间的导数关系，建立起被测结构应变场与位移场之间的关系，得到应变-位移转换矩阵，再输入由应变传感器测得的应变值，便可得到整个结构的应变场。

模态叠加原理进行有限元分析时，应变场的反演精度会受到有限元模型建模以及传感器数量的影响。截取不同的模态分析阶数会影响重构效果，当截取的模态数过少时，低阶模态对结构变形影响很大，缺少模态信息造成很大误差；当截取的模态阶数过多时，由于高阶模态计算结果会产生较大误差，因此也会影响到结构应变场的重构效果。研究发现，随着光纤FBG传感器数量的增加，得到的应变信息越多，重构误差越小。因此，针对模态叠加法的不足，在保证反演精度的同时减少FBG传感器的使用，使测试系统更加简洁。为此，本发明提出了基于改进型四节点逆有限元理论的结构应变场反演方法。

发明内容

技术问题：本发明所要解决的技术问题是提供一种用于工程板壳类结构的应变场反演方法。

该方法采用光纤FBG传感器测量结构不同点处的应变响应信息，将采集到的应变信息作为逆有限元输入计算得到结构的位移响应。由于FBG传感器是单点测量，因此得到的位移信息是离散不连续的，因此采用高阶函数来拟合重构得到的离散位移数据，得到结构表面挠度函数。最后根据欧拉-伯努利梁理论的应变与挠度之间的函数关系，对拟合得到的位移函数求二阶导数，便可得到结构表面应变函数。

技术方案：为解决上述技术问题，本发明的所采用的技术方案包括以下步骤：

1、一种基于改进型四节点逆有限元理论的壁板结构应变场反演方法，其特征在于，包括以下步骤：

步骤一：确定结构单元网格划分与传感器布局

选用四节点逆壳单元对结构进行离散，沿着结构长度方向离散为N个四节点逆壳单元，根据结构几何形状特点，曲率大的地方单元划分更密；每个单元节点按逆时针方向编号1～4，设置节点1所在点为单元局部坐标系原点O，沿结构长度方向为x轴方向，沿结构宽度方向为y轴方向，z轴与xy平面垂直；在每个逆壳单元上、下表面分别布置n组传感器，结构表面N×2n组传感器构成传感测量网络，监测区域为整个结构表面；传感器测点位置根据实际情况选择，每组设置三个传感器，编号为FBG1、FBG2和FBG3，按照0°、90°和45°方向粘贴在结构表面，测得结构该点x方向、y方向以及与x轴方向呈45°夹角方向的应变。

步骤二：根据壁板结构单元网格划分计算各个单元形函数矩阵

根据经典有限元理论分析，单元内任意一点的位移可由该单元内4个节点的挠度和转角线性叠加表示，表达式为：

式(1)中，u表示单元内一点沿x方向的位移，v表示单元内一点沿y方向的位移，ω表示单元内一点变形挠度，θ_x表示单元内一点沿x轴转角，θ_y表示单元内一点沿y轴转角，二维平面壳单元，在只考虑横向、纵向载荷的作用下不考虑沿z轴的转角θ_z；u_i、v_i、ω_i、θ_xi、θ_yi分别表示每个单元内各个节点的位移分量。N_i、L_i、M_i为四节点单元形函数，i(i＝1,2,3,4)为每个单元节点编号。

每个单元包含4个节点，每个节点的位移向量可以表示为：

因此，单元内4个节点的位移向量可以表示为：

实际应用当中，结构单元的划分不能保证全是规则的矩形单元，因此需要将一般形状的四边形单元等参成规则的矩形单元。引入基础坐标系s-t，s-t坐标系与x-y坐标系转换关系式为：

x2、x3、y3、y4为单元节点坐标。

则四节点单元形函数具体表达式为：

L₁＝y₁₄S₄-y₂₁S₁，M₁＝x₄₁S₄-x₁₂S₁

L₂＝y₂₁S₁-y₃₂S₂，M₂＝x₁₂S₁-x₂₃S₂

L₃＝y₃₂S₂-y₄₃S₃，M₃＝x₂₃S₂-x₃₄S₃

L₄＝y₄₃S₃-y₁₄S₄， M₄＝x₃₄S₃-x₄₁S₄ (6)

其中，

式(11)中：

x_ij＝x_i-x_j

y_ij＝y_i-y_j(i,j＝1,2,3,4) (8)

式(9)中，x_i、x_j、y_i、y_j为单元节点坐标，则每个单元的形函数矩阵N如下式所示：

N＝[N₁ N₂ N₃ N₄] (9)

其中N_i具体表达式为：

因此各个单元的形函数为6×24阶矩阵。

对于二维平面四节点单元，我们只考虑结构受横向载荷和纵向载荷的作用，即不考虑结构绕z轴方向的转角，因此单元形函数最后一行元素全为零。

结合式(3)和式(10)，可以求出结构各个单元内任一点位移，单元内任意一点位移为：

步骤三：求解结构表面应变

结构表面应变εe可以表示为面内拉压应变e(u^e)与弯曲应变k(u^e)的线性组合：

εe＝e(u^e)+z₀·k(u^e) (12)

其中，z₀表示壁板表面到中性层之间的距离。

逆有限元法需要结构表面三个应变分量：ε_x、ε_y、ε_45°，分别表示x方向、y方向以及45°方向的应变；结构表面拉压应变和弯曲应变用测得的三个方向应变表示为：

式中，“+”表示结构上表面应变，“-”表示结构效表面应变；γ_xy表示xy面内的切应变，可由测得的三个应变分量表示为：

γ_xy＝ε_x+ε_y-2·ε_45° (14)

式(13)用四节点单元形函数偏导矩阵B和单元节点位移u^e进一步表示为：

ε^e＝e(u^e)+z₀·k(u^e)＝B^mu^e+z₀·B^bu^e (15)

四节点单元形函数偏导矩阵B分为

和/>

具体表达式为：

四节点单元形函数的偏导矩阵可完整表示为：

步骤四：构造误差函数

将结构表面实验测得的应变值记为ε，包括拉压应变e^ε、弯曲应变k^ε和横向剪切应变g^ε，则结构实际应变值与理论应变值的误差函数为：

Φ(ε^e,ε)＝||e(u^e)-eε||²+||k(u^e)-k^ε||²+λ||g(u^e)-g^ε||² (18)

式中，g(u^e)表示结构表面理论横向剪切应变，通常为0；λ表示应变测量数据与理论结果之间相关程度的罚参数(0<λ<1)。

误差函数Φ对单元节点位移向量u^e求偏导并使其为0，求解微分方程得到误差函数的极小值，结果如式(20)所示：

通过计算得到如式(20)所示的矩阵方程：

k^eu^e＝f^e (20)

式中，k^e、f^e可由下式计算：

上式中，积分区域为单元面积A_e。

将式(21)和(22)代入式(20)，便可求出结构单元节点位移向量u^e，将u^e结果代回到式(11)便可求出结构内任一点位移分量。

步骤五：重构结构应变场

根据欧拉-伯努利理论，结构表面应变ε和变形挠度y之间的关系式为：

式中，z₀为结构表面到中性层的距离。

根据步骤四求得单元内任一点位移分量，选取各单元中心点处的挠度分量，共N个离散位移，用MATLAB拟合工具拟合出挠度与x方向坐标的关系式，对于平面四节点单元，位移函数选为：

y＝a₁xⁿ+a₂x^n-1+a₃x^n-2+...+a_nx+a_n+1(n＝3,4,5,...) (24)

将式(24)代入式(23)，便可得到壁板结构表面应变函数为：

ε(x)＝z₀[n(n-1)a₁x^n-2+(n-1)(n-2)a₂x^n-4+...+a_n-1] (25)

式(25)表示结构表面x方向应变分布规律，代入结构表面不同点坐标就可以重构出整个结构表面x方向应变场。

本发明的优点是：

本发明为一种基于改进型四节点逆有限元理论的结构应变场反演方法，该方法通过结构上、下表面布置的光纤FBG传感网络测得的离散应变响应来实现结构的位移场、应变场反演。本发明适用于板壳类结构的变形、应变监测和反演等工程应用领域，所具备的优点是：最少传感器情况下仅需采用4组12个光纤FBG传感器构成传感网络，相较于传统传感方式，具有线路布置简单、抗电磁干扰能力强等优点。其次，本发明通过用高阶位移场函数代替原来四边形单元的位移函数，将四边形常应变单元转换成应变值随测点坐标变化的单元，不需要考虑结构的材料属性，简化了位移、应变重构过程并且反演精度较高。

附图说明

图1是结构单元网格划分以及光纤FBG传感器布局；

图2是壁板结构挠度重构值与仿真值对比图；

图3是壁板结构应变重构值与仿真值对比图；

图4是改进型四节点逆有限元法流程图；

图5是重构得到结构的位移和应变结果。

具体实施方式

为了使本技术领域的人员能更好地理解本发明方案，下面将结合本发明实例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

如图1所示，本发明提供了一种基于逆有限元法和曲线拟合方法的结构变形、应变监测及重构方法，包括如下步骤：

S1、选用合适的逆壳单元对结构进行离散。

S2、基于中厚板理论计算每个单元的拉压应变、弯曲应变以及剪切应变。

S3、在逆壳单元上下表面选取应变测量点，并在应变测量点上粘贴光纤FBG传感器，实时测量应变，得到应变测量数据。

S4、基于得到的所述应变测量数据，计算每个逆壳单元的膜应变和曲率。

S5、基于最小二乘方法构造误差函数，对节点自由度进行求导，得到逆壳单元的伪刚度矩阵和伪载荷列阵并进行组装，最后结合边界条件计算得出结构的位移场。

S6、根据得到的结构离散的位移场，选择合适的高阶函数来拟合离散位移。

S7、对拟合得到的位移函数进行二次求导，再乘以壁板厚度的1/2就得到结构的应变函数。

具体实施时，作为本发明优选的实施方式，所述步骤S1中的结构包括所有板壳构成的复杂几何模型且均为等厚度板。如图1所示，为铝合金壁板单元划分、传感器网络布局示意图。具体实施时，做为本发明优选的实施方式，所述步骤S2的具体实现过程如下：

S21、计算每个逆壳单元的拉压应变，计算公式如下：

e(u^e)＝B^mu^e (1)

其中，u^e表示每个逆壳单元节点位移向量，B^m表示包含形函数导数的矩阵。

S22、计算每个逆壳单元的弯曲应变，计算公式如下：

k(u^e)＝B^bu^e (2)

其中，B^b表示包含形函数导数的矩阵。

S23、计算每个逆壳单元的剪切应变，计算公式如下：

g(u^e)＝B^su^e (3)

其中，B^s表示包含形函数导数的矩阵。

具体实施时，作为本发明优选的实施方式，所述步骤S3中的应变传感器实时测量的应变分别为：

其中，i表示第i个逆壳单元，n表示逆壳单元内的应变测量数目，ε_x表示x方向正应变，ε_y表示y方向的正应变，ε_45°表示结构表面45°方向的应变；“+”表示结构上表面应变，“-”表示结构下表面应变。

具体实施时，作为本发明优选的实施方式，所述应变传感器选用的是光纤FBG传感器。

具体实施时，作为本发明优选的实施方式，如图1所示，所述光纤FBG传感器沿单轴或应变花的形式设置在所述逆壳单元上下表面。

具体实施时，作为本发明优选的实施方式，所述步骤S4的具体实施过程如下：

S41、根据得到的应变实时测量数据，计算每个逆壳单元的拉压应变，计算公式如下：

S42、根据得到的应变实时测量数据，计算每个逆壳单元的弯曲应变，计算公式如下：

其中，z₀表示单元厚度的1/2。

具体实施时，作为本发明优选的实施方式，所述步骤S5的具体实施过程如下：

S51、对于每个逆壳单元，采用最小二乘误差函数，函数表达式为：

Φ(ε^e,ε＝‖e(u^e)-‖e^ε‖²+‖k(u^e)-k^ε‖²+λ‖g(u^e)-g^ε‖²(7)

式中，e(u^e)、k(u^e)、g(u^e)分别表示每个逆壳单元内的理论拉压应变、弯曲应变和剪切应变，e^ε、k^ε、g^ε分别表示每个逆壳单元实时测量的拉压应变、弯曲应变和剪切应变。λ表示应变测量数据与理论结果之间相关程度的罚参数(0<λ<1)。

S52、对结构的节点自由度进行求导，得到每个逆壳单元的伪刚度矩阵和伪载荷列阵。

k^eu^e=f^e(8)

上式中，k^e表示逆壳单元的伪刚度矩阵，f^e表示逆壳单元的伪载荷列阵，其中：

上式中，积分区域为单元面积A_e。

S53、将逆壳单元的伪刚度矩阵和伪载荷列阵按照有限元程序进行组装，得到离散结构的整体线性方程组，表达式如下：

K^eU^e＝F^e (10)

式中，K^e表示逆有限元方法整体的伪刚度矩阵，F^e表示逆有限元方法整体的伪载荷矩阵。S54、结合结构的边界条件，对整体伪刚度矩阵K^e进行修正，得到结构的位移场如下：

式中，

表示施加边界条件后的伪刚度矩阵，是正定矩阵。/>

具体实施时，作为本发明优选的实施方式，所述步骤S6的具体实施过程如下：

S61、根据欧拉-伯努利理论，结构表面应变ε和变形挠度y之间的关系式为：

式中，z₀为壁板表面到中性层的距离。

S62、根据步骤S54求得离散挠度分量，用MATLAB拟合工具拟合出挠度与x方向坐标的关系式。对于本发明实施例位移函数选为：

y＝a₁x³+a₂x²+a₃x+a₄ (13)

S63、将式(13)代入式(12)，实施例厚度为5mm，便可得到壁板结构表面应变场函数为：

ε(x)＝z₀(6a₁x+a₂x) (14)

重构得到结构的位移和应变结果，如图5所示。

以上各实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员依然可以对前述各实施例所记载的技术方案进行修改。

Claims (3)

1.一种基于改进型四节点逆有限元理论的结构应变场反演方法，其特征在于，包括以下步骤：

步骤一：确定结构单元网格划分与传感器布局

选用四节点逆壳单元对结构进行离散，沿着结构长度方向离散为N个四节点逆壳单元，根据结构几何形状特点，曲率大的地方单元划分更密；每个单元节点按逆时针方向编号为1～4，设置节点1所在点为单元局部坐标系原点O，沿结构长度方向为x轴方向，沿结构宽度方向为y轴方向，z轴与xy平面垂直；在每个逆壳单元上、下表面分别布置n组传感器，结构表面N×2n组传感器构成传感测量网络，监测区域为整个结构表面；传感器测点位置根据实际情况选择，每组设置三个传感器，编号为FBG1、FBG2和FBG3，按照0°、90°和45°方向粘贴在结构表面，测得结构该点x方向、y方向以及与x轴方向呈45°夹角方向的应变；

步骤二：根据结构单元网格划分计算各个单元形函数矩阵

根据经典有限元理论分析，单元内任意一点的位移可由该单元内4个节点的挠度和转角线性叠加表示，表达式为：

式(1)中，u表示单元内一点沿x方向的位移，v表示单元内一点沿y方向的位移，ω表示单元内一点变形挠度，θ_x表示单元内一点沿x轴转角，θ_y表示单元内一点沿y轴转角，二维平面壳单元，在只考虑横向、纵向载荷的作用下不考虑沿z轴的转角θ_z；u_i、v_i、ω_i、θ_xi、θ_yi、θ_zi分别表示每个单元内各个节点的位移分量；N_i、L_i、M_i为四节点单元形函数，具体计算方法见后续式(6)，i(i＝1,2,3,4)为每个单元节点编号；

每个单元包含4个节点，每个节点的位移向量可以表示为：

因此，单元内4个节点的位移向量可以表示为：

实际应用当中，结构单元的划分不能保证全是规则的矩形单元，因此需要将一般形状的四边形单元等参成规则的矩形单元；引入基础坐标系s-t，s-t坐标系与x-y坐标系转换关系式为：

x2、x3、y3、y4为单元节点坐标；

则四节点单元形函数具体表达式为：

L₁＝y₁₄S₄-y₂₁S₁，M₁＝x₄₁S₄-x₁₂S₁

L₂＝y₂₁S₁-y₃₂S₂，M₂＝x₁₂S₁-x₂₃S₂

L₃＝y₃₂S₂-y₄₃S₃，M₃＝x₂₃S₂-x₃₄S₃

L₄＝y₄₃S₃-y₁₄S₄，M₄＝x₃₄S₃-x₄₁S₄ (6)

其中，

式(11)中：

x_ij＝x_i-x_j

y_ij＝y_i-y_j(i,j＝1,2,3,4) (8)

式(9)中，x_i、x_j、y_i、y_j为单元节点坐标，则每个单元的形函数矩阵N如下式所示：

N＝[N₁ N₂ N₃ N₄] (9)

其中N_i具体表达式为：

因此各个单元的形函数为6×24阶矩阵；

对于二维平面四节点单元，我们只考虑结构受横向载荷和纵向载荷的作用，即不考虑结构绕z轴方向的转角，因此单元形函数最后一行元素全为零；

结合式(3)和式(9)，可以求出结构各个单元内任一点位移，单元内任意一点位移为：

步骤三：求解结构表面应变

结构表面应变ε^e可以表示为面内拉压应变e(u^e)与弯曲应变k(u^e)的线性组合：

ε^e＝e(u^e)+z₀·k(u^e) (12)

其中，z₀表示壁板表面到中性层之间的距离；

逆有限元法需要结构表面三个应变分量：ε_x、ε_y、ε_45°，分别表示x方向、y方向以及45°方向的应变；结构表面拉压应变和弯曲应变用测得的三个方向应变表示为：

式中，“+”表示结构上表面应变，“-”表示结构效表面应变；γ_xy表示xy面内的切应变，可由测得的三个应变分量表示为：

γ_xy＝ε_x+ε_y-2·ε_45° (14)

式(12)用四节点单元形函数偏导矩阵B和单元节点位移u^e进一步表示为：

ε^e＝e(u^e)+z₀·k(u^e)＝B^mu^e+z₀·B^bu^e (15)

四节点单元形函数偏导矩阵B分为

和/>

具体表达式为：

四节点单元形函数的偏导矩阵可完整表示为：

步骤四：构造误差函数

将结构表面实验测得的应变值记为ε，包括拉压应变e^ε、弯曲应变k^ε和横向剪切应变g^ε，则结构实际应变值与理论应变值的误差函数为：

Φ(ε^e,ε)＝||e(u^e)-e^ε||²+||k(u^e)-k^ε||²+λ||g(u^e)-g^ε||² (18)

式中，g(u^e)表示结构表面理论横向剪切应变，通常为0；λ表示应变测量数据与理论结果之间相关程度的罚参数(0<λ<1)；

误差函数Φ对单元节点位移向量u^e求偏导并使其为0，求解微分方程得到误差函数的极小值，结果如式(19)所示：

k^e表示结构的伪刚度矩阵，f^e表示结构的伪载荷列阵；

通过计算得到如式(20)所示的矩阵方程：

k^eu^e＝f^e (20)

式中，k^e、f^e可由下式计算：

上式中，积分区域为单元面积A_e；

将式(21)和(22)代入式(20)，便可求出结构单元节点位移向量u^e，将u^e结果代回到式(11)便可求出结构内任一点位移分量；

步骤五：重构结构应变场

根据欧拉-伯努利理论，结构表面应变ε和变形挠度y之间的关系式为：

式中，z₀为结构表面到中性层的距离；

根据步骤四求得单元内任一点位移分量，选取各单元中心点处的挠度分量，共N个离散位移，用MATLAB拟合工具拟合出挠度与x方向坐标的关系式，对于平面四节点单元，位移函数选为：

y＝a₁xⁿ+a₂x^n-1+a₃x^n-2+...+a_nx+a_n+1(n＝3,4,5,...) (24)

将式(24)代入式(23)，便可得到结构表面应变函数为：

ε(x)＝z₀[n(n-1)a₁x^n-2+(n-1)(n-2)a₂x^n-4+...+a_n-1] (25)

式(25)表示结构表面x方向应变分布规律，代入结构表面不同点坐标就可以重构出整个结构表面x方向应变场。

2.根据权利要求1所述的一种基于改进型四节点逆有限元理论的结构应变场反演方法，其特征在于：所述传感器为FBG光纤传感器。

3.根据权利要求1所述的一种基于改进型四节点逆有限元理论的结构应变场反演方法，其特征在于：所述传感器按照应变花形式设置在结构上下表面同一位置处。

Images