CN113468667A - 一种基于逆有限元与有限元方法的结构状态监测及载荷识别方法 - Google Patents

Abstract

本发明提供一种基于逆有限元与有限元方法的结构状态监测及载荷识别方法，包括：选用适配的逆壳单元对板壳结构进行离散；基于mindlin板理论计算每个逆壳单元的膜应变、弯曲应变以及剪切应变；在逆壳单元的上下表面选取应变测量点，并在应变测量点上粘贴应变传感器实时测量应变，得到应变测量数据；基于应变测量数据，计算每个逆壳单元的膜应变以及曲率；基于最小二乘方法构造泛函，对节点自由度求导，得到逆壳单元的类刚度矩阵和载荷矩阵并组装，赋予恰当的边界条件，计算结构的位移场；利用有限元方法构建的位移场与加载力向量之间的关系，并利用Tikhonov降低应变采集和位移重构产生的误差，实现对结构的状态信息的监测和载荷信息识别。

Description

一种基于逆有限元与有限元方法的结构状态监测及载荷识别 方法

技术领域

本发明涉及航空航天、车辆、建筑等领域的结构实时变形监测及载荷识别等技术领域，具体而言，尤其涉及一种基于逆有限元与有限元方法的结构状态监测及载荷识别方法。

背景技术

结构在长期的使用过程之中会不可避免的发生损伤。而损伤会引起结构的刚度退化甚至对结构的承载性能造成毁灭性的影响。在航空航天、车辆、建筑等领域，每年都要投入大量的精力来对结构进行检测和维修以保证结构的功能。因此，基于结构健康监测技术实时监测结构的状态对提升结构性能、保障结构安全、降低维修成本等起到巨大作用。结构健康监测技术是一种多学科的前沿技术，其中包含各种传感、重构方法。在各种状态传感方法中，逆有限元法(iFEM)充分考虑了边界条件和结构拓扑的复杂性，而且计算不需要载荷条件以及材料信息等众多先验信息，是一种强大的在线健康监测工具。结构的变形、应变等状态信息可以借助逆有限元方法实时准确得到。此外，与微元动态响应方法等损伤检测方法联合使用，可以进一步明确结构的损伤信息。而外载荷作为表征结构状态的重要一环，对结构损伤的扩展存在至关重要的影响，但是逆有限元方法在计算过程之中并不考虑外载荷的影响。

发明内容

根据上述提出的技术问题，提供一种基于逆有限元与有限元方法的结构状态监测及载荷识别方法。本发明通过实时测量板壳结构的表面应变，并采用逆有限元方法重构结构的变形等状态信息，再将重构得到的位移数据作为有限元方法的输入来反算结构的外载荷，在计算过程中，通过引入Tikhonov正则化方法来降低测量噪声以及逆有限元方法重构误差的影响，从而准确识别结构的外载荷的大小，进而实现对结构的状态监测以及载荷识别。

本发明采用的技术手段如下：

一种基于逆有限元与有限元方法的结构状态监测及载荷识别方法，包括如下步骤：

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

S2、基于mindlin板理论计算每个逆壳单元的膜应变、弯曲应变以及剪切应变；

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

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

S5、基于最小二乘方法构造泛函，对节点自由度求导，得到逆壳单元的类刚度矩阵和载荷矩阵并组装，最后赋予恰当的边界条件，计算结构的位移场；

S6、根据重构的位移信息构建经典有限元静力学平衡方程，推导得到重构的节点位移与未知的载荷大小之间的关系；

S7、利用Tikhonov正则化方法，降低应变采集和逆有限元重构过程中存在的误差，重构出结构外载荷的大小。

进一步地，所述步骤S1中的结构包括所有由板壳构成的复杂几何模型，且均为等厚度板。

进一步地，所述步骤S2的具体实现过程如下：

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

e(u^e)＝B^eu^e

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

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

k(u^e)＝B^ku^e

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

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

g(u^e)＝B^gu^e

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

进一步地，所述步骤S3中的应变传感器实时测量的应变分别为：

其中，i表示第i个逆壳单元，n表示逆壳单元内的应变测量点数目，ε表示1、2方向的正应变，γ表示1-2方向的切应变。

