基于单元信息矩阵组装的逆有限元变形重构方法

1.本发明属于航空航天结构健康监测技术领域，具体涉及一种基于单元信息矩阵组装的逆有限元变形重构方法。


背景技术：

2.空间大型阵列天线在轨运行中会出现复杂热变形，该变形与一般航空航天结构相比更加难以预测，原因在于热激励是由剧烈的高低温变化和非均匀的温度分布所引起，此外还需考虑与天线框架或卫星其他结构产生的复杂热力耦合情况。为确保空间大型天线能够正常工作，需要对因结构变形而产生的阵列幅相误差进行校正，对阵元波束进行反馈控制，而实时可靠的变形监测是校正和补偿的基础。
3.在基于应变信息的结构变形监测方法中，光纤布拉格光栅传感器具有重量轻、体积小、抗干扰能力强、电磁绝缘性能好等特点，适用于密布电磁敏感元件的在轨天线结构的变形监测。传感器获取实时应变，并结合先进的信号处理和变形重构算法，可实现航空航天结构在线实时感知。在已发展的基于应变信息的变形重构算法中，逆向有限元算法具有不需要结构材料和载荷特性且重构精度高等优势获得广泛关注。
4.然而，目前逆向有限元法在单元组装部分常采用的方法是首先分别确定每个单元的等效刚度矩阵和等效外力矩阵所有元素在全局等效刚度矩阵和等效外力矩阵中的行和列，再输送到对应位置。当单元数量大且对组装时间有要求时，该方法比较费时且占用较大数据内存，影响变形重构的实时性。


技术实现要素：

