船舶力学

单元数量/ N×N	参考解	最大总位移/m		最大误差/ %
6×6	2.65×10⁻³	2.293×10⁻³	2.068×10⁻³		12.98	21.52
8×8	2.401×10⁻³	2.172×10⁻³		8.88	17.57
10×10	2.479×10⁻³	2.260×10⁻³		5.92	14.23
12×12	2.532×10⁻³	2.339×10⁻³		3.91	11.23
14×14	2.567×10⁻³	2.400×10⁻³		2.58	8.92
16×16	2.580×10⁻³	2.445×10⁻³		1.98	7.21
18×18	2.583×10⁻³	2.465×10⁻³		1.97	6.45
20×20	2.583×10⁻³	2.480×10⁻³		1.97	5.88
22×22	2.583×10⁻³	2.482×10⁻³		1.97	5.80
24×24	2.583×10⁻³	2.482×10⁻³		1.97	5.80

单元数量/ N×N	参考解	最大总位移/m		最大误差/ %
6×6	2.65×10⁻³	2.293×10⁻³	2.068×10⁻³		12.98	21.52
8×8	2.401×10⁻³	2.172×10⁻³		8.88	17.57
10×10	2.479×10⁻³	2.260×10⁻³		5.92	14.23
12×12	2.532×10⁻³	2.339×10⁻³		3.91	11.23
14×14	2.567×10⁻³	2.400×10⁻³		2.58	8.92
16×16	2.580×10⁻³	2.445×10⁻³		1.98	7.21
18×18	2.583×10⁻³	2.465×10⁻³		1.97	6.45
20×20	2.583×10⁻³	2.480×10⁻³		1.97	5.88
22×22	2.583×10⁻³	2.482×10⁻³		1.97	5.80
24×24	2.583×10⁻³	2.482×10⁻³		1.97	5.80

数据来源	最大总位移/m	最大误差/%
参考解	7.362×10⁻⁴	0
改进方法	7.232×10⁻⁴	1.77
iQS4法	6.952×10⁻⁴	5.57

数据来源	最大总位移/m	最大误差/%
参考解	7.362×10⁻⁴	0
改进方法	7.232×10⁻⁴	1.77
iQS4法	6.952×10⁻⁴	5.57

高精度四边形逆壳单元的算法研究

PDF下载

闫宏生 ¹^,² , 唐蒋拼 ²

船舶力学 | 结构力学 2025,29(12): 1930-1942

收起

船舶力学 |结构力学 2025 , 29 (12) : 1930 -1942

高精度四边形逆壳单元的算法研究

全屏

[Author(id=1281326978879959715, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, orderNo=0, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=hongsheng@tju.edu.cn, emailSecond=null, emailThird=null, correspondingAuthor=1, authorType=1, ext={EN=AuthorExt(id=1281326978951262886, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, authorId=1281326978879959715, language=EN, stringName=Hong-sheng YAN, firstName=Hong-sheng, middleName=null, lastName=YAN, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=¹^,², address=^1.State Key Laboratory of Water Conservancy Engineering Simulation and Safety, Tianjin University, Tianjin 300350
^2.School of Architecture and Civil Engineering, Tianjin University, 300350, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1281326979018371751, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, authorId=1281326978879959715, language=CN, stringName=闫宏生, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=¹^,², address=^1.天津大学水利工程仿真与安全国家重点实验室，天津 300072
^2.天津大学建筑工程学院，天津 300350, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1281326978682827420, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, xref=1., ext=[AuthorCompanyExt(id=1281326978691216029, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, companyId=1281326978682827420, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^1.State Key Laboratory of Water Conservancy Engineering Simulation and Safety, Tianjin University, Tianjin 300350), AuthorCompanyExt(id=1281326978699604638, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, companyId=1281326978682827420, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^1.天津大学水利工程仿真与安全国家重点实验室，天津 300072)]), AuthorCompany(id=1281326978800267935, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, xref=2., ext=[AuthorCompanyExt(id=1281326978808656544, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, companyId=1281326978800267935, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^2.School of Architecture and Civil Engineering, Tianjin University, 300350), AuthorCompanyExt(id=1281326978812850849, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, companyId=1281326978800267935, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^2.天津大学建筑工程学院，天津 300350)])]), Author(id=1281326979077092009, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, orderNo=1, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=strugglesad@tju.edu.cn, emailSecond=null, emailThird=null, correspondingAuthor=0, authorType=1, ext={EN=AuthorExt(id=1281326979148395179, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, authorId=1281326979077092009, language=EN, stringName=Jiang-pin TANG, firstName=Jiang-pin, middleName=null, lastName=TANG, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=², address=^2.School of Architecture and Civil Engineering, Tianjin University, 300350, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1281326979207115436, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, authorId=1281326979077092009, language=CN, stringName=唐蒋拼, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=², address=^2.天津大学建筑工程学院，天津 300350, bio={"content":"

唐蒋拼（2000–），男，硕士，E-mail：strugglesad@tju.edu.cn

"}, bioImg=null, bioContent=

唐蒋拼（2000–），男，硕士，E-mail：strugglesad@tju.edu.cn

, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1281326978800267935, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, xref=2., ext=[AuthorCompanyExt(id=1281326978808656544, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, companyId=1281326978800267935, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^2.School of Architecture and Civil Engineering, Tianjin University, 300350), AuthorCompanyExt(id=1281326978812850849, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281326957195407835, companyId=1281326978800267935, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^2.天津大学建筑工程学院，天津 300350)])])]

闫宏生¹^,², 唐蒋拼²

作者信息

^1.天津大学水利工程仿真与安全国家重点实验室，天津 300072

^2.天津大学建筑工程学院，天津 300350

通讯作者:

闫宏生（1973–），男，博士，副教授，通讯作者，E-mail：hongsheng@tju.edu.cn

作者简介:

唐蒋拼（2000–），男，硕士，E-mail：strugglesad@tju.edu.cn

Algorithm for high-precision quadrilateral inverse shell elements

Hong-sheng YAN¹^,², Jiang-pin TANG²

Affiliations

^1.State Key Laboratory of Water Conservancy Engineering Simulation and Safety, Tianjin University, Tianjin 300350

^2.School of Architecture and Civil Engineering, Tianjin University, 300350

出版时间: 2025-12-15 doi: 10.3969/j.issn.1007-7294.2025.12.010

文章导航

摘要

收起

船舶结构会因极端或周期性载荷、腐蚀和侵蚀而失效，因此有必要对船舶进行结构健康监测（SHM），以确保船舶结构的安全性、可靠性和完整性，以避免重大事故的发生。位移重构是结构健康监测的主要目标之一，逆有限元法（iFEM）是一种用于板壳结构全场位移重构的结构健康监测方法。现有的逆有限元算法应用于大曲率曲面壳体，会出现单元翘曲问题而影响拟合精度。本文提出了一种改进逆有限元算法，该算法不改变单元本身的位移–应变关系，在不增加网格和节点数量的前提下，通过在单元高斯积分点上建立切平面坐标系，对单元刚度积分进行改进，提高逆壳单元的计算精度。数值算例结果表明，改进的逆壳单元能有效降低逆有限元重构后的位移误差。

关键词

逆有限元法 / 结构健康监测 / 位移重构 / 板壳结构

Abstract

收起

Ship structures can fail due to extreme or cyclic loading, corrosion and erosion, so it is necessary to perform structural health monitoring (SHM) on ships to ensure the safety, reliability and integrity of ship structures to avoid major accidents. Displacement reconstruction is one of the main objectives of structural health monitoring. The inverse finite element method (iFEM) is a structural health monitoring method for real-time reconstruction of full-field displacement in plate and shell structures. In the case of shells with large curvature, existing iFEM algorithm will suffer from the warping problem, which affects the fitting accuracy. In this paper, an improved iFEM algorithm is proposed, which does not change the displacement-strain relationship of the element itself. The new algorithm improves the calculation accuracy of the shell element by estblishing tangent plane coordinate system at the element’s Gaussian integration points and improving the element stiffness integration, without increasing the mesh and node quantities. The result of numerical examples shows that the improved inverse shell element can effectively reduce the displacement error after inverse finite element reconstruction.

Key words

inverse finite element method / structural health monitoring / displacement reconstruction / plate and shell structures

引用本文

闫宏生, 唐蒋拼. 高精度四边形逆壳单元的算法研究. 船舶力学, 2025 , 29 (12) : 1930 -1942 . DOI: 10.3969/j.issn.1007-7294.2025.12.010

Hong-sheng YAN, Jiang-pin TANG. Algorithm for high-precision quadrilateral inverse shell elements[J]. Journal of Ship Mechanics, 2025 , 29 (12) : 1930 -1942 . DOI: 10.3969/j.issn.1007-7294.2025.12.010

正文

收起

0　引言

收起

如今，随着海洋资源的不断开发，船舶在油气开采、远洋捕捞和科学研究等诸多方面发挥着重要作用。然而，船舶经常航行在恶劣的海洋环境中，其结构可能会因极端或周期性负载、腐蚀和侵蚀而失效，而结构失效可能导致重大事故，造成船员或乘客丧生、海洋环境污染以及高昂的维护或维修成本。因此，有必要利用结构健康监测（SHM）系统对船舶进行实时监测，以避免重大事故的发生。结构健康监测（SHM）是一个跨学科的工程，通过将传感系统集成到结构上，实时处理从传感系统收集到的数据，提供整体或局部结构的实时信息。逆有限元法（iFEM）是结构健康监测方法的一种，其在实际工程应用中，通过采集结构上离散点的应变数据实时重构结构整体的位移场，实现从应变到位移的转换，从而得到结构的整体变形信息，进而对结构状态进行实时的安全评估。

Tessler等^[1−2]开发了逆有限元法（iFEM）用于工程结构的位移和应力监测。与其他已开发的SHM方法相比，iFEM 普遍适用于任何类型的结构拓扑和边界条件，因为使用逆梁、板和壳单元可以对物理域进行有效离散^[3−5]。由于 iFEM 的加权最小二乘函数仅通过应变–位移关系来确定，因此 iFEM 可以在不预先了解材料特性和加载信息的情况下重构结构的变形形状。另外，与其他的 SHM 系统不同，通过 iFEM 获得的位移重构结果的稳定性和准确性，与结构的材料属性以及作用在结构上的任何类型的静态或动态负载无关。一旦获得了结构的变形形状，就可以根据应变–位移关系计算出全场位移，进而评估结构的全场应力。研究^[6-8]表明，iFEM 是一种精确、稳定、高效的变形和应力监测算法。Tessler等^[9]采用一阶运动学假设和剪切变形理论，开发了一种三节点逆壳单元（iMIN3）。Quach等^[10]通过实验测量实时应变数据，证实了 iMIN3 单元的稳定性。Kefal等^[11−12]利用一阶运动学假设和横向剪切变形理论开发了四节点四边形逆壳单元 iQS4 和八节点四边形逆壳单元 iCS8，对海洋结构进行了位移和应力监测^[13−16]。国内也有学者开展了逆有限元法的研究。周富^[17]通过机床立柱变形重构实验，将逆有限元法重构的机床立柱变形位移场与接触式位移传感器测量值相比较，检验逆有限元法对机床立柱变形场监测的正确性。赵雪晨^[18]设计了一种基于逆有限元法的损伤判别方法。付书山等^[19]针对机翼蒙皮变形动态监测需求，研究基于逆有限元法和应变信息的蒙皮变形重构方法。徐浩等^[20]基于蜂窝板模型以及等效的夹层板模型，利用不同的应变传感器布设方案，研究了最佳的蜂窝板状态重构策略。张科等^[21]以真实复杂变形机翼的主承力结构——鱼骨为研究对象，提出了一种基于逆有限元法与位移分段叠加思想结合的变形监测方法。闫宏生等^[22]对结构变形重构技术进行了研究，利用基于最小二乘变分原理的逆有限元法对多种案例的船体加筋板结构进行位移重构。

2016年，Kefal^[11]首次将四节点四边形逆壳单元（iQS4）应用到薄壁圆筒等三维结构上，初步展示了 iQS4 单元在三维结构中的应用潜力。然而，使用传统的iQS4单元分析曲面结构，并不能保证单元的所有节点都在同一个平面上，这就造成了单元的翘曲问题，不能准确模拟曲面几何形状。为此，需要将网格加密以逼近真实的曲面形状。特别是对于大曲率的曲面壳体结构，需要较密的网格才能达到良好的精度，导致计算量增大。

各种大型海洋结构，包括船体、海工结构和潜艇等，都是具有复杂曲面的海洋结构。为了解决 iQS4 单元应用于以上结构时可能出现的翘曲问题，以提高 iQS4 单元在这些结构中的计算精度，本文提出一种可以提高大曲率曲面壳体位移重构精度的改进逆有限元算法，在单元的高斯积分点处建立切平面坐标系，然后对单元刚度积分进行改进，最后结合算例对改进方法的准确性和有效性进行了证明。

1　iQS4单元的改进

收起

iQS4 单元是基于Mindlin–Ressiner板理论，通过加权最小二乘公式进行推导的。本文提出一种改进逆有限元算法，开发了适用于大曲率曲面的改进四节点四边形逆壳单元。板壳结构的逆有限元方法通过将结构先离散化再进行叠加来重构结构变形，本节首先对 iQS4 单元的单元局部坐标系（如图1所示）进行改进，然后对iQS4单元的刚度积分进行改进，最后给出改进的算法流程。

1.1　改进iQS4局部坐标系

应用iQS4单元对曲面结构进行网格划分的过程中，特别是对于大曲率的曲面壳体结构，单元的四个节点往往不在同一个平面上。对于这样的翘曲单元，利用传统方法建立的单元局部坐标系不能很好地描述单元的位置信息，特别是在单元刚度矩阵由局部坐标系向整体坐标系转换的时候，将会造成较大的误差。如图2所示，由于单元的四个节点不在同一个平面上，单元的局部坐标系所在的平面和实际的单元面并不能完全重合，两者偏差的程度可用两者在节点处相差的距离d来表示。

针对传统iQS4单元的局限性，本文引入一种单元切平面建立局部坐标系的方法，如图3所示。建立单元切平面局部坐标系的关键是确定逆壳单元在原点处的切平面，首先求出四节点四边形单元的双线性等参形函数 $ N(\xi ,\eta ) $

。然后求出参数坐标系下的坐标 $ {x}_{1}(\xi ,\eta ) $

、$ {x}_{2}(\xi ,\eta ) $

、$ {x}_{3}(\xi ,\eta ) $

，可以采用整体坐标系下单元各个节点的坐标表示，即

$ {x}_{{i}}\left(\xi ,\eta \right)=\sum _{\mathrm{a}=1}^{\mathrm{M}}N\left(\xi ,\eta \right){x}_{{i}}^{\mathrm{a}}\left({i}=1\sim 3\right)$

(1)

式中：$ N(\xi ,\eta ) $

为单元的形函数，M为单元节点个数，对于iQS4单元M=4。

建立参数坐标系之后，下一步即可建立单元切平面局部坐标系，三个坐标轴的方向向量$ {\boldsymbol{r}}_{\mathrm{\xi }} $

、$ {\boldsymbol{r}}_{\mathrm{\eta }} $

和 $ \boldsymbol{n} $

的表达式为^[23]

$ {\boldsymbol{r}}_{\mathrm{\xi }}=\frac{\partial {x}_{1}(\mathrm{\xi },\mathrm{\eta })}{\partial \mathrm{\xi }}{\boldsymbol{i}}+\frac{\partial {x}_{2}(\mathrm{\xi },\mathrm{\eta })}{\partial \mathrm{\xi }}{\boldsymbol{j}}+\frac{\partial {x}_{3}(\mathrm{\xi },\mathrm{\eta })}{\partial \mathrm{\xi }}{\boldsymbol{k}} $

(2a)

$ {\boldsymbol{r}}_{\mathrm{\eta }}=\frac{\partial {x}_{1}(\xi ,\eta )}{\partial \eta }{\boldsymbol{i}}+\frac{\partial {x}_{2}(\xi ,\eta )}{\partial \eta }{\boldsymbol{j}}+\frac{\partial {x}_{3}(\xi ,\eta )}{\partial \eta }{\boldsymbol{k}} $

(2b)

$ {\boldsymbol{n}}={\boldsymbol{r}}_{\mathrm{\xi }}\times {\boldsymbol{r}}_{\mathrm{\eta }} $

(2c)

根据文献[25]，可以得到切平面局部坐标系的坐标转换矩阵 $ {\boldsymbol{T}}^{{\mathrm{e}}} $

，从而得到从局部坐标系到整体坐标系的转换矩阵 $ \boldsymbol{T} $

，对于iQS4单元 $ \boldsymbol{T} $

为24×24的矩阵。通过转换矩阵 $ \boldsymbol{T} $

将单元刚度矩阵 $ {\boldsymbol{k}}^\text{e} $

转换到整体坐标系得到 $ \boldsymbol{K} $

，即

$ K={\boldsymbol{T}}^{\mathrm{T}}\left(\sum _{{i}=1}^{{N}}\sum _{{j}=1}^{{N}}{{\text{φ}}}_{{i}{j}}\right)\boldsymbol{T} $

(3)

式中：N表示高斯积分点数，$ {{\text{φ}}}_{{i}{j}} $

为单元刚度矩阵在髙斯点处的计算值。

前文提到iQS4单元的局部坐标系所在的平面和实际的单元面并不能完全重合，其偏差的程度可用两者在节点处相差的距离d来表示，要进一步提高计算精度，就要减小距离d。

根据文献[24]，高斯积分求和过程只在高斯积分点取值，在每个高斯积分点建立不同的局部坐标系，可以有效保证每次求和过程和刚度转换过程的精度。由于每个局部坐标系都建立在单元面高斯点处的切平面上，多个局部坐标系的应用便能够有效反映曲面壳体结构，减小距离d，如图4所示。

在每个高斯点处的局部坐标系下，依据式（3）便可得到从该局部坐标系到整体坐标系的坐标转换矩阵 $ \boldsymbol{T} $

。在积分求和之前，将每个高斯积分点处的刚度值转换到整体坐标系。根据局部坐标系所在的切平面位置不同，实现对转换矩阵的及时修正，从而能够更加有效模拟曲面壳体的真实形状。改进算法的整体坐标系下，逆壳单元刚度矩阵的表达式应修正为

$ K=\sum _{{i}=1}^{{N}}\sum _{{j}=1}^{{N}}{\boldsymbol{T}}_{{i}{j}}^\text{T}{{\text{φ}}}_{{i}{j}}{\boldsymbol{T}}_{{i}{j}} $

(4)

式中：N表示高斯积分点数；$ {\boldsymbol{T}}_{{i}{j}} $

表示从建立在不同高斯积分点处的局部坐标系到整体坐标系的单元刚度转换矩阵；$ {{\text{φ}}}_{{i}{j}} $

为单元刚度矩阵在髙斯点处的计算值。

1.2　改进iQS4单元刚度积分

根据Kefal^[11]等推导的公式，iQS4单元的单元方程如下

$ {\boldsymbol{k}}^{\mathrm{e}}{\boldsymbol{u}}^{\mathrm{e}}={\boldsymbol{f}}^{\mathrm{e}} $

(5)

式中： $ {\boldsymbol{k}}^{\mathrm{e}} $

为单元刚度矩阵， $ {\boldsymbol{f}}^{\mathrm{e}} $

为单元载荷列阵，均由式（6）计算得到。

$ {\boldsymbol{k}}^{\mathrm{e}}=\iint \left({w}^{m}\right({\boldsymbol{B}}^{m}{)}^\text{T}{\boldsymbol{B}}^{m}+{w}^{b}\left(2h{)}^{2}\right({\boldsymbol{B}}^{b}{)}^\text{T}{\boldsymbol{B}}^{b}+{w}^{g}\left({\boldsymbol{B}}^{g}{)}^\text{T}{\boldsymbol{B}}^{g}\right){\mathrm{d}}x{\mathrm{d}}y $

(6a)

$ {\boldsymbol{f}}^{\mathrm{e}}=\iint \left({w}^{m}\right({\boldsymbol{B}}^{m}{)}^\text{T}{\boldsymbol{e}}_{i}^{\boldsymbol{\varepsilon }}+{w}^{b}\left(2h{)}^{2}\right({\boldsymbol{B}}^{b}{)}^\text{T}{\boldsymbol{k}}_{i}^{\boldsymbol{\varepsilon }}+{w}^{g}\left({\boldsymbol{B}}^{g}{)}^\text{T}{\boldsymbol{g}}_{i}^{\boldsymbol{\varepsilon }}\right){\mathrm{d}}x{\mathrm{d}}y $

(6b)

式中：$ {\boldsymbol{e}}_{i}^{\boldsymbol{\varepsilon }} $

、$ {\boldsymbol{k}}_{i}^{\boldsymbol{\varepsilon }} $

和$ {\boldsymbol{g}}_{i}^{\boldsymbol{\varepsilon }} $

为应变传感器测量到的表面应变，w^m、w^b和w^g为加权常数，$ {\boldsymbol{B}}^{m} $

、$ {\boldsymbol{B}}^{b} $

和 $ {\boldsymbol{B}}^{g} $

为单元应变矩阵。要得到单元刚度矩阵 $ {\boldsymbol{k}}^{\mathrm{e}} $

，需要先计算单元形函数对局部坐标的导数。传统逆有限元方法需要推导出偏导数 $ \dfrac{\partial }{\partial x},\dfrac{\partial }{\partial y} $

与 $ \dfrac{\partial }{\partial \xi },\dfrac{\partial }{\partial \eta } $

之间的变换关系，并计算雅可比行列式 $ \left|\boldsymbol{J}\right| $

，根据复合函数求导法则有

$ \left[\begin{array}{c}\dfrac{\partial }{\partial {x}}\\ \dfrac{\partial }{\partial {y}}\end{array}\right]={\left[\boldsymbol{J}\right]}^{-1}\left[\begin{array}{c}\dfrac{\partial }{\partial \mathrm{\xi }}\\ \dfrac{\partial }{\partial \mathrm{\eta }}\end{array}\right] $

(7)

$ \left[\boldsymbol{J}\right]=\left[\begin{array}{cc}\dfrac{\partial x}{\partial \mathrm{\xi }}& \dfrac{\partial y}{\partial \mathrm{\xi }}\\ \dfrac{\partial x}{\partial \mathrm{\eta }}& \dfrac{\partial y}{\partial \mathrm{\eta }}\end{array}\right]=\left[\begin{array}{cc}\displaystyle\sum _{{i}=1}^{4}\dfrac{\partial {N}_{{i}}}{\partial \mathrm{\xi }}{x}_{{i}}& \displaystyle\sum _{{i}=1}^{4}\dfrac{\partial {N}_{{i}}}{\partial \mathrm{\xi }}{y}_{{i}}\\ \displaystyle\sum _{{i}=1}^{4}\dfrac{\partial {N}_{{i}}}{\partial \mathrm{\eta }}{x}_{{i}}& \displaystyle\sum _{{i}=1}^{4}\dfrac{\partial {N}_{i}}{\partial \mathrm{\eta }}{y}_{{i}}\end{array}\right] $

(8)

式中：$ {x}_{{i}} $

与$ {y}_{{i}} $

为将单元的整体坐标转换到二维平面局部坐标系后得到的坐标值。iQS4单元的误差源自于局部坐标系所在的平面和实际单元面在节点处相差的距离，为了减小上述误差，本文引入基于三维整体坐标系下的形函数对局部坐标的导数 $ {\boldsymbol{N}}_{{i},{x}'} $

、$ {\boldsymbol{N}}_{{i},{y}'} $

和雅可比行列式 $ \left|\boldsymbol{J}\right| $

的求法，不需要将单元的整体坐标转换到二维平面局部坐标系即可求得，更适用于三维壳体曲面的积分，具有更高的精度。由Gao和Davies^[25]推导的公式可知，从整体坐标系到参数坐标系积分转换的雅可比行列式为

$\left|\boldsymbol{J}\right|=\left|{\boldsymbol{r}}_{\xi }\times {\boldsymbol{r}}_{\eta }\right| $

(9)

式中：$ {\boldsymbol{r}}_{\mathrm{\xi }} $

和$ {\boldsymbol{r}}_{\mathrm{\eta }} $

为参数坐标轴$ (\xi ,\eta ) $

分别在原点相切的向量。由复合函数求导法则可知

$\binom{\boldsymbol{N}_{{i}, {x}^{\prime}}}{\boldsymbol{N}_{{i}, {y}^{\prime}}}=\binom{\boldsymbol{N}_{{i}, {x}_1}}{\boldsymbol{N}_{{i}, {x}_2}}=\left[\begin{array}{cc}\dfrac{\partial \xi}{\partial \boldsymbol{x}_1^{\prime}} & \dfrac{\partial \eta}{\partial \boldsymbol{x}_1^{\prime}} \\\dfrac{\partial \xi}{\partial \boldsymbol{x}_2^{\prime}} & \dfrac{\partial \eta}{\partial \boldsymbol{x}_2^{\prime}}\end{array}\right]\binom{\boldsymbol{N}_{{i}, \xi}}{\boldsymbol{N}_{{i}, \eta}}$

(10)

由Lachat和Becher^[26-27]推导的公式如下

$ {\boldsymbol{x}}'_{1}=\left|{m}_{1}\right|\xi +\left|{m}_{2}\right|\eta {\mathrm{cos}}\theta $

(11a)

$ {\boldsymbol{x}}'_{2}=\left|{m}_{2}\right|\eta {\mathrm{sin}}\theta $

(11b)

式中：$ \theta $

为图3中所示$ {x}'_{1} $

轴和$ \eta $

轴的夹角，且有

$ {\mathrm{cos}}\theta =\frac{1}{\left|{m}_{1}\right|\left|{m}_{2}\right|}\frac{\partial {x}_{{i}}}{\partial \xi }\frac{\partial {x}_{{i}}}{\partial \eta }(i=1\sim 3) $

(12)

整理可得

$\frac{\partial \xi }{\partial {\boldsymbol{x}}'_{1}}=\frac{1}{\left|{m}_{1}\right|}$

(13a)

$ \frac{\partial \xi }{\partial {\boldsymbol{x}}'_{2}}=\frac{-{\mathrm{cos}}\theta }{\left|{m}_{1}\right|{\mathrm{sin}}\theta } $

(13b)

$ \frac{\partial \eta }{\partial {\boldsymbol{x}}'_{1}}=0 $

(13c)

$ \frac{\partial \eta }{\partial {\boldsymbol{x}}'_{2}}=\frac{1}{\left|{m}_{1}\right|{\mathrm{sin}}\theta } $

(13d)

由式（10）可以得到单元刚度矩阵 $ {\boldsymbol{B}}^{m} $

、$ {\boldsymbol{B}}^{b} $

和 $ {\boldsymbol{B}}^{g} $

，按前文中组装整体刚度矩阵的方法合成总体刚度矩阵和总体逆载荷列阵，进而可得到改进方法的结构整体位移。

1.3　改进方法的算法流程

本文提出的改进方法实际上是一种新的逆有限元算法，其有效性可以通过数值计算进行验证。为了编写计算机程序，我们提出了改进方法的算法流程，具体实现过程如图5所示。其中应变数据通过有限元法（FEM）计算获得。只需要保证有限元模型的网格足够密集，以获得近似平滑的应变场，单元高斯积分点上的应变即可由此获得。这种方法可以确定任何高斯积分点的应变，并将其用作测量应变。

2　数值算例分析

收起

在本节中将通过两个数值算例验证改进方法在大曲率曲面壳体位移重构中的性能，所有算例均使用 ABAQUS 软件进行计算。首先提取表面上的应变，将其作为结构变形后应变传感系统测得的应变，然后将测量到的应变输入到使用改进方法编写的程序中，以实现结构的全场位移重构，所有算例均采用2×2的完全积分。为了验证改进方法的准确性，我们将把改进方法与 iQS4 单元的计算结果进行比较，并将ABAQUS的计算结果作为参考解。另外，为了确保网格的细化不会过度干扰计算结果，本文引入了网格独立性验证。

2.1　大曲率曲面壳体位移重构

如图6所示大曲率双曲曲壳，曲壳底边半径为3 m，从底边到顶点高为12 m，壳体厚度为20 mm，材料为线弹性的、均匀且各向同性的，弹性模量为2×10¹¹ Pa，泊松比为0.3。曲壳的顶点A为刚性固定，曲壳B、D两点施加一对指向圆心的径向集中力P₁=1000 N；曲壳C、E两点施加一对远离圆心的径向集中力P₂=1000 N。对曲壳进行不规则网格划分，AB等四条侧边都被划分为N个单元，底边BCDE也被划分为N个单元。为便于说明，我们将单元数量表示为 N×N。

翘曲度是用来衡量一个单元偏离平面程度的指标，翘曲发生在单元面的节点不共面的时候，若将四边形单元沿其对角线分成两个三角形，这两个三角形单元的法向量之间的夹角（锐角）就定义为该四边形单元的翘曲度。一般来说，翘曲度小于 5°时，可以认为较为理想。

图7（a）展示了曲壳在不同单元数量下的平均翘曲度，平均翘曲度定义为所有单元翘曲度的平均值。图7（b）展示了曲壳在不同单元数量下的最大翘曲度，最大翘曲度定义为所有单元翘曲度的最大值。可以发现，由于改进算法较好地解决了单元翘曲问题，改进方法的平均翘曲度和最大翘曲度都明显小于iQS4单元。此外，即使在网格较为稀疏时，改进方法的平均翘曲度已经小于5°，而iQS4单元则需要更为密集的网格，才能保证平均翘曲度小于5°。当有足够的网格密度时，改进算法的最大翘曲度能够保证小于5°，而iQS4单元的最大翘曲度在本文工况范围内始终不能满足小于5°的要求。

定义总位移 $ {U}_{\mathrm{M}\mathrm{a}\mathrm{g}\mathrm{n}{i}\mathrm{t}\mathrm{u}\mathrm{d}\mathrm{e}}=\sqrt{{U}_{1}^{2}+{U}_{2}^{2}+{U}_{3}^{2}} $

，其中$ {U}_{1} $

、$ {U}_{2} $

和$ {U}_{3} $

分别为x、y和${{z}} $

方向上的位移，ABAQUS计算得到的总位移为参考总位移，相似地，逆有限元法重构得到的总位移即为重构总位移。最大总位移为所有节点总位移中的最大值，最大百分比误差为所有百分比误差中的最大值。

为了直观地说明单元数量对重构精度的影响，图8（a）展示了曲壳在不同单元数量下的最大总位移，图8（b）展示了曲壳在不同单元数量下，改进方法和iQS4单元相对于参考解的最大百分比误差。从图8（a）中可以明显看出，改进方法在单元数量为16×16时收敛，而 iQS4单元则在单元数量为22×22时收敛。从图8（b）可以发现，无论单元数量多少，改进方法的最大误差都明显小于iQS4。

为了更精确地说明重构精度的变化，表1展示了曲壳在不同单元数量下的最大总位移和最大百分比误差的具体数据。从表中可以发现，当改进方法收敛，即单元数量为16×16时，改进方法的最大总位移为2.580×10⁻³ m，而 iQS4 的对应值为 2.445×10⁻³ m。参考解的最大总位移为2.653×10⁻³ m，相较于参考解，改进方法的最大误差为1.97%， iQS4的最大误差则为 7.21%。对比分析表明，改进方法不仅比 iQS4 收敛得更快，而且重构精度远高于 iQS4。

图9 展示了参考解的总位移云图，图10展示了单元数量为6×6、12×12、18×18和24×24时重构解的总位移云图。可以看出，随着单元数量的增加，改进方法的总位移云图与参考解越来越相似。此外，相较于iQS4单元法，改进方法的重构精度明显更高。当单元数量为12×12时，改进方法的总位移云图已十分光滑，和参考解基本一致，而此时iQS4单元由于存在翘曲，总位移云图仍较为粗糙。

为进一步地定量分析曲壳的整体变形情况，须对曲壳进行节点位移分析。由于曲壳的几何对称性，选取图11中的13个节点进行分析，其中A点$z $

轴坐标为12 m，C点$z $

轴坐标为0 m。画出单元数量分别为6×6、12×12、18×18和24×24时参考解和重构解的$z $

轴坐标–总位移曲线，如图12所示。从图12中可以明显看出，相对于iQS4单元法，改进方法的总位移曲线更接近参考解。此外，随着单元数量的增加，改进方法的总位移曲线与参考解越来越吻合。当单元数量为24×24时，iQS4单元发生了单元翘曲，尽管单元数量已不少，但总位移重构结果仍然存在一定误差，而改进方法由于较好地解决了单元翘曲问题，其总位移曲线已经和参考解基本一致。

2.2　典型船体结构的三维位移重构

为了验证改进方法在船舶结构中的实用性，本文选取Kefal^[12]提出的加筋圆柱壳进行位移重构。实际上，加筋圆柱壳可视为船舶结构（如加筋的舭部）的代表。因此，本部分的分析有助于提高对曲面海洋结构的实时SHM的能力。如图13所示，圆柱壳的均匀厚度为 30 mm，长度、宽度和半径均为6 m。圆柱壳内表面加固了11个等间距的纵骨，高度为0.2 m，厚度均为20 mm。此外，圆柱壳还由另外三根等间距的横桁横向支撑，这些横桁呈L形，高度为0.4 m，翼缘尺寸为0.2 m，厚度为25 mm。整体结构的材料弹性模量为210 GPa，泊松比为0.3。圆柱壳外表面承受1 MPa的均匀法向压力，压力方向指向圆柱壳内表面。圆柱壳的四个边缘均为刚性固定。考虑到工程实用性，需要尽可能减少逆单元数量，因此将加筋圆柱壳划分为246个逆单元，如图14所示。

图15和图16 分别展示了加筋圆柱壳内外侧的总位移云图，表2展示了加筋圆柱壳的最大总位移和最大百分比误差。可以看出，参考解的最大位移为7.362×10⁻⁴ m，改进方法的最大位移为7.232×10⁻⁴ m，百分比误差为1.77%。而iQS4法的最大位移为6.952×10⁻⁴ m，百分比误差为5.57%。显然，改进方法的重构结果更为精确，而且改进方法的位移云图变化趋势与参考解更为相似。

由于加筋圆柱壳是一个包含多个构件的复杂结构，需要进一步分析加筋圆柱壳各个构件的变形情况。一般来说，工程上最为关心的是关键构件的变形情况，关键构件由于受载较大或厚度较小等原因，往往是最容易发生结构失效的部位，因此我们只取变形较大的三个关键构件作为研究对象，如图17所示，包括圆柱壳面、横桁和纵骨。

图18展示了三个关键构件的具体位置和节点选取。三个关键构件都选取了13个节点。其中圆柱壳面和横桁的所有节点从下至上编号为1–13，纵骨的所有节点从左至右编号为1–13。画出参考解和重构解的节点序号–总位移曲线，如图19~21所示。

图19、图20和图21分别展示了三个关键构件的节点总位移和节点百分比误差。可以看出，对于三个关键构件，改进方法与参考解十分接近，百分比误差在2%左右。另外，在位移较大的节点（如节点7）处，iQS4法与参考解差别较大，百分比误差远大于改进方法。显然，对于关键构件而言，改进方法的重构结果更为精确。

3　结论

收起

本文提出了一种高精度逆有限元法，以提高大曲率曲面壳体位移重构的精度。通过建立单元切平面坐标系并在高斯积分点上建立单元局部坐标系，然后改进单元刚度积分，改进了iQS4单元的算法流程，数值实例证明，改进的逆壳单元降低了由逆有限元法重构引起的位移误差。具体结论如下：

（1）对于大曲率曲面壳体，iQS4单元的四个节点往往不在同一个平面上。对于这样的翘曲单元，在高斯积分点建立切平面局部坐标系，可以实现对转换矩阵的及时修正，能够更加有效地模拟曲面壳体的真实形状，从而有效解决单元翘曲问题。

（2）对于iQS4单元，为了得到单元刚度矩阵，需要计算二维局部坐标系下的形函数导数，当曲面曲率较大时，得到的单元刚度矩阵误差也较大，而基于三维空间给出的形函数导数更适合三维壳体曲面的积分，精度更高。

（3）对于大曲率曲面壳体，本文提出的改进方法在计算精度和收敛速度方面均优于iQS4单元法，能有效实现位移重构。此外，改进方法被用于加筋圆柱壳的位移重构，数值结果表明，对于一些多构件组成的船舶结构，改进方法仍然有能力实现较为精确的位移重构。

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

Tessler A, Spangler J L. A least–squares variational method for full–field reconstruction of elastic deformations in shear–deformable plates and shells[J]. Computer Methods in Applied Mechanics & Engineering, 2005, 194(2–5): 327–339. DOI:10.1016/j.cma.2004.03.015.

Tessler A, Spangler J L. A variational principle for reconstruction of elastic deformations in shear deformable plates and shells: NASA/TM–2003–212445[R]. Hampton: NASA, 2003.

黄辉. 基于逆有限元方法的板梁结构变形重构的研究[D]. 湖北: 华中科技大学, 2021.

Huang H. Research on deformation reconstruction of plate and beam structure based on inverse finite element method[D]. Hubei: Huazhong University of Science and Technology, 2021. (in Chinese)

Kefal A, Oterkus E. Shape sensing of aerospace structures by coupling isogeometric analysis and inverse finite element method[C]//58th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. 2017. DOI:10.2514/6.2017–0427.

Kefal A, Oterkus E. Isogeometric iFEM analysis of thin shell structures[J]. Sensors, 2020, 20(9): 2685. DOI:10.3390/s20092685.

Abdollahzadeh M A, Belur M Y, Basoglu M F, et al, Shape sensing of beam–like structures using the robust iFEM–iQS4 inverse shell element[J]. IEEE Transactions on Instrumentation and Measurement, 2023, 72: 1–9.

Kefal A, Tessler A, Oterkus E. An efficient inverse finite element method for shape and stress sensing of laminated composite and sandwich plates and shells: NASA/TP–2018–220079[R]. Hampton: NASA, 2018.

Abdollahzadeh M A, Kefal A, Yildiz M. A comparative and review study on shape and stress sensing of flat/curved shell geometries using C0–continuous family of iFEM elements[J]. Sensors, 2020, 20(14): 3808.

Tessler A, Hughes T J R. A three–node mindlin plate element with improved transverse shear[J]. Computer Methods in Applied Mechanics & Engineering, 1985, 50(1): 71–101.

Quach C C, Vazquez S L, Tessler A, et al. Structural anomaly detection using fiber optic sensors and inverse finite element method[C]//Proceedings of AIAA Guidance, Navigation, and Control Conference and Exhibit, San Francisco, California, 2005.

Kefal A, Oterkus E, Tessler A, et al. A quadrilateral inverse–shell element with drilling degrees of freedom for shape sensing and structural health monitoring[J]. Engineering Science & Technology, an International Journal, 2016, 19(3): 1299–1313.

Kefal A. An efficient curved inverse–shell element for shape sensing and structural health monitoring of cylindrical marine structures[J]. Ocean Engineering, 2019, 188, 106262.

Li M, Kefal A, Cerik B C, et al. Structural health monitoring of submarine pressure hull using inverse finite element method[C]//7th International Conference on Marine Structures MARSTRUCT 2019, 2019.

Kefal A, Oterkus E. Shape and stress sensing of offshore structures by using inverse finite element method[C]//MARSTRUCT 2017, the 6th International Conference on Marine Structures. 2017.

Kefal, Adnan, Oterkus, et al. Displacement and stress monitoring of a chemical tanker based on inverse finite element method[J]. Ocean Engineering, 2016, 112(15): 33–46.

Kefal A, Oterkus E. Displacement and stress monitoring of a Panamax containership using inverse finite element method[J]. Ocean Engineering, 2016, 119: 16–29.

周富. 基于逆有限元法的重型数控机床立柱变形重构研究[D]. 湖北: 武汉理工大学, 2021.

Zhou F. Research on deformation reconstruction of heavy duty CNC machine tool column based on inverse finite element method[D]. Hubei: Wuhan University of Technology, 2021. (in Chinese)

赵雪晨. 基于逆有限元的板梁结构的损伤判别[D]. 陕西: 西安电子科技大学, 2022.

Zhao X C. Damage discrimination of plate and beam structures based on inverse finite element[D]. Shanxi: Xi'an University of Electronic Science and Technology, 2022. (in Chinese)

付书山, 孙广开, 何彦霖, 等. 基于逆有限元的机翼蒙皮变形监测方法仿真研究[J]. 航空制造技术, 2022, 65(6): 107‒114.

Fu S S, Sun G K, He Y L, et al. Simulation study of wing skin deformation monitoring method based on inverse finite element[J]. Aviation Manufacturing Technology, 2022, 65(6): 107‒114. (in Chinese)

徐浩, 王中枢, 马寅魏, 等. 基于逆有限元的蜂窝板结构形状传感研究[J]. 沈阳航空航天大学学报, 2022, 39(4): 44‒50.

Xu H, Wang Z S, Ma Y W, et al. A study on shape sensing of honeycomb panel structure based on inverse finite element[J]. Journal of Shenyang University of Aeronautics and Astronautics, 2022, 39(4): 44‒50. (in Chinese)

张科, 袁慎芳, 任元强, 等. 基于逆向有限元法的变形机翼鱼骨的变形重构[J]. 航空学报, 2020, 41(8): 244‒254.

Zhang K, Yuan S F, Ren Y Q, et al. Deformation reconstruction of deformed wing fishbone based on inverse finite element method[J]. Journal of Aeronautics, 2020, 41(8): 244‒254. (in Chinese)

闫宏生, 白超迪, 贾同宇, 等. 逆有限元法在船体加筋板结构变形重构中的应用[J]. 中国造船, 2023, 64(1): 168‒179.

Yan H S, Bai C D, Jia T Y, et al. Application of inverse finite element method in deformation reconstruction of hull stiffened plate structure[J]. China Shipbuilding, 2023, 64(1): 168‒179. (in Chinese)

Burke W L. Applied differential geometry[M]. Cambridge: Cambridge University Press, 1985.

刘云飞. 一种提高平板型壳单元计算精度的改进算法研究[D]. 大连: 大连理工大学, 2016.

Liu Y F. Research on an improved algorithm to increase the calculation accuracy of flat plate–shaped shell unit [D]. Dalian: Dalian University of Technology, 2016. (in Chinese)

Gao X W, Davies T G. Boundary Element Programming in Mechanics[M]. Cambridge: Cambridge University Press, 2002.

Lachat J C. A further development of the boundaryintegral technique for elastostatics[D]. Southampton: University of Southampton, 1975.

Becker A A. The boundary element method in engineering: A complete course[M]. New York: McGraw–Hill Companies, 1992.

2025年第29卷第12期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.3969/j.issn.1007-7294.2025.12.010

接收时间：2024-09-12
首发时间：2026-07-07
出版时间：2025-12-15

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-09-12

基金

作者信息

^1.天津大学水利工程仿真与安全国家重点实验室，天津 300072

^2.天津大学建筑工程学院，天津 300350

通讯作者:

闫宏生（1973–），男，博士，副教授，通讯作者，E-mail：hongsheng@tju.edu.cn

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/cblx/CN/10.3969/j.issn.1007-7294.2025.12.010

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT