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

权利要求书

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方向应变场。