CN114692446A - 一种矢量转动球坐标参数化的非线性壳有限元建模方法 - Google Patents

一种矢量转动球坐标参数化的非线性壳有限元建模方法 Download PDF

Info

Publication number
CN114692446A
CN114692446A CN202210233409.2A CN202210233409A CN114692446A CN 114692446 A CN114692446 A CN 114692446A CN 202210233409 A CN202210233409 A CN 202210233409A CN 114692446 A CN114692446 A CN 114692446A
Authority
CN
China
Prior art keywords
shell
vector
deformation
finite element
nonlinear
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210233409.2A
Other languages
English (en)
Inventor
王介龙
张发
徐吉峰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Commercial Aircraft Corp of China Ltd
Beijing Aeronautic Science and Technology Research Institute of COMAC
Original Assignee
Commercial Aircraft Corp of China Ltd
Beijing Aeronautic Science and Technology Research Institute of COMAC
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Commercial Aircraft Corp of China Ltd, Beijing Aeronautic Science and Technology Research Institute of COMAC filed Critical Commercial Aircraft Corp of China Ltd
Priority to CN202210233409.2A priority Critical patent/CN114692446A/zh
Publication of CN114692446A publication Critical patent/CN114692446A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种矢量转动球坐标参数化的非线性壳有限元建模方法，涉及非线性结构动力学、非线性有限元建模技术领域。该方法包括以下步骤：对薄壁三维结构进行二维中面壳抽取；将二维中面壳离散成规整四边形壳单元，所述四边形壳单元采用九节点二阶精度四边形壳单元；所述九节点二阶壳单元形函数中的五个参数在单元离散之后，在拉格朗日坐标系和矢量转动球坐标系下初始化各值；对所述壳单元建立结构运动控制方程。本申请提供了一种矢量转动球坐标参数化的非线性壳有限元建模方法。本申请技术方案将会定义中面法向矢量并直接推导壳的大变形与运动控制方程，建立新型几何精确壳单元。

Description