进一步地，所述应变传感器包括电阻应变片传感器、光纤布拉格光栅传感器以及分布式光纤传感器的一种或两种以上的结合。

进一步地，所述应变传感器沿单轴或应变花的形式设置在所述逆壳单元上下表面。

进一步地，所述步骤S4的具体实现过程如下：

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

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

其中，h表示单元厚度的一半。

进一步地，所述步骤S5的具体实现过程如下：

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

Φ_e(u^e)＝w_e||e(u^e)-e^ε||²+w_k||k(u^e)-k^ε||²+w_g||g(u^e)-g^ε||²

其中w_e、w_k和w_g分别表示与每个逆壳单元的膜应变、弯曲应变以及剪切应变相关的权重系数，若能通过应变测量值计算得到膜应变、弯曲应变以及剪切应变，则相应的权重系数为1；否则权重系数调整为一个小值；

S52、对结构的节点自由度取变分，得到每个逆壳单元的类刚度矩阵和载荷矩阵：

上式中，

表示逆壳单元的类刚度矩阵，

表示逆壳单元的类载荷矩阵；其中：

上式中，A^e表示逆壳单元的面积；

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

K_iU_i＝F_i

上式中，

其中，T^e表示坐标转换矩阵，K_i表示逆有限元方法全局的类刚度矩阵，F_i表示逆有限元方法全局的类载荷矩阵；

S54、施加合适的边界条件，计算得到结构的位移场，如下：

K_RiU_i＝F_i

其中，K_Ri表示施加边界条件后的类刚度矩阵，是正定阵。

进一步地，所述步骤S6的具体实现过程如下：

S61、根据经典有限元理论，将静态有限元求解方程改写为如下形式：

K_fU_f＝F_f

其中，K_f表示有限元方法中的刚度阵，U_f表示有限元方法中的位移向量，F_f表示有限元方法中的载荷向量；

S62、根据结构边界条件的不同，将上述方程表达为分块矩阵的形式，如下：

其中，下标1表示施加位移边界条件区域，下标2表示位移已知区域，下标3表示待求外载区域，下标4表示已知外载或者内部无外载区域；1、2区域的节点的位移是已知的，4区域的载荷是已知的，而测点选择时一般避开加载区域，因此区域2外载荷为0；

S63、以固定边界条件为例，即U₁＝0，将分块矩阵简化为如下形式：

S64、通过逆有限元计算得到的U₂，计算未知的作用力F_c：

其中，

进一步地，所述步骤S7的具体实现过程如下：

S71、根据Tikhonov正则化方法，将拟合误差表示为如下形式：

为了减小拟合误差

引入成本函数：

其中，λ表示正则化参数，H表示厄米特转置；为了使成本函数最小，对力向量F_c的一阶导数必须为零，计算得到：

S72、利用普通交叉验证法、广义交叉验证法或者L曲线法确定正则化参数λ的值，最后计算得到外载荷的大小。

较现有技术相比，本发明具有以下优点：

1、本发明提供的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，采用应变传感器测量得到结构表面应变，采用逆有限元方法实时重构结构的变形，再利用有限元方法重构结构的载荷大小。

2、本发明方法在精确重构结构状态的基础上，能够准确的指示载荷的大小，为后续结构的损伤扩展以及寿命评估提供依据。

3、与目前损伤识别方法相比，本发明方法具有精度高、准确、抗噪能力强等优点。

基于上述理由本发明可在航空航天、车辆、建筑等领域的结构实时的变形监测及载荷识别等领域广泛推广。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做以简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明方法流程图。

图2为本发明中铝合金板大小及加载位置示意图；

图3为本发明实施例提供的结构上下表面张贴传感器位置图

图中：1、加载1号位置；2、加载2号位置；3、结构上表面传感器布设位置；4、结构下表面传感器布设位置

具体实施方式

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

需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、***、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

如图1所示，本发明提供了一种基于逆有限元与有限元方法的结构状态监测及载荷识别方法，包括如下步骤：

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