5.针对于上述现有技术的不足，本发明的目的在于提供一种基于单元信息矩阵组装的逆有限元变形重构方法，以解决现有技术中基于单元逐个组装的逆向有限元法重构实时性较差的问题。
6.为达到上述目的，本发明采用的技术方案如下：
7.本发明的一种基于单元信息矩阵组装的逆有限元变形重构方法，步骤如下：
8.1)构造子单元信息矩阵：将被监测结构离散为若干个带单元结点的逆向单元，分别对所有逆向单元和所有单元结点进行编号；构造每个逆向单元解析应变场和实测应变的最小二乘误差函数并进行寻优，得到每个逆向单元的等效刚度矩阵和等效外力矩阵，然后根据每个逆向单元的结点编号将得到的等效刚度矩阵中的各系数整理成子单元等效刚度信息矩阵，同理将等效外力矩阵中的各系数整理成子单元等效外力信息矩阵，形成n个子单元等效刚度信息矩阵和n个子单元等效外力信息矩阵；
9.2)整体组装并求解全局重构方程组：分别根据逆向单元编号依次排列n个子单元等效刚度信息矩阵和n个子单元等效外力信息矩阵；相邻逆向单元存在公共单元结点，将单元结点编号相同的信息进行合并，得到以稀疏矩阵形式存储的全单元等效刚度信息矩阵和全单元等效外力信息矩阵，再将稀疏矩阵转换为全矩阵即可构造全局等效刚度矩阵和全局
等效外力矩阵实现整体组装，最后建立全局重构方程组得到全局单元结点位移实现变形重构。
10.进一步地，所述步骤1)中的构造子单元信息矩阵过程具体包括：
11.11)将被监测结构离散为n个逆向单元，根据逆向单元类型共形成m个单元结点，分别对逆向单元和单元结点进行编号，构造每个逆向单元解析应变场和实测应变的最小二乘误差函数并进行寻优，得到当解析应变与实测应变误差最小时每个逆向单元的结点位移列阵ae所满足的方程组：
12.ken·ae
＝fen,n＝1,2...n
13.其中，ken和fen分别为第n个逆向单元的等效刚度矩阵和等效外力矩阵；
14.12)根据单元结点编号生成每个逆向单元结点自由度的全局编号；每个逆向单元具有多个单元结点，假设第n(n＝1
……
n)个逆向单元的第j个单元结点的编号为m(m＝1
……
m)，每个单元结点自由度为dof，则该单元结点的所有自由度编号为dof
·
(m-1) 1、dof
·
(m-1) 2
……
dof
·
m，遍历该逆向单元的所有单元结点得到第n(n＝1
……
n)个逆向单元全部单元结点的自由度编号，作为第n个逆向单元的等效刚度矩阵ken和等效外力矩阵fen的元素行和列的位置索引；
15.13)将等效刚度矩阵ken所有元素对应的下标存成行索引向量i
ke-n
和列索引向量j
ke-n
，以及等效刚度矩阵ken数据向量v
ke-n
，[i
ke-n j
ke-n v
ke-n
]即为子单元等效刚度信息矩阵；同理将等效外力矩阵fen所有元素对应的下标存成行索引向量i
fe-n
和列索引向量j
fe-n
，以及等效外力矩阵fen数据向量v
fe-n
，[i
fe-n j
fe-n v
fe-n
]即为子单元等效外力信息矩阵；遍历所有逆向单元，得到n个子单元等效刚度信息矩阵和n个子单元等效外力信息矩阵。
[0016]
进一步地，所述步骤2)中的整体组装并求解全局重构方程组过程具体包括：
[0017]
21)按照逆向单元编号依次排列所有子单元等效刚度信息矩阵，得到全单元等效刚度信息拼接矩阵[i
ke j
ke v
ke
]，同理按照逆向单元编号依次排列所有子单元等效外力信息矩阵，得到全单元等效外力信息拼接矩阵[i
fe j
fe v
fe
]：
[0018][0019][0020]
22)将全单元等效刚度信息拼接矩阵[i
ke j
ke v
ke
]中具有相同行索引i
ke
和列索引j
ke
位置的v
ke
上的元素进行累加，构造稀疏矩阵形式的全单元等效刚度信息矩阵s
ke
，同理将全单元等效外力信息拼接矩阵[i
fe j
fe v
fe
]中具有相同行索引i
fe
和列索引j
fe
位置的v
fe
上的元素进行累加，构造稀疏矩阵形式的全单元等效外力信息矩阵s
fe
；
[0021]
23)全单元等效刚度信息矩阵s
ke
和全单元等效外力信息矩阵s
fe
为稀疏矩阵，实现了对全单元等效刚度信息和全单元等效外力信息的压缩存储；通过补充零元素将稀疏矩阵转变为全矩阵，得到全局等效刚度矩阵k和全局等效外力矩阵f，并建立全局重构方程组：
[0022]k·
u＝f
[0023]
补充结构约束条件解方程得到全局单元结点位移u，进一步采用形函数插值法得到被监测结构内任一点变形情况，从而实现变形重构。
[0024]
本发明的有益效果：
[0025]
本发明通过构造每个逆向单元的子单元信息矩阵，并利用稀疏矩阵实现所有逆向单元的子信息矩阵的组装和存储得到全单元信息矩阵，再将稀疏矩阵形式的全单元信息矩阵转换为全矩阵的方法，得到全局等效刚度矩阵和全局等效外力矩阵，最终建立全局变形重构方程组，求解该方程组得到全局单元结点位移，进一步可以采用形函数插值法得到结构内任一点变形情况从而实现变形重构，数据存储量小，重构速度快，能够用于实时变形监测，在满足航空航天结构变形监测需求的同时大幅提高了重构速度。
附图说明
[0026]
图1为本发明方法的流程图。
[0027]
图2为实施例采用的复合材料蜂窝夹芯板结构及夹具位置示意图。
[0028]
图3为逆向单元划分示意图。
[0029]
图4a为作为逆有限元法的算法输入的蜂窝夹芯板上表面沿x轴方向正应变的仿真应变场示意图。
[0030]
图4b为作为逆有限元法的算法输入的蜂窝夹芯板上表面沿y轴方向正应变的仿真应变场示意图。
[0031]
图4c为作为逆有限元法的算法输入的蜂窝夹芯板上表面在xy平面内剪应变的仿真应变场示意图。
[0032]
图5a为作为逆有限元法的算法输入的蜂窝夹芯板下表面沿x轴方向正应变的仿真应变场示意图。
[0033]
图5b为作为逆有限元法的算法输入的蜂窝夹芯板下表面沿y轴方向正应变的仿真应变场示意图。
[0034]
图5c为作为逆有限元法的算法输入的蜂窝夹芯板下表面在xy平面内剪应变的仿真应变场示意图。
[0035]
图6a为蜂窝夹芯板在垂直于板面方向的仿真位移场示意图。
[0036]
图6b为采用逆有限元法进行蜂窝夹芯板变形重构得到的垂直于板面方向的重构位移场示意图。
[0037]
图6c为重构位移场与仿真位移场的误差分布示意图。
具体实施方式
[0038]
为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。
[0039]
参照图1所示，本发明的一种基于单元信息矩阵组装的逆向有限元变形重构方法，步骤如下：
[0040]
1)构造子单元信息矩阵：将被监测结构离散为若干个带单元结点的逆向单元，分别对所有逆向单元和所有单元结点进行编号；构造每个逆向单元解析应变场和实测应变的最小二乘误差函数并进行寻优，得到每个逆向单元的等效刚度矩阵和等效外力矩阵，然后
根据每个逆向单元的结点编号将得到的等效刚度矩阵中的各系数整理成子单元等效刚度信息矩阵，同理将等效外力矩阵中的各系数整理成子单元等效外力信息矩阵，形成n个子单元等效刚度信息矩阵和n个子单元等效外力信息矩阵；
[0041]
具体地，所述步骤1)中的构造子单元信息矩阵过程具体包括：
[0042]
11)将被监测结构离散为n个逆向单元，根据逆向单元类型共形成m个单元结点，分别对逆向单元和单元结点进行编号，构造每个逆向单元解析应变场和实测应变的最小二乘误差函数并进行寻优，得到当解析应变与实测应变误差最小时每个逆向单元的结点位移列阵ae所满足的方程组：
[0043]
ken·ae
＝fen,n＝1,2...n
[0044]
其中，ken和fen分别为第n个逆向单元的等效刚度矩阵和等效外力矩阵；
[0045]
12)根据单元结点编号生成每个逆向单元结点自由度的全局编号；每个逆向单元具有多个单元结点，假设第n(n＝1
……
n)个逆向单元的第j个单元结点的编号为m(m＝1
……
m)，每个单元结点自由度为dof，则该单元结点的所有自由度编号为dof
·
(m-1) 1、dof
·
(m-1) 2
……
dof
·
m，遍历该逆向单元的所有单元结点得到第n(n＝1
……
n)个逆向单元全部单元结点的自由度编号，作为第n个逆向单元的等效刚度矩阵ken和等效外力矩阵fen的元素行和列的位置索引；
[0046]
13)将等效刚度矩阵ken所有元素对应的下标存成行索引向量i
ke-n
和列索引向量j
ke-n
，以及等效刚度矩阵ken数据向量v
ke-n
，[i
ke-n j
ke-n v
ke-n
]即为子单元等效刚度信息矩阵；同理将等效外力矩阵fen所有元素对应的下标存成行索引向量i
fe-n
和列索引向量j
fe-n
，以及等效外力矩阵fen数据向量v
fe-n
，[i
fe-n j
fe-n v
fe-n
]即为子单元等效外力信息矩阵；遍历所有逆向单元，得到n个子单元等效刚度信息矩阵和n个子单元等效外力信息矩阵。
[0047]
2)整体组装并求解全局重构方程组：分别根据逆向单元编号依次排列n个子单元等效刚度信息矩阵和n个子单元等效外力信息矩阵；相邻逆向单元存在公共单元结点，将单元结点编号相同的信息进行合并，得到以稀疏矩阵形式存储的全单元等效刚度信息矩阵和全单元等效外力信息矩阵实现整体组装，再将稀疏矩阵转换为全矩阵即可构造全局等效刚度矩阵和全局等效外力矩阵，最后建立全局重构方程组得到全局单元结点位移实现变形重构；
[0048]
具体地，所述步骤2)中的整体组装并求解全局重构方程组过程具体包括：
[0049]
21)按照逆向单元编号依次排列所有子单元等效刚度信息矩阵，得到全单元等效刚度信息拼接矩阵[i
ke j
ke v
ke
]，同理按照逆向单元编号依次排列所有子单元等效外力信息矩阵，得到全单元等效外力信息拼接矩阵[i
fe j
fe v
fe
]：
[0050][0051][0052]
22)将全单元等效刚度信息拼接矩阵[i
ke j
ke v
ke
]中具有相同行索引i
ke
和列索引jke
位置的v
ke
上的元素进行累加，构造稀疏矩阵形式的全单元等效刚度信息矩阵s
ke
，同理将全单元等效外力信息拼接矩阵[i
fe j
fe v
fe
]中具有相同行索引i
fe
和列索引j
fe
位置的v
fe
上的元素进行累加，构造稀疏矩阵形式的全单元等效外力信息矩阵s
fe
；
[0053]
23)全单元等效刚度信息矩阵s
ke
和全单元等效外力信息矩阵s
fe
为稀疏矩阵，实现了对全单元等效刚度信息和全单元等效外力信息的压缩存储；通过补充零元素将稀疏矩阵转变为全矩阵，得到全局等效刚度矩阵k和全局等效外力矩阵f，并建立全局重构方程组：
[0054]k·
u＝f
[0055]
补充结构约束条件解方程得到全局单元结点位移u，进一步采用形函数插值法得到被监测结构内任一点变形情况，从而实现变形重构。
[0056]
本实施例的变形重构对象为复合材料蜂窝夹芯板结构，如图2所示，尺寸为3125mm
×
1500mm
×
20mm(长
×
宽
×
厚)，下表面设置有加热元件，在加热过程中厚度方向存在温度梯度，使得蜂窝夹芯板产生翘曲变形。将该蜂窝夹芯板按照仿真位移的梯度1mm划分成29
×
23个四节点逆向壳单元，如图3所示，每个单元有4个结点，每个结点有6个自由度。蜂窝夹芯板上表面和下表面正应变和剪应变的应变场分布如图4a-图4c和图5a-图5c所示，提取每个逆向单元的仿真应变与解析应变场构造最小二乘误差函数，得到每个逆向单元的等效刚度矩阵ke和等效外力矩阵fe，其中，等效刚度矩阵ke的矩阵大小为24
×
24，等效外力矩阵fe的矩阵大小为24
×
1。
[0057]
根据单元结点编号生成每个逆向单元结点自由度的全局编号；假设第n个逆向单元的四个单元结点编号为i、j、k、s，则该逆向单元结点自由度的全局编号为(6
×
i-5
…6×
i，6
×
j-5
…6×
j，6
×
k-5
…6×
k，6
×
s-5
…6×
s)，分别将等效刚度矩阵ken所有元素对应的下标存成行索引向量i
ke-n
和列索引向量j
ke-n
，以及数据向量v
ke-n
，此处行索引向量、列索引向量以及数据向量的大小均为576
×
1，得到子单元等效刚度信息矩阵[i
ke-n j
ke-n v
ke-n
]；同理将等效外力矩阵fen所有元素对应的下标存成行索引向量i
fe-n
和列索引向量j
fe-n
，以及数据向量v
fe-n
，此处行索引向量、列索引向量以及数据向量的大小均为24
×
1，得到子单元等效外力信息矩阵[i
fe-n j
fe-n v
fe-n
]。
[0058]
按照逆向单元编号依次排列全部667个逆向单元的子单元等效刚度信息矩阵得到全单元等效刚度信息拼接矩阵[i
ke j
ke v
ke
]，该矩阵大小为384192
×
3；依次排列全部667个逆向单元的子单元等效外力信息矩阵得到全单元等效外力信息拼接矩阵[i
fe j
fe v
fe
]，该矩阵大小为16008
×
3。
[0059]
将具有相同行索引i
ke
和列索引j
ke
位置处的数据向量v
ke
的元素进行累加，得到全单元等效刚度信息矩阵s
ke
，同理将具有相同行索引i
fe
和列索引j
fe
位置处的数据向量v
fe
的元素进行累加得到全单元等效外力信息矩阵s
fe
。
[0060]ske
和s
fe
分别是以稀疏矩阵形式实现了对所有逆向单元的等效刚度矩阵和等效外力矩阵中各元素及其对应位置的压缩存储和单元整体组装。
[0061]
在s
ke
和s
fe
空余位置补充零元素实现稀疏矩阵到全矩阵的转换，分别得到大小为4320
×
4320的全局等效刚度矩阵k和大小为4320
×
1的全局等效外力矩阵f并建立全局重构方程组：
[0062]k·
u＝f
[0063]
补充结构约束条件并求解该方程组即可得到全局单元结点位移u，进一步可以采
用形函数插值法得到结构内任一点变形情况，从而实现变形重构，重构结果如图6a-图6c所示。
[0064]
比较基于单元信息矩阵组装与单元逐个叠加的重构精度和组装速度，结果表明两种单元组装方式计算的重构位移在0.01mm精度上结果相等，最大误差位置重构误差均为8.86％，基于单元信息矩阵组装时间为28.43ms，而单元逐个叠加组装时间为52.94ms。因此基于单元信息矩阵整体组装能够提速超过40％，有效改善变形重构实时性。
[0065]
本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。