一种矢量转动球坐标参数化的非线性壳有限元建模方法
技术领域
本申请涉及非线性结构动力学、非线性有限元建模技术领域，具体涉及一种基于矢量转动球坐标系参数化有限元壳单元，可应用于飞机薄壁结构大变形非线性建模。
背景技术
为了减轻结构重量，民用飞机设计广泛采用了薄壁结构。当利用有限元技术对民机结构进行应力分析和强度校核时，壳单元非常适合于建立薄壁结构的有限元模型。因为处理复杂几何外形和边界条件时，壳单元具有很强的灵活性。过去几十年以来，壳理论研究受到了很大关注。根据采用的基础理论不同，现有的壳单元建模技术方案可分为两类，其一为实体单元降维建模，其二为直接建模。
技术方案一充分利用发展成熟的三维实体建模理论，通过将连续介质力学中的三维基本方程降阶，进而推导出二维壳单元控制方程，然后直接应用常规有限元插值技术获得相应的壳单元数值模型。该理论首先由Ahmad等人提出，之后获得了快速发展。代表性的工作包括Hughes等人以及Bathe和Dvorkin。商用程序NASTRAN中的线性壳单元QUAD4也属于这类技术方案。另外，Bathe，Hughes和Crisfield对壳的实体单元降维建模理论也进行了广泛深入的讨论。该技术方案背后的基础理论相对简单，非常适合于解决线性问题，但是计算效率低。该方案中，基于有限转动的参数化理论，采用两个自由度描述壳单元法向矢量的转动。例如，Hughes和Liu采用两个近似的Rodrigues参数描述转动。由于引入了小角度转动假设，该壳单元无法处理大旋转。
技术方案二与方案一不同，不采用降阶的三维方程而是定义垂直于壳单元中性面的法向矢量并直接推导控制方程描述壳等薄壁结构的变形与运动，形成了全新的壳理论。基于这一理论，可以构建更先进的有限元模型。这类技术方案首先由Argyris及其合作者提出，被称之为直接方法。Eriksen和Truesdell采用该方法，首次将壳简化为包含两个切向方向矢量的几何表面，发展了相应的壳理论。Simo等人开展了相似的工作，开发了基于应力-合力的几何精确壳模型。几何精确壳模型这一概念是指在推导壳的连续力学基本方程时，除最初引入的运动假设外，在推导过程中不再引入任何其他假设，并且在后续的推导过程中严格遵守引入的运动假设。因此，技术方案二中建立的壳单元精度非常高，适合于处理大变形和大旋转。Simo等人采用转动的Rodrigues表达式精确描述了有限转动。另外，Crisfiled和Moita[4]给出了共旋转表达式中对有限转动的处理。转动与线性位移不同，需要用旋转张量描述，因此转动的参数矢量化精确描述相当复杂，会耗费大量的计算机资源，降低计算效率，目前仍然是一个尚未彻底解决的问题，严重制约了技术方案二的应用。
现有技术方案一降阶连续介质力学中的三维基本方程，获得二维壳单元的控制方程，然后直接应用有限元插值技术获得相应的壳单元数值模型。建模理论相对简单，但是计算效率低，计算精度不高。通常采用两个转动自由度描述壳单元法向矢量的小角度转动，仅适合于解决线性小转动问题，无法处理大旋转。该技术方案无法解决薄壁结构在承载中的大变形问题。按照该技术方案设计的结构，有可能会过于保守，也可能造成灾难性的危险。
技术方案二直接推导控制方程描述壳等薄壁结构的变形与运动，是全新的壳理论。但是建模精度的提高使得理论推导相当复杂，计算效率极低。壳法向矢量转动的参数化描述进一步使问题复杂化。转动需要通过旋转张量精确描述，即使参数化后，仍然与位移矢量不同，不能直接采用有限元插值，需要进一步特殊处理。技术方案二的低效率限制了该方案的应用。
发明内容
针对以上问题，本申请提供了一种矢量转动球坐标参数化的非线性壳有限元建模方法。本申请技术方案将会定义中面法向矢量并直接推导壳的大变形与运动控制方程，建立新型几何精确壳单元。与常规处理方案不同，不对法向矢量的转动进行参数化处理，而是引入球坐标描述该矢量，极大简化了壳理论，提高了计算效率。法向矢量的球坐标描述未引入小角度转动假设，两个球坐标参数能够精确定义经过任意转动后的法向矢量。壳体结构的任意大旋转由变形构型中的法向矢量球坐标参数扣除初始值的相对量来描述，解决了技术方案一无法处理大旋转的问题。三个平动位移矢量分量与两个法向矢量的球坐标分量结合，精确描述了壳体单元的位移和转动。根据格林-拉格朗日应变张量定义，严格推导了全新的应变与平动位移和转动球坐标之间的全新表达式，保留了挠度应变二阶项。该新型应变-位移和球坐标表达式为严格几何非线性描述，能够处理大变形问题。另外，法向矢量的球坐标的引入，无须再对其进行进一步转动参数化处理，其本身就可以严格描述壳体结构的转动，节约了计算成本，提高了计算效率。因此，本技术方案解决了技术方案二转动参数复杂化、计算效率低的问题。
本发明是通过以下技术方案实现的：
根据本发明技术方案的第一方面，提供一种矢量转动球坐标参数化的非线性壳有限元建模方法，所述方法包括以下步骤：
(1)对薄壁三维结构进行二维中面壳抽取；
(2)将二维中面壳离散成规整四边形壳单元，所述四边形壳单元采用九节点二阶精度四边形壳单元；
(3)所述九节点二阶壳单元形函数中的五个参数在单元离散之后，在拉格朗日坐标系和矢量转动球坐标系下初始化各值；
(4)对所述壳单元建立结构运行控制方程。
进一步地，步骤(2)中，所述九节点二阶精度四边形壳单元形函数表示成
Figure BDA0003539412460000031
其中Hi为单元插值函数，ui为各节点坐标自由度ui＝[u1 u2 u3 θ φ]T，其中，u1、u2和u3是线位移三个分量，φ和θ为球坐标描述的转动位移两个分量。
进一步地，步骤(3)具体包括：
(31)所述九节点二阶壳单元插值函数Hi采用拉格朗日二次单元表达节点位移运动，将各节点线位移采用拉格朗日坐标系描述；
(32)所述九节点二阶壳单元插值函数Hi采用矢量转动球坐标系表达节点角位移大转动变形，引入球坐标定义法向矢量。
进一步地，步骤(31)中，在初始参考坐标系B中给定壳的厚度h和中面面积Ω，在参考构型中，给定壳中面上任意材料点的位置坐标x0和沿中面法向矢量
Figure BDA0003539412460000032
的材料坐标ζ，则壳中任意材料点的位置坐标为
Figure BDA0003539412460000033
其中，α1和α2是壳中面的材料点坐标，坐标α1，α2和ζ称之为曲线坐标，则变形构型中材料点的位置矢量为
Figure BDA0003539412460000034
式中，x0是壳中面材料点位置矢量，u(α1,α2)壳中面位移矢量，
Figure BDA0003539412460000035
为壳变形后的中面法向矢量。
进一步地，步骤(32)中法向矢量
Figure BDA0003539412460000041
定义为
Figure BDA0003539412460000042
如图3所示，法向矢量
Figure BDA0003539412460000043
在球坐标系中正交分解后，可以由球坐标变量θ和φ的方向正弦函数和方向余弦函数描述。
进一步地，步骤(4)中，所述壳单元的运动控制方程为
Figure BDA0003539412460000044
Figure BDA0003539412460000045
其中，符号
Figure BDA0003539412460000046
代表对时间的一阶导数，(·),1和(·),2代表针对坐标α1和α2的空间导数，(·)T表示矩阵转置；
Figure BDA0003539412460000047
和
Figure BDA0003539412460000048
代表惯性力贡献项，N 1、N 2和N 3为变形引起的弹性力，M 1和M 2为变形引起的弹性力矩；方程右边为外载荷贡献项，f为外力，m为外力矩，符号(·)*表示外力距m在局部坐标系中描述；Ef矩阵为球坐标转换矩阵,EB为Ef矩阵截断旋转矢量第三个分量对应列后的减缩矩阵。
壳理论中沿中面法向方向，即沿着厚度方向壳体没有变形，因此N 3不存在导数，只出现在力矩平衡方程中。另外，对应力矩M 3＝0，意味着忽略了平面内旋转自由度后仅保留5自由度，3个线位移2个角位移。上述方程中的下标1、2和3表示局部坐标系三个主轴方向，下标1和2在中性层内，下标3代表法向方向。
忽略上述方程中的惯性力，即获得静力学变形控制方程。
进一步地，步骤(4)中，给出如下应变-位移关系
Figure BDA0003539412460000049
Figure BDA00035394124600000410
Figure BDA00035394124600000411
Figure BDA00035394124600000412
Figure BDA00035394124600000413
上述应变-位移关系式适合于薄壁结构弹性体发生任意位移、大旋转和大应变的情况，其中，εij，i,j＝1,2,3为格林-拉格朗日应变张量非零分量，
Figure BDA00035394124600000414
和
Figure BDA00035394124600000415
为无量纲化的板壳中面正交切向矢量，R1和R2为薄壳单元的曲率半径。另外，ζ2对应的是挠度应变二阶项。
进一步地，步骤(4)中，所述壳单元的结构运行控制方程由Hamilton原理推导，具体推导过程如下：
在时间域内对任意材料点的位置矢量进行求导，得速度矢量
Figure BDA0003539412460000051
其中，矢量X为薄壳中沿中面法向任意材料的变形后在惯性系中的位置矢量。矢量u为薄壳中面参考点的位移。由此，壳的动能K可写为
Figure BDA0003539412460000052
其中，ρ为材料密度。动能变分可写为
Figure BDA0003539412460000053
式中h和
Figure BDA0003539412460000054
分别代表线性动量矩矢量和角动量矩矢量，具体表达式可以从动能变分过程中获得。
壳的应变能变分是
Figure BDA0003539412460000055
(9)
其中τ*ij(i,j＝1,2,3)为随体坐标系中描述的Cauchy应力张量；将应变分量带入应变能变分表达式并沿厚度积分可得应变能变分的等价形式
Figure BDA0003539412460000056
其中，N 1、N 2和N 3为在薄壳中面坐标系中描述的弹性力，M 1和M 2为在薄壳中面坐标系中描述的弹性力矩。外力虚功表示为
δW＝∫Ω(δu T f+δψ T m)dΩ (11)
其中，δψ代表虚角位移，f和m代表壳中面单位面积上作用的力和力矩，将动能变分、应变能变分和外力虚功带入Hamilton原理，最终得到壳的运动控制方程
Figure BDA0003539412460000057
Figure BDA0003539412460000058
忽略上述方程中的惯性力，即可获得静力学变形控制方程。
进一步地，步骤(4)中的壳体结构运动控制方程分为两种情况，一是壳体动力学常微分方程，而是壳体结构静力学非线性几何方程。均通过线化方程迭代法求解，每次迭代预测变形增量。变形未满足收敛要求时，需要更新刚度、变形和弹性载荷，重复以上步骤总迭代次数最多10次；变形迭代结果满足收敛要求时，即线性迭代收敛，结构任意大变形响应为随求结果，结束以上求解过程。
根据本发明技术方案的第二方面，提供一种根据以上任一方面所述的非线性壳有限元建模方法应用于建模构建基于矢量转动球坐标系参数化有限元壳单元的用途。
根据本发明技术方案的第三方面，提供一种基于矢量转动球坐标系参数化有限元壳单元，其特征在于，所述基于矢量转动球坐标系参数化有限元壳单元采用根据以上任一方面所述的非线性壳有限元建模方法进行建模构建。
本发明有益效果是：
1.该单元除引入初始小挠度假设外，未引入其他任何假设，因此可以直接处理大变形问题。
2.未对法向矢量的转动进行参数化处理，而是在球坐标系中定义该矢量，采用两个球坐标参数描述了法向矢量的任意转动。
3.相对于转动的张量描述，转动的球坐标参数化极为简化、高效。
4.本申请提案的技术方案兼顾效率与精度。
附图说明
图1示出本申请提案的技术方案流程图；
图2示出9节点2阶四边形壳单元；
图3示出法向矢量的球坐标描述；
图4示出扭曲平板构型；
图5示出扭曲平板的变形；
图6示出扭曲平板自由端变形的时间历程；
图7示出半圆柱形薄壁壳体以及有限元网格；
图8示出半圆柱形薄壁壳体加载点位移的时间历程；
图9示出半圆柱形薄壁壳体形心曲率应变的时间历程。
以下结合附图对本发明的实施方式作进一步说明。
具体实施方式
另外，下面结合附图描述的本发明的实施方式是示例性的，仅用于解释本发明的实施方式，而并非对本发明的限制。
本申请技术方案将详细介绍新型壳单元的理论细节及利用该技术方案完成的薄壁结构大变形和非线性振动分析。其中，本申请建模方法提供一种保证精度前提下的简单高效求解几何大变形有限元壳单元，在球坐标系下利用三个位移分量和两个球坐标参数描述壳单元大变形，解决技术方案一中传统壳单元降阶后无法处理大旋转的问题，解决技术方案二中大规模参数化描述单元转动变形导致的计算效率低的问题。此单元是一种新型高精度、高效率的壳单元，其建模技术流程如图1所示。
本发明是通过以下技术方案实现的：一种矢量转动球坐标参数化的非线性壳有限元建模方法，方法包括以下步骤：
(1)对薄壁三维结构进行二维中面壳抽取；
(2)将中面壳离散成规整四边形壳单元，单元数量可基于薄壁结构面内尺寸而定，所述四边形壳单元为提高分析精度，采用采用常规的九节点二阶精度四边形壳单元，如图2所示。单元形函数可以表示成
Figure BDA0003539412460000071
其中Hi为单元插值函数，ui为各节点坐标自由度
Figure BDA0003539412460000072
(3)九节点二阶壳单元形函数中的五个参数在单元离散之后，在相应的坐标系下初始化各值：
(31)九节点二阶壳单元插值函数Hi采用拉格朗日二次单元表达节点位移运动：各节点线位移采用拉格朗日坐标系描述。
在初始参考坐标系B中给定壳的厚度h和中面面积Ω。在参考构型中，给定壳中面上任意材料点的位置坐标x0和沿中面法向矢量
Figure BDA0003539412460000073
的材料坐标ζ。壳中任意材料点的位置坐标可写为
Figure BDA0003539412460000074
其中，α1和α2是壳中面的材料点坐标。坐标α1，α2和ζ称之为曲线坐标。相应的，变形构型中材料点的位置矢量可写为
Figure BDA0003539412460000075
式中，x0是壳中面材料点位置矢量，u(α1,α2)壳中面位移矢量。
Figure BDA0003539412460000076
为壳变形后的中面法向矢量。方程(2)表明任意材料点的变形与位移可分解为中面材料点位移u(α1,α2)和对应法向矢量的转动
Figure BDA0003539412460000081
(32)九节点二阶壳单元插值函数Hi采用矢量转动球坐标表达节点角位移大转动变形重新构建。
本申请提案的技术方案引入球坐标定义法向矢量，并将其作为未知量直接求解。避免了转动张量的参数化，极大简化了问题，同时能够直接应用有限元插值。法向矢量
Figure BDA0003539412460000082
可定义为
Figure BDA0003539412460000083
相应的几何描述参见图3。
球坐标φ和θ作为未知量通过壳的控制方程来求解。在参考构型中，这两个坐标的初始值为φ0和θ0，定义了初始法向矢量
Figure BDA0003539412460000084
现有技术方案均采用转动张量的参数化表达式来描述法向矢量的转动。本发明解决了现有技术方案一采用线性小转动假设，引入两个转动参数，不能处理大旋转；和现有技术方案二采用Rodrigues表达式等转动参数化形式描述有限转动，使转动问题进一步复杂化，不能直接应用有限元插值等问题。
(4)完成形函数各值坐标系定义，对所述壳单元建立结构运行控制方程方程。所述壳单元是几何非线性结构有限元单元；将基于上述转动的球坐标描述结构大转动变形，具有较强非线性特征，重新推导壳单元基本方程，最终开发出全新的非线性几何精确壳单元。
壳单元在应变定义中，引入Mindlin/Reissner假设：与壳中面垂直的材料线变形后仍然保持直线，不发生拉伸变形。在该假设中未约束壳变形后材料线仍然垂直中面，因此会发生剪切变形。该假设适用于中厚度壳。依据弹性力学理论中的格林-拉格朗日应变张量的定义，本申请提案的技术方案给出了如下应变-位移关系
Figure BDA0003539412460000085
Figure BDA0003539412460000086
Figure BDA0003539412460000091
Figure BDA0003539412460000092
Figure BDA0003539412460000093
上述应变-位移关系式由于引入球坐标矢量，因此是一种全新表达式，适合于任意位移、大旋转和大应变。其中，ζ2对应的是挠度应变二阶项。现有技术方案二认为该项为高阶小量，将其忽略。所得应变-位移关系式仅仅适合于小应变情况。本申请提案的技术方案利用从初始时间ti到当前时间tf数值积分过程中，应变的增量变化很小这一事实，将ζ2对应的挠度应变二阶项由其在时间ti时的初始值近似，并在每次积分逐步更新。因此，新型壳单元能够解决大变形问题。
壳单元的运动控制方程由Hamilton原理推导。在时间域内对任意材料点的位置矢量，方程(2)，进行求导，可得速度矢量
Figure BDA0003539412460000094
由此，壳的动能为
Figure BDA0003539412460000095
动能变分可写为
Figure BDA0003539412460000096
式中h和
Figure BDA0003539412460000097
分别代表线性动量矩矢量和角动量矩矢量。
壳的应变能变分是
δV＝∫Ω∫h(τ*11δε11+τ*22δε22+τ*12δε12+τ*13δε13+τ*23δε23)dζdΩ
(9)
其中τ*ij(i，j＝1，2，3)为随体坐标系中描述的Cauchy应力张量；将应变分量带入应变能变分表达式并沿厚度积分可得应变能变分的等价形式
Figure BDA0003539412460000098
外力虚功可表示为
δW＝∫Ω(δu T f+δψ T m)dΩ (11)
其中f和m代表壳中面单位面积上作用的力和力矩。将动能变分、应变能变分和外力虚功带入Hamilton原理，最终得到壳的运动控制方程
Figure BDA0003539412460000101
Figure BDA0003539412460000102
忽略上述方程中的惯性力，即可获得静力学变形控制方程。
基于以上结构单元控制节点变形方程，遵循经典有限元分析流程，包括赋予有限元单元质量属性、单元法向向量、材料向量和单元刚度矩阵、载荷向量及边界约束、和阻尼特性；基于Newton-Raphson迭代算法的变形计算、计算结果的后处理获得结构变形、应变和应力分布。
运动控制方程采用线化方程迭代法，当变形不满足收敛性要求时，计算出当前步骤的变形增量后，更新刚度、变形和弹性力矢量，进一步迭代，重复以上步骤总迭代次数最多10次；如变形迭代满足收敛条件，即迭代收敛，结构任意大变形响应为随求结果。结束以上求解过程。
步骤(1)～(3)，依据上述推导出的公式并结合有限元理论，本申请提案的技术方案完成了九节点二阶精度四边形壳单元的开发，该单元不仅具有静力学分析功能，而且还具备模态分解和动力学分析功能。
实施例
根据本实施例的一种基于矢量转动球坐标参数化的非线性有限元建模方法具体实施方式包含三个位移分量和两个球坐标参数的新型5自由度壳单元可以直接进行有限元插值，是一种新型高精度、高效率的壳单元。该单元建模的技术流程图如图1所示。
扭曲平板的大变形弯曲分析
矩形截面的悬臂平板沿展向扭转90°，并在自由端承受集中载荷，如图4所示。扭曲平板的长度为L＝0.3048米，宽度b＝0.0279米，厚度0.00813米，离散成2×8单元网格(宽度方向2个壳单元，展向8个壳单元)。该矩形板的材料属性设置为：杨氏模量E＝200×109帕斯卡，泊松比ν＝0.22。
为了验证所开发壳单元的精度以及处理大变形能力，对该单元进行了动力学分析。预测结果和几何精确梁单元[5]模型的计算结果进行了对比验证。在固支边界条件下，采用当前开发的几何精确壳单元和公开发表的几何精确梁单元两种不同的单元预测了扭曲矩形板的变形。在自由端，P＝2.224×104牛顿集中剪切载荷作用下扭曲平板发生变形。变形后的构型如图5所示。该剪切载荷在5秒内从零线性增加到2.224×104牛顿，然后保持该值不变。
结果表明，所开发的几何精确壳单元具有预测大变形的能力。与几何精确梁单元模型计算结果的对比验证可参见图6。图中描述了自由端中点处变形的时间历程，实线为梁模型结果，虚线为壳模型结果，二者计算结果基本保持一致。
半圆柱形薄壁壳体应力计算
半圆柱形薄壁壳体一端固支，一端自由，半径为1.2米，长度为2米，其受载情况如图7所示。集中载荷大小为F1＝F2＝F3＝1.0×102牛顿。该薄壁结构离散成4×8的有限元网格(展向4个壳单元，弧向8个壳单元)。材料属性：杨氏模量E＝73×109帕斯卡，泊松比ν＝0.3，壁板厚度为0.006米。
采用新型壳单元预测了该薄壁结构在集中载荷作用下的动力响应。
加载点位移和转动的时间历程可参考图8，薄壁壳体几何形心的曲率应变可参考图9。图8和9中，本申请提案的技术方案计算结果与有限元程序DYMORE[6]的计算结果进行了比较。图8和9中实线为DYMORE计算结果，虚线为新型几何非线性壳单元计算结果。二者的计算结果一致，表明新型几何非线性壳单元具有很好的精度。
上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (10)