S2、基于mindlin板理论计算每个逆壳单元的膜应变、弯曲应变以及剪切应变；

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

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

S5、基于最小二乘方法构造泛函，对节点自由度求导，得到逆壳单元的类刚度矩阵和载荷矩阵并组装，最后赋予恰当的边界条件，计算结构的位移场；

S6、根据重构的位移信息构建经典有限元静力学平衡方程，推导得到重构的节点位移与未知的载荷大小之间的关系；

S7、利用Tikhonov正则化方法，降低应变采集和逆有限元重构过程中存在的误差，重构出结构外载荷的大小。

具体实施时，作为本发明优选的实施方式，所述步骤S1中的结构包括所有由板壳构成的复杂几何模型，且均为等厚度板。如图2所示，为铝合金板大小及加载位置示意图。

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

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

e(u^e)＝B^eu^e

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

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

k(u^e)＝B^ku^e

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

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

g(u^e)＝B^gu^e

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

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

其中，i表示第i个逆壳单元，n表示逆壳单元内的应变测量点数目，ε表示1、2方向的正应变，γ表示1-2方向的切应变。

具体实施时，作为本发明优选的实施方式，所述应变传感器包括电阻应变片传感器、光纤布拉格光栅传感器以及分布式光纤传感器的一种或两种以上的结合。

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

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

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

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

其中，h表示单元厚度的一半。

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

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

Φ_e(u^e)＝w_e||e(u^e)-e^ε||²+w_k||k(u^e)-k^ε||²+w_g||g(u^e)-g^ε||²

其中w_e、w_k和w_g分别表示与每个逆壳单元的膜应变、弯曲应变以及剪切应变相关的权重系数，若能通过应变测量值计算得到膜应变、弯曲应变以及剪切应变，则相应的权重系数为1；否则权重系数调整为一个小值；

S52、对结构的节点自由度取变分，得到每个逆壳单元的类刚度矩阵和载荷矩阵：

上式中，

表示逆壳单元的类刚度矩阵，

表示逆壳单元的类载荷矩阵；其中：

上式中，A^e表示逆壳单元的面积；

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

K_iU_i＝F_i

上式中，

其中，T^e表示坐标转换矩阵，K_i表示逆有限元方法全局的类刚度矩阵，F_i表示逆有限元方法全局的类载荷矩阵；

S54、施加合适的边界条件，计算得到结构的位移场，如下：

K_RiU_i＝F_i

其中，K_Ri表示施加边界条件后的类刚度矩阵，是正定阵。

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

S61、根据经典有限元理论，将静态有限元求解方程改写为如下形式：

K_fU_f＝F_f

其中，K_f表示有限元方法中的刚度阵，U_f表示有限元方法中的位移向量，F_f表示有限元方法中的载荷向量；

S62、根据结构边界条件的不同，将上述方程表达为分块矩阵的形式，如下：

其中，下标1表示施加位移边界条件区域，下标2表示位移已知区域，下标3表示待求外载区域，下标4表示已知外载或者内部无外载区域；1、2区域的节点的位移是已知的，4区域的载荷是已知的，而测点选择时一般避开加载区域，因此区域2外载荷为0；

S63、以固定边界条件为例，即U₁＝0，将分块矩阵简化为如下形式：

S64、通过逆有限元计算得到的U₂，计算未知的作用力F_c：

其中，

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

S71、根据Tikhonov正则化方法，将拟合误差表示为如下形式：

为了减小拟合误差

引入成本函数：

其中，λ表示正则化参数，H表示厄米特转置；为了使成本函数最小，对力向量F_c的一阶导数必须为零，计算得到：

S72、利用普通交叉验证法、广义交叉验证法或者L曲线法确定正则化参数λ的值，最后计算得到外载荷的大小。重构得到的载荷与施加载荷结果，如表1所示：

表1计算结果

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (10)

1.一种基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，包括如下步骤：

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

S2、基于mindlin板理论计算每个逆壳单元的膜应变、弯曲应变以及剪切应变；

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

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

S5、基于最小二乘方法构造泛函，对节点自由度求导，得到逆壳单元的类刚度矩阵和载荷矩阵并组装，最后赋予恰当的边界条件，计算结构的的位移场；

S6、根据重构的位移信息构建经典有限元静力学平衡方程，推导得到重构的节点位移与未知的载荷大小之间的关系；

S7、利用Tikhonov正则化方法，降低应变采集和逆有限元重构过程中存在的误差，重构出结构外载荷的大小。

2.根据权利要求1所述的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，所述步骤S1中的结构包括所有由板壳构成的复杂几何模型，且均为等厚度板。

3.根据权利要求1所述的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，所述步骤S2的具体实现过程如下：

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

e(u^e)＝B^eu^e

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

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

k(u^e)＝B^ku^e

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

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

g(u^e)＝B^gu^e

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

4.根据权利要求1所述的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，所述步骤S3中的应变传感器实时测量的应变分别为：

其中，i表示第i个逆壳单元，n表示逆壳单元内的应变测量点数目，ε表示1、2方向的正应变，γ表示1-2方向的切应变。

5.根据权利要求4所述的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，所述应变传感器包括电阻应变片传感器、光纤布拉格光栅传感器以及分布式光纤传感器的一种或两种以上的结合。

6.根据权利要求5所述的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，所述应变传感器沿单轴或应变花的形式设置在所述逆壳单元上下表面。

7.根据权利要求1所述的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，所述步骤S4的具体实现过程如下：

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

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

其中，h表示单元厚度的一半。

8.根据权利要求1所述的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，所述步骤S5的具体实现过程如下：

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

Φ_e(u^e)＝w_e||e(u^e)-e^ε||²+w_k||k(u^e)-k^ε||²+w_g||g(u^e)-g^ε||²

其中w_e、w_k和w_g分别表示与每个逆壳单元的膜应变、弯曲应变以及剪切应变相关的权重系数，若能通过应变测量值计算得到膜应变、弯曲应变以及剪切应变，则相应的权重系数为1；否则权重系数调整为一个小值；

S52、对结构的节点自由度取变分，得到每个逆壳单元的类刚度矩阵和载荷矩阵：

上式中，

表示逆壳单元的类刚度矩阵，f_i ^e表示逆壳单元的类载荷矩阵；其中：

上式中，A^e表示逆壳单元的面积；

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

K_iU_i＝F_i

上式中，

其中，T^e表示坐标转换矩阵，K_i表示逆有限元方法全局的类刚度矩阵，F_i表示逆有限元方法全局的类载荷矩阵；

S54、施加合适的边界条件，计算得到的结构位移场，如下：

K_RiU_i＝F_i

其中，K_Ri表示施加边界条件后的类刚度矩阵，是正定阵。

9.根据权利要求1所述的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，所述步骤S6的具体实现过程如下：

S61、根据经典有限元理论，将静态有限元求解方程改写为如下形式：

K_fU_f＝F_f

其中，K_f表示有限元方法中的刚度阵，U_f表示有限元方法中的位移向量，F_f表示有限元方法中的载荷向量；

S62、根据结构边界条件的不同，将上述方程表达为分块矩阵的形式，如下：

其中，下标1表示施加位移边界条件区域，下标2表示位移已知区域，下标3表示待求外载区域，下标4表示已知外载或者内部无外载区域；1、2区域的节点的位移是已知的，4区域的载荷是已知的，而测点选择时一般避开加载区域，因此区域2外载荷为0；

S63、以固定边界条件为例，即U₁＝0，将分块矩阵简化为如下形式：

S64、通过逆有限元计算得到的U₂，计算未知的作用力F_c：

其中，

10.根据权利要求1所述的基于逆有限元与有限元方法的结构状态监测及载荷识别方法，其特征在于，所述步骤S7的具体实现过程如下：

S71、根据Tikhonov正则化方法，将拟合误差表示为如下形式：

为了减小拟合误差

引入成本函数：

其中，λ表示正则化参数，H表示厄米特转置；为了使成本函数最小，对力向量F_c的一阶导数必须为零，计算得到：

S72、利用普通交叉验证法、广义交叉验证法或者L曲线法确定正则化参数λ的值，最后计算得到外载荷的大小。

Images