1.一种矢量转动球坐标参数化的非线性壳有限元建模方法，其特征在于，所述方法包括以下步骤：
(1)对薄壁三维结构进行二维中面壳抽取；
(2)将二维中面壳离散成规整四边形壳单元，所述四边形壳单元采用九节点二阶精度四边形壳单元；
(3)所述九节点二阶壳单元形函数中的五个参数在单元离散之后，在拉格朗日坐标系和矢量转动球坐标系下初始化各值；
(4)对所述壳单元建立壳体结构运动控制方程。
2.根据权利要求1所述的非线性壳有限元建模方法，其特征在于，所述九节点二阶精度四边形壳单元形函数表示成
Figure FDA0003539412450000011
其中Hi为单元插值函数，ui为各节点坐标自由度ui＝[u1 u2 u3 θ φ]T，其中，u1、u2和u3是线位移三个分量，φ和θ为球坐标描述的转动位移两个分量。
3.根据权利要求1所述的非线性壳有限元建模方法，其特征在于，步骤(3)具体包括：
(31)所述九节点二阶壳单元插值函数Hi采用拉格朗日二次单元表达节点位移运动，将各节点线位移采用拉格朗日坐标系描述；
(32)所述九节点二阶壳单元插值函数Hi采用矢量转动球坐标系表达节点角位移大转动变形，引入球坐标定义法向矢量。
4.根据权利要求3所述的非线性壳有限元建模方法，其特征在于，步骤(31)中，在初始参考坐标系B中给定壳的厚度h和中面面积Ω，在参考构型中，给定壳中面上任意材料点的位置坐标x 0和沿中面法向矢量
Figure FDA0003539412450000012
的材料坐标ζ，则壳中任意材料点的位置坐标为
Figure FDA0003539412450000013
其中，α1和α2是壳中面的材料点坐标，坐标α1，α2和ζ称之为曲线坐标，则变形构型中材料点的位置矢量为
Figure FDA0003539412450000014
式中，x 0是壳中面材料点位置矢量，u(α1，α2)壳中面位移矢量，
Figure FDA0003539412450000015
为壳变形后的中面法向矢量。
5.根据权利要求3所述的非线性壳有限元建模方法，其特征在于，步骤(32)中法向矢量
Figure FDA0003539412450000021
定义为
Figure FDA0003539412450000022
6.根据权利要求1所述的非线性壳有限元建模方法，其特征在于，步骤(4)中，所述壳单元的运动控制方程为
Figure FDA0003539412450000023
Figure FDA0003539412450000024
其中，符号
Figure FDA0003539412450000025
代表对时间的一阶导数，(·)，1和(·)，2代表针对坐标α1和α2的空间导数，(·)T表示矩阵转置；
Figure FDA0003539412450000026
和
Figure FDA0003539412450000027
代表惯性力贡献项，N 1、N 2和N 3为变形引起的弹性力，M 1和M 2为变形引起的弹性力矩；方程右边为外载荷贡献项，f为外力，m为外力矩，符号(·)*表示外力距m在局部坐标系中描述；Ef矩阵为球坐标转换矩阵，EB为Ef矩阵截断旋转矢量第三个分量对应列后的减缩矩阵。
7.根据权利要求1所述的非线性壳有限元建模方法，其特征在于，步骤(4)中，包括应变-位移关系：
Figure FDA0003539412450000028
Figure FDA0003539412450000029
Figure FDA00035394124500000210
Figure FDA00035394124500000211
Figure FDA00035394124500000212
上述应变-位移关系式适合于薄壁结构弹性体发生任意位移、大旋转和大应变的情况，其中，εij，i，j＝1，2，3为格林-拉格朗日应变张量非零分量，
Figure FDA00035394124500000213
和
Figure FDA00035394124500000214
为无量纲化的板壳中面正交切向矢量，R1和R2为薄壳单元的曲率半径，ζ2对应的是挠度应变二阶项。
8.根据权利要求1所述的非线性壳有限元建模方法，其特征在于，步骤(4)中，所述壳单元的结构运行控制方程由Hamilton原理推导，具体推导过程如下：
在时间域内对任意材料点的位置矢量进行求导，得速度矢量
Figure FDA0003539412450000031
其中，矢量X为薄壳中沿中面法向任意材料点变形后在惯性系中定义的位置矢量，矢量u为薄壳中面参考点的位移，壳的动能K可写为
Figure FDA0003539412450000032
其中，ρ为材料密度，动能变分可写为
Figure FDA0003539412450000033
式中h和
Figure FDA0003539412450000034
分别代表线性动量矩矢量和角动量矩矢量，具体表达式可以从动能变分过程中获得，壳的应变能变分是
δV＝∫Ω∫h(τ*11δε11+τ*22δε22+τ*12δε12+τ*13δε13+τ*23δε23)dζdΩ
其中τ*ij(i，j＝1，2，3)为随体坐标系中描述的Cauchy应力张量；将应变分量带入应变能
变分表达式并沿厚度积分可得应变能变分的等价形式
Figure FDA0003539412450000035
其中，N 1、N 2和N 3为在薄壳中面坐标系中描述的弹性力，M 1和M 2为在薄壳中面坐标系中描述的弹性力矩，外力虚功表示为
δW＝∫Ω(δu T f+δψ T m)dΩ
其中f和m代表壳中面单位面积上作用的力和力矩，将动能变分、应变能变分和外力虚功带入Hamilton原理，最终得到壳的运动控制方程
Figure FDA0003539412450000036
Figure FDA0003539412450000037
忽略上述方程中的惯性力，即可获得静力学变形控制方程。
9.根据权利要求1所述的非线性壳有限元建模方法，其特征在于，步骤(4)中，壳体结构动力学平衡方程为常微分方程，壳体结构静力学为非线性几何方程，均通过线化方程迭代法求解，根据迭代结果预测变形增量，变形未满足收敛要求时，更新刚度、变形和弹性载荷，重复以上步骤总迭代次数最多10次；变形迭代结果满足收敛要求时，即线性迭代收敛，结构任意大变形响应为随求结果，结束以上求解过程。
10.一种基于矢量转动球坐标系参数化有限元壳单元，其特征在于，所述基于矢量转动球坐标系参数化有限元壳单元采用根据权利要求1至9中任一项所述的非线性壳有限元建模方法进行建模构建。
CN202210233409.2A 2022-03-09 2022-03-09 一种矢量转动球坐标参数化的非线性壳有限元建模方法 Pending CN114692446A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210233409.2A CN114692446A (zh) 2022-03-09 2022-03-09 一种矢量转动球坐标参数化的非线性壳有限元建模方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210233409.2A CN114692446A (zh) 2022-03-09 2022-03-09 一种矢量转动球坐标参数化的非线性壳有限元建模方法

Publications (1)

Publication Number Publication Date
CN114692446A true CN114692446A (zh) 2022-07-01

Family

ID=82138400

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210233409.2A Pending CN114692446A (zh) 2022-03-09 2022-03-09 一种矢量转动球坐标参数化的非线性壳有限元建模方法

Country Status (1)

Country Link
CN (1) CN114692446A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115879279A (zh) * 2022-11-18 2023-03-31 上海交通大学 换流变有载分接开关应力薄弱点分析方法及系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115879279A (zh) * 2022-11-18 2023-03-31 上海交通大学 换流变有载分接开关应力薄弱点分析方法及系统
CN115879279B (zh) * 2022-11-18 2024-03-08 上海交通大学 换流变有载分接开关应力薄弱点分析方法及系统

Similar Documents

Publication Publication Date Title
Lang et al. Multi-body dynamics simulation of geometrically exact Cosserat rods
Yu et al. Generalized Timoshenko theory of the variational asymptotic beam sectional analysis
Mikkola et al. A non-incremental finite element procedure for the analysis of large deformation of plates and shells in mechanical system applications
Mian et al. Numerical investigation of structural geometric nonlinearity effect in high-aspect-ratio wing using CFD/CSD coupled approach
CN104281730B (zh) 一种大转动变形的板壳结构动响应的有限元分析方法
Park et al. Prediction of damping coefficients using the unsteady Euler Equations
Almeida et al. Corotational nonlinear dynamic analysis of laminated composite shells
Sun et al. Topology optimization based on level set for a flexible multibody system modeled via ANCF
Bekhoucha et al. Nonlinear forced vibrations of rotating anisotropic beams
CN112580241B (zh) 一种基于结构降阶模型的非线性气动弹性动响应分析方法
Ren et al. An accurate and robust geometrically exact curved beam formulation for multibody dynamic analysis
Pan et al. Geometric nonlinear dynamic analysis of curved beams using curved beam element
Fan et al. A novel numerical manifold method with derivative degrees of freedom and without linear dependence
Funes et al. An efficient dynamic formulation for solving rigid and flexible multibody systems based on semirecursive method and implicit integration
CN114692446A (zh) 一种矢量转动球坐标参数化的非线性壳有限元建模方法
Fang et al. An efficient radial basis functions mesh deformation with greedy algorithm based on recurrence Choleskey decomposition and parallel computing
Patil et al. A scalable time-parallel solution of periodic rotor dynamics in X3D
Zhang et al. An asynchronous parallel explicit solver based on scaled boundary finite element method using octree meshes
CN108182330A (zh) 一种基于b样条计算柔性矩形薄板刚柔耦合动力学响应的方法
Hodges et al. Geometrically-exact, intrinsic theory for dynamics of moving composite plates
Mavriplis et al. Recent advances in high-fidelity multidisciplinary adjoint-based optimization with the NSU3D flow solver framework
Ren et al. An adaptive triangular element of absolute nodal coordinate formulation for thin plates and membranes
Yuan et al. A general nonlinear order-reduction method based on the referenced nodal coordinate formulation for a flexible multibody system
CN104091003A (zh) 一种基础运动时柔性壳结构大变形响应的有限元建模方法
Le et al. A thin-walled composite beam model for light-weighted structures interacting with fluids

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination