基于高阶复合曲面单元的非均匀大薄板装配偏差计算模型

doi:10.16183/j.cnki.jsjtu.2020.333

基于高阶复合曲面单元的非均匀大薄板装配偏差计算模型

林张鹏, 余海东^,, 袁可

上海交通大学上海市复杂薄板结构数字化制造重点实验室,上海 200240

Model of Assembly Deviation of Irregular Large Thin-Walled Structures Based on Higher-Order Composite Shell Element

LIN Zhangpeng, YU Haidong^,, YUAN Ke

Shanghai Key Laboratory of Digital Manufacture for Thin-Walled Structures, Shanghai Jiao Tong University, Shanghai 200240, China

通讯作者: 余海东,男,教授;电话(Tel.):021-34206542,E-mail:hdyu@sjtu.edu.cn.

责任编辑: 孙伟

收稿日期: 2020-10-20

基金资助:

国家重点研发计划(2019YF0709000)
国家自然科学基金(51775345)

Received: 2020-10-20

作者简介 About authors

林张鹏(1996-),男,澳门特别行政区人,硕士生,主要研究方向为薄壁结构装配偏差预测与控制.

摘要

结构的刚度非均匀分布易导致非规则大型薄壁在加工和装配过程中产生翘曲大变形,其结构建模和变形精确计算是装配偏差预测及控制的关键.基于绝对节点坐标法,引入高阶斜率坐标并改进单元插值函数,考虑单元接触界面间连续性,构建一种非规则高阶曲面复合单元模型;采用该非规则高阶曲面单元建立复合单元组合离散非规则曲面结构,解决不同厚度单元间节点连接和变形耦合问题,并有效地描述非规则大型薄壁结构的大变形,基于连续介质力学建立非规则大型薄壁结构装配协调变形计算模型.对重型火箭贮箱箱底边缘区域加厚的超大型非规则薄壁结构进行装配过程变形计算,研究加厚区域尺寸与瓜瓣装配偏差以及结构刚度之间的关联关系,为未来重型火箭箱底超大非规则曲面结构设计提供指导.

关键词： 非规则薄壁结构; 绝对节点坐标法; 复合单元; 装配偏差; 协调变形

Abstract

The non-uniform stiffness distribution of irregular large thin-walled structures leads to warping deformation in the process of machining and assembly. The structure modeling and accurate calculation of the large deformation is the key to predict and control the deviation. A novel irregular higher-order composite shell element based on the absolute nodal coordinate formulation (ANCF) was proposed to accurately calculate the deformation of irregular large thin-walled structures with higher-order slope coordinates and improved shape function, considering the continuity of elements at the contact interface. A discrete method of irregular large thin-walled structures was proposed to solve the problem of node connection and deformation coupling between elements with different thicknesses, and to describe the large deformation of irregular large thin-walled structures effectively by using the combination of irregular higher-order elements and composite elements. An assembly deformation calculation model of irregular large thin-walled structures was proposed based on continuum mechanics. The influence of the thickened area on assembly deviation and stiffness of scalloped segment structures of heavy rocket tank bottom was studied with the new elements and calculation model, so as to provide guidance for structure design.

Keywords： irregular thin-walled structure; absolute nodal coordinate formulation (ANCF); composite element; assembly deviation; compatible deformation

PDF (1933KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

林张鹏, 余海东, 袁可. 基于高阶复合曲面单元的非均匀大薄板装配偏差计算模型[J]. 上海交通大学学报, 2022, 56(5): 584-593 doi:10.16183/j.cnki.jsjtu.2020.333

LIN Zhangpeng, YU Haidong, YUAN Ke. Model of Assembly Deviation of Irregular Large Thin-Walled Structures Based on Higher-Order Composite Shell Element[J]. Journal of Shanghai Jiaotong University, 2022, 56(5): 584-593 doi:10.16183/j.cnki.jsjtu.2020.333

厚度非均匀分布的非规则大型薄壁结构件具有重量轻、强度高等优点,广泛运用于航空航天、船舶等领域^[1].非规则大型薄壁结构件的几何特征造成其结构刚度的非均匀分布,在加工和装配过程中易产生变形以及转动耦合的非线性翘曲变形,导致结构存在较大偏差.结构的复杂性使得利用解析法求解其非线性变形十分困难.而在数值仿真中,为满足精度要求,传统有限元法受单元线性假设限制,离散结构需要大量单元而导致计算效率低.因此,针对非规则大型薄壁结构,构建合理的高阶单元模型以及采用有效的结构离散方法是描述其翘曲变形并提高计算效率的关键.

非均匀大型薄壁结构的变形和应力的高精度分析一直是其制造偏差精确计算的难点.现有有限元法主要采用三维实体单元模型或基于板壳理论的二维壳单元模型对结构进行离散.三维实体单元能够精确计算结构变形及力学特性,但对存在凸台和承受弯曲变形区域的计算效率低;二维壳单元的单元自由度和数量较少,计算效率较高,但忽略了截面变形且用于厚度非均匀结构的建模存在局限.

利用等效单元法^[2]和复合单元法^[3]可以解决壳单元离散的局限性并减少单元数量.等效单元法通过内力等效将附加厚度结构的作用等效至主体结构,但无法描述两者之间的变形耦合^[4];复合单元法将主体及附加厚度结构等效为单元叠加,通过单元间变形协调及连续性条件构建复合单元模型.针对板类结构,Carrera^[5]总结了两种常用的板壳叠层结构建模方法, 其中假定各层单元均以单层板形式变形的单一板模型忽略了柔性结构的非线性大变形特征,Layer-wise模型(一种预训练语言模型)通过在各单层板接触界面处插入高阶边界条件提高计算精度,但结构自由度的增加导致计算效率低,且该模型常用于复合材料结构建模及变形计算^[6].目前的复合单元主要基于小变形和小转动假设,在针对大变形和大旋转的非线性问题方面,Shabana等^[7-8]提出绝对节点坐标法(ANCF),用斜率坐标替代传统的转角来精确描述单元截面变形,并提出非线性弹性力和大变形的计算方法;Matikainen等^[9]通过在厚度方向引入二阶位移插值构建高阶ANCF板单元,克服了弯曲产生的泊松锁定;Ebel等^[10]通过增加单元节点消除前者产生的附加剪切锁定现象,提高了计算精度.Abbas等^[11]通过修正板单元的形函数构建变厚度的板单元模型,消除了ANCF对规则单元模型的限制,并利用模态分析以及ANSYS软件的对比结果验证了单元模型的正确性.基于复合单元法,Yu等^[12]考虑板单元间接触界面的变形协调条件,构建全参数的板/板单元模型来离散含筋板结构,消除了筋板穿透现象且具有较高的计算精度,但该单元在描述弯曲单元时会产生泊松锁闭问题.现有的ANCF单元主要用于离散厚度非均匀的规则薄壁结构或单一厚度的非规则薄壁结构,而针对厚度非均匀且形状非规则薄壁结构离散的高阶单元模型研究较少.

基于ANCF,构建非规则高阶曲面复合单元模型并计算单元非线性刚度矩阵,精确描述结构的非线性翘曲变形;利用非规则高阶曲面单元及复合单元组合离散非规则大型薄壁结构,解决不同厚度单元间节点连接和变形耦合问题;基于连续介质力学建立非规则薄壁结构装配协调变形计算模型.对重型火箭箱底瓜瓣结构进行装配偏差变形计算仿真,并研究边缘加厚区域尺寸对瓜瓣装配偏差和结构刚度的影响.

1 非规则高阶曲面复合单元建模

1.1 非规则大型薄壁结构离散方法

工程中复杂的薄壁曲面结构常由对边不等且各区域厚度不同的非规则大型薄壁加工装配而成.针对非规则的形状特征,为保证模型准确性并减少单元数量,引入高阶插值函数并改进单元形函数,构建边长可独立变化的高阶曲面单元.而采用同一类型的曲面单元离散存在厚度变化的曲面结构时,不同厚度区域单元间的共用节点不在同一中性面上,导致单元间无法使用共同节点进行连接.因此,需要根据厚度分布将非规则大型薄壁结构划分为多个区域进行离散.最小厚度区域使用非规则高阶曲面单元进行离散,而其余区域均使用非规则高阶曲面复合单元进行离散.不同厚度单元间的耦合建模如图1所示,其中(上)下标P表示下层单元.上层单元与最小区域单元的厚度相同,并通过公共节点A(B_I)和D(C_I) 连接, 而下层单元节点通过转换矩阵由上层单元表示因此曲面结构变形可以通过最小区域单元节点(A_I~D_I)和复合单元的上层单元节点(A~D)表示.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 非规则高阶曲面单元及其复合单元耦合建模

Fig.1 Coupling of irregular higher-order shell element and composite element

以贮箱箱底瓜瓣结构为例,其离散过程如图2所示.其中,N为节点数,j₁为非规则单元,其4个节点编号分别为k、k+1、l和l+1,j₂为复合单元.其中间区域与边缘区域的厚度不同,因此分别采用非规则高阶曲面单元与复合单元进行离散.非规则高阶曲面单元和复合单元组合离散非规则大型薄壁结构能够有效解决单层单元无法离散厚度非均匀结构的问题.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 非规则大型薄壁结构区域离散方法

Fig.2 Meshing of irregular large thin-walled structure

1.2 非规则高阶曲面单元模型

基于ANCF,构建非规则高阶曲面单元模型,如图3所示.其中,a_i、b_i和t_i(i=1, 2, 3, 4)是可独立变化的单元边长,r和r₀分别为变形前后单元上任意一点的广义位置坐标,点P为单元上任意一点.单元边长可独立变化且由于引入高阶插值函数,各边均由高阶曲线描述,并能用较少数量单元描述复杂的非规则曲面结构.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 具有可变边长的非规则高阶ANCF单元模型

Fig.3 Higher-order ANCF element with variable side lengths

非规则高阶曲面单元的4个节点A、B、C和D由位置坐标和斜率坐标两部分组成.沿单元厚度z方向增加代表二次斜率坐标的高阶项能够减弱泊松锁闭问题并提高数值仿真准确性,因此单元每个节点坐标包含3个位置坐标和12个斜率坐标,单元共有60个节点坐标,定义为

e^PE60=

$[\begin{array}{l} {e^{T}}_{A} & {e^{T}}_{B} & {e^{T}}_{C} & {e^{T}}_{D} \end{array}]$

(1)

式中:

e_i=

$[\begin{array}{l} {r^{T}}_{i} & \frac{\partial {r^{T}}_{i}}{\partial x} & \frac{\partial {r^{T}}_{i}}{\partial y} & \frac{\partial {r^{T}}_{i}}{\partial z} & \frac{\partial^{2} {r^{T}}_{i}}{\partial z \partial z} \end{array}]$
i=A, B, C, D

(2)

故单元上任意一点坐标可以表示为形函数矩阵与节点广义坐标的乘积:

r=Se

(3)

式中:

S=[S₁I S₂I S₃I S₄I ··· S ₁₈I S₁₉I S₂₀I]

(4)

式中:I为60×60的单位矩阵,且

S₁=(2ξ+1)(ξ-1)²(2η+1)(η-1)²

S₂=a₁ξ(ξ-1)²(2η+1)(η-1)²

S₃=b₁η(ξ-1)²(2ξ+1)(η-1)²

S₄=t₁ζ(ξ-1)(η-1)

S₅=

$\frac{t_{1}^{2} ζ^{2} (η - 1) (ξ - 1)}{2}$

S₆=-ξ²(2ξ-3)(2η+1)(η-1)²

S₇=a₁ξ²(ξ-1)(2η+1)(η-1)²

S₈=-b₂ηξ²(2ξ-3)(η-1)²

S₉=-t₂ζξ(η-1), S₁₀=

$\frac{t_{2}^{2} ζ^{2} ξ (1 - η)}{2}$

S₁₁=η²ξ²(2ξ-3)(2η-3)

S₁₂=-a₂η²ξ²(ξ-1)(2η-3)

S₁₃=-b₂η²ξ²(η-1)(2ξ-3)

S₁₄=t₃ζξη, S₁₅=

$\frac{t_{3}^{2} ζ^{2} η ξ}{2}$

S₁₆=-η²(2ξ+1)(ξ-1)²(2η-3)

S₁₇=-a₂ξη²(ξ-1)²(2η-3)

S₁₈=b₁η²(ξ-1)²(2ξ+1)(η-1)

S₁₉=-t₄ηζ(ξ-1), S₂₀=

$\frac{t_{4}^{2} ζ^{2} η (1 - ξ)}{2}$

(5)

式中:ξ=x/a₁、η=y/b₁和ζ=z/t₁为无量纲变量.

1.3 非规则高阶曲面复合单元模型

基于ANCF,构建非规则高阶曲面复合单元模型,如图4所示.其中,(上)下标M表示上层单元,l₁和l₂为上层单元中性面边长的一半.在上层单元与下层单元间引入高阶斜率约束条件,使得单元间的几何连续性达到C²连续,消除单元间变形耦合出现的穿透问题,保证复合单元的合理性并提高数值计算精度.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 上层单元中性面及复合单元接触横截面

Fig.4 Neutral surface of main element and cross-sections of composite element

定义下层单元与上层单元接触点分别为E_P和E、F_P和F、G_P和G以及H_P和H,A_PE_P交A_PE于上层板单元中性面的点A',同理可得点B'、C'和D'.

以点E_P和E为例,在上层单元与下层单元间的接触点处添加协调变形条件,保证几何连续性.接触点坐标需满足:

$\begin{array}{l} [r_{E}]_{z = - t / 2} = [r_{E_{P}}]_{z = t_{P} / 2} \\ {[\frac{\partial {r^{T}}_{E}}{\partial x}]}_{z = - t / 2} = {[\frac{\partial {r^{T}}_{E_{P}}}{\partial x}]}_{z = t_{P} / 2} \\ {[\frac{\partial {r^{T}}_{E}}{\partial y}]}_{z = - t / 2} = {[\frac{\partial {r^{T}}_{E_{P}}}{\partial y}]}_{z = t_{P} / 2} \\ {[\frac{\partial {r^{T}}_{E}}{\partial z}]}_{z = - t / 2} = {[\frac{\partial {r^{T}}_{E_{P}}}{\partial z}]}_{z = t_{P} / 2} \\ {[\frac{\partial^{2} {r^{T}}_{E}}{\partial z \partial z}]}_{z = - t / 2} = {[\frac{\partial^{2} {r^{T}}_{E_{P}}}{\partial z \partial z}]}_{z = t_{P} / 2} \end{array}\}$

(6)

式中:r_E和 $r_{E_{P}}$ 分别为点E和E_P的位置坐标.根据式(3),得到点E和E_P的全局坐标:

$[\begin{array}{l} r_{E} \\ \frac{\partial r_{E}}{\partial x} \\ \frac{\partial r_{E}}{\partial y} \\ \frac{\partial r_{E}}{\partial z} \\ \frac{\partial^{2} r_{E}}{\partial z \partial z} \end{array}]$ =

$[\begin{array}{l} S^{M} (l_{1} - \frac{a_{P}^{M}}{2}, l_{2} - \frac{b_{P}^{M}}{2}, - \frac{t}{2}) e^{P 1} \\ {S^{M}}_{x} (l_{1} - \frac{a_{P}^{M}}{2}, l_{2} - \frac{b_{P}^{M}}{2}, - \frac{t}{2}) e^{P 1} \\ {S^{M}}_{y} (l_{1} - \frac{a_{P}^{M}}{2}, l_{2} - \frac{b_{P}^{M}}{2}, - \frac{t}{2}) e^{P 1} \\ {S^{M}}_{z} (l_{1} - \frac{a_{P}^{M}}{2}, l_{2} - \frac{b_{P}^{M}}{2}, - \frac{t}{2}) e^{P 1} \\ {S^{M}}_{z z} (l_{1} - \frac{a_{P}^{M}}{2}, l_{2} - \frac{b_{P}^{M}}{2}, - \frac{t}{2}) e^{P 1} \end{array}]$ = T_Ee^M

(7)

$[\begin{array}{l} r_{E_{P}} \\ \frac{\partial r_{E_{P}}}{\partial x} \\ \frac{\partial r_{E_{P}}}{\partial y} \\ \frac{\partial r_{E_{P}}}{\partial z} \\ \frac{\partial^{2} r_{E_{P}}}{\partial z \partial z} \end{array}]$ =

$[\begin{array}{l} S^{P} (0,0, \frac{t_{P}}{2}) e^{P} \\ {S^{P}}_{x} (0,0, \frac{t_{P}}{2}) e^{P} \\ {S^{P}}_{y} (0,0, \frac{t_{P}}{2}) e^{P} \\ {S^{P}}_{z} (0,0, \frac{t_{P}}{2}) e^{P} \\ {S^{P}}_{z z} (0,0, \frac{t_{P}}{2}) e^{P} \end{array}]$ =

$T_{E_{P}}$ e^P

(8)

式中: ${S^{M}}_{k}$ =∂S^M/∂k, ${S^{P}}_{k}$ =∂S^P/∂k, k=x, y, z; ${S^{M}}_{z z}$ =∂²S^M/∂z∂z, ${S^{P}}_{z z}$ =∂²S^P/∂z∂z.结合式(6)~(8)可得:

$T_{E_{P}}$ e^P=T_Ee^M

(9)

对于点F_P和F、G_P和G以及H_P和H,同理可得:

$\begin{array}{l} T_{F_{P}} e^{P} = T_{F} e^{M} \\ T_{G_{P}} e^{P} = T_{G} e^{M} \\ T_{H_{P}} e^{P} = T_{H} e^{M} \end{array}\}$

(10)

矩阵形式为

$[\begin{array}{l} T_{E_{P}} \\ T_{F_{P}} \\ T_{G_{P}} \\ T_{H_{P}} \end{array}]$ e^P=

$[\begin{array}{l} T_{E} \\ T_{F} \\ T_{G} \\ T_{H} \end{array}]$ e^M

(11)

式中:[ ${T^{T}}_{E_{P}}$ ${T^{T}}_{F_{P}}$ ${T^{T}}_{G_{P}}$ ${T^{T}}_{H_{P}}$ ]^T为60×60的非奇异矩阵,因此得到下层单元与上层单元间节点的转换关系:

e^P=

${[\begin{array}{l} T_{E_{P}} \\ T_{F_{P}} \\ T_{G_{P}} \\ T_{H_{P}} \end{array}]}^{- 1} [\begin{array}{l} T_{E} \\ T_{F} \\ T_{G} \\ T_{H} \end{array}]$ e^M=T_Pe^M

(12)

2 非规则大型薄壁结构装配变形协调模型

2.1 非规则高阶曲面复合单元刚度矩阵

基于非线性连续介质力学理论,推导非规则高阶曲面单元的刚度矩阵.ANCF中的单元计算定义均基于当前构型,单元刚度矩阵随构型变化,因此可以准确描述结构的非线性大变形.定义高阶曲面单元的变形梯度为

$\frac{\partial r}{\partial x} {(\frac{\partial r_{0}}{\partial x})}^{- 1}$ =

$\frac{\partial (S e)}{\partial x} {[\frac{\partial (S e_{0})}{\partial x}]}^{- 1}$ =

$[\begin{array}{l} S_{1 x} e & S_{1 y} e & S_{1 z} e \\ S_{2 x} e & S_{2 y} e & S_{2 z} e \\ S_{3 x} e & S_{3 y} e & S_{3 z} e \end{array}] J_{0}^{- 1}$

(13)

式中:x为单元坐标;e₀为变形前单元的节点广义坐标;S_ij=∂S_i/∂j, i=1, 2, 3, j=x, y, z,S_i为单元形函数的第i行; $J_{0}^{- 1}$ 为初始构型下变形梯度的逆,且有

$J_{0}^{- 1}$

${(\frac{\partial r_{0}}{\partial x})}^{- 1}$

${[\begin{array}{l} S_{1 x} e_{0} & S_{1 y} e_{0} & S_{1 z} e_{0} \\ S_{2 x} e_{0} & S_{2 y} e_{0} & S_{2 z} e_{0} \\ S_{3 x} e_{0} & S_{3 y} e_{0} & S_{3 z} e_{0} \end{array}]}^{- 1}$

$[\begin{array}{l} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array}]$

(14)

基于相对构型,利用变形梯度推导格林-拉格朗日应变张量为

ε_G-L=

$\frac{1}{2}$ (J^TJ-I)=

$[\begin{array}{l} ε_{1} & ε_{4} & ε_{5} \\ ε_{4} & ε_{2} & ε_{6} \\ ε_{5} & ε_{6} & ε_{2} \end{array}]$

(15)

式中:

ε₁=

$\frac{1}{2}$ (

$a_{11}^{2}$ e^TS_ae+

$a_{12}^{2}$ e^TS_be+

$a_{13}^{2}$ e^TS_ce+2a₁₁a₁₂e^TS_de+2a₁₁a₁₃e^TS_ee+2a₁₂a₁₃e^TS_fe-1)

ε₂=

$\frac{1}{2}$ (

$a_{21}^{2}$ e^TS_ae+

$a_{22}^{2}$ e^TS_be+

$a_{23}^{2}$ e^TS_ce+2a₂₁a₂₂e^TS_de+2a₂₁a₂₃e^TS_ee+2a₂₂a₂₃e^TS_fe-1)

ε₃=

$\frac{1}{2}$ (

$a_{31}^{2}$ e^TS_ae+

$a_{32}^{2}$ e^TS_be+

$a_{33}^{2}$ e^TS_ce+2a₃₁a₃₂e^TS_de+2a₃₁a₃₃e^TS_ee+2a₃₂a₃₃e^TS_fe-1)

ε₄=

$\frac{1}{2}$ (a₁₁a₂₁e^TS_ae+a₁₂a₂₂e^TS_be+a₁₃a₂₃e^TS_ce+(a₂₁a₁₂+a₁₁a₂₂)e^TS_de+(a₂₁a₁₃+a₁₁a₂₃)e^TS_ee+(a₂₂a₁₃+a₁₂a₂₃)e^TS_fe)

ε₅=

$\frac{1}{2}$ (a₁₁a₃₁e^TS_ae+a₁₂a₃₂e^TS_be+a₁₃a₃₃e^TS_ce+(a₃₁a₁₂+a₁₁a₃₂)e^TS_de+(a₃₁a₁₃+a₁₁a₃₃)e^TS_ee+(a₃₂a₁₃+a₁₂a₃₃)e^TS_fe)

ε₆=

$\frac{1}{2}$ (a₂₁a₃₁e^TS_ae+a₂₂a₃₂e^TS_be+a₂₃a₃₃e^TS_ce+(a₂₁a₃₂+a₃₁a₂₂)e^TS_de+(a₂₁a₃₃+a₃₁a₂₃)e^TS_ee+(a₂₂a₃₃+a₃₂a₂₃)e^TS_fe)

(16)

式中:

$\begin{array}{l} S_{a} = (S_{1 x})^{T} S_{1 x} + (S_{2 x})^{T} S_{2 x} + (S_{3 x})^{T} S_{3 x} \\ S_{b} = (S_{1 y})^{T} S_{1 y} + (S_{2 y})^{T} S_{2 y} + (S_{3 y})^{T} S_{3 y} \\ S_{c} = (S_{1 z})^{T} S_{1 z} + (S_{2 z})^{T} S_{2 z} + (S_{3 z})^{T} S_{3 z} \\ S_{d} = (S_{1 x})^{T} S_{1 y} + (S_{2 x})^{T} S_{2 y} + (S_{3 x})^{T} S_{3 y} \\ S_{e} = (S_{1 x})^{T} S_{1 z} + (S_{2 x})^{T} S_{2 z} + (S_{3 x})^{T} S_{3 z} \\ S_{f} = (S_{1 y})^{T} S_{1 z} + (S_{2 y})^{T} S_{2 z} + (S_{3 y})^{T} S_{3 z} \end{array}\}$

(17)

应变张量具有对称性,因此应变可表示为包含6个独立应变分量的向量:

ε=

${[\begin{array}{l} ε_{1} & ε_{2} & ε_{3} & 2 ε_{4} & 2 ε_{5} & 2 ε_{6} \end{array}]}^{T}$

(18)

与应变张量对应,选取基于相对构型的应力张量,即第二类皮奥拉-基尔霍夫应力张量(σ_P-K2)表示曲面单元应力.在弹性变形阶段,σ_P-K2与ε_G-L遵循广义胡克定律,因此两者关系为

σ_P-K2=|J|J^-1σJ^-^T=D: ε _G-L

(19)

式中:D为曲面板单元的弹性矩阵,可由拉梅常数λ和μ表示为

$[\begin{array}{l} λ + 2 μ & λ & λ & 0 & 0 & 0 \\ λ & λ + 2 μ & λ & 0 & 0 & 0 \\ λ & λ & λ + 2 μ & 0 & 0 & 0 \\ 0 & 0 & 0 & μ & 0 & 0 \\ 0 & 0 & 0 & 0 & μ & 0 \\ 0 & 0 & 0 & 0 & 0 & μ \end{array}]$

(20)

式中:

λ=

$\frac{E ν}{(1 + ν) (1 - 2 ν)}$ , μ=

$\frac{E}{2 (1 + ν)}$

(21)

E和ν分别为材料的弹性模量和泊松比.根据广义胡克定律,曲面单元的弹性应变能为

$\frac{1}{2} \int_{V_{0}}$ σ_P-K2: ε _G-LdV=

$\frac{1}{2} \int_{V_{0}}$ ε^TDεdV

(22)

式中:V₀为初始构型中单元的体积.令式(22)对e求偏导,得到单元广义弹性力:

${(\frac{\partial U}{\partial e})}^{T}$ =

$\int_{V_{0}} {(\frac{\partial ε}{\partial e})}^{T}$ DεdV=Ke

(23)

从而得到与当前构型相关的刚度矩阵:

K=(λ+2μ)K₁+λK₂+4μK₃

(24)

式中:

K₁= $\int_{V_{0}} [(\frac{\partial ε_{1}}{\partial e}) ε_{1} + (\frac{\partial ε_{2}}{\partial e}) ε_{2} + (\frac{\partial ε_{3}}{\partial e}) ε_{3}]$ dV₀

K₂= $\int_{V_{0}} [(\frac{\partial ε_{2}}{\partial e} + \frac{\partial ε_{3}}{\partial e})$ ε₁+ $(\frac{\partial ε_{1}}{\partial e} + \frac{\partial ε_{3}}{\partial e})$ ε₂+ $(\frac{\partial ε_{1}}{\partial e} + \frac{\partial ε_{2}}{\partial e})$ ε₃]dV₀

K₃= $\int_{V_{0}} [(\frac{\partial ε_{4}}{\partial e}) ε_{4} + (\frac{\partial ε_{5}}{\partial e}) ε_{5} + (\frac{\partial ε_{6}}{\partial e}) ε_{6}]$ dV₀

针对构建的非规则高阶曲面复合单元模型,根据式(24)可分别得到上下层单元的刚度矩阵K_M和K_P.结合式(12),利用转换矩阵T_P将下层单元刚度转换至上层单元局部坐标系内:

${\bar{K}}_{P}$ =

$T_{P}^{T}$ K_PT_P

(25)

则非规则高阶曲面复合单元的刚度矩阵表示为

K_M-P=K_M+

${\bar{K}}_{P}$

(26)

由式(23)及(24)可知,单元刚度矩阵为节点坐标的函数,与当前构型相关,为非线性的刚度矩阵.因此,利用单元刚度矩阵能够计算非线性大变形.

2.2 非规则大型薄壁结构刚度及广义弹性力

对非规则大型薄壁结构进行离散,分别得到非规则高阶曲面单元和非规则高阶曲面复合单元的非线性刚度矩阵,复合单元的下层单元节点通过转换矩阵由上层单元表示,进而下层单元刚度矩阵可由上层单元节点求得.利用图1所示的单元间耦合方法,非规则结构的刚度矩阵可由上层单元和单层单元节点推导得到.

典型非规则大型薄壁结构的离散同图2所示.设结构I被划分为n₁个单层单元和n₂个复合单元,单层单元和上层单元共有N个节点.定义j₁的节点与结构I所有节点之间的对应关系为

$e_{j_{1}}$ =

${B^{I}}_{j_{1}}$ e_I

(27)

式中: $e_{j_{1}}$ 和e_I分别为j₁的节点和结构I中N个节点的广义坐标; ${B^{I}}_{j_{1}}$ 为 $e_{j_{1}}$ 相对于e_I的变换矩阵,且

${B^{I}}_{j_{1}}$

${[\begin{array}{l} 0_{30 \times 15 (k - 1)} & [\begin{array}{l} I_{15 \times 15} & 0_{15 \times 15} \\ 0_{15 \times 15} & I_{15 \times 15} \end{array}] & 0_{30 \times 15 (N - k - 1)} \\ 0_{30 \times 15 (l - 1)} & [\begin{array}{l} I_{15 \times 15} & 0_{15 \times 15} \\ 0_{15 \times 15} & I_{15 \times 15} \end{array}] & 0_{30 \times 15 (N - l - 1)} \end{array}]}_{60 \times 15 N}$

(28)

同理得到j₂的上层单元与结构I所有节点之间的对应关系为

$e_{j_{2}}$ =

${B^{I}}_{j_{2}}$ e_I

(29)

利用转换矩阵进行刚度组装,得到结构I的刚度矩阵为

K_I=

$\overset{n_{1}}{\sum_{j_{1} = 1}}$ (

${B^{I}}_{j_{1}}$ )^T

${K^{I}}_{j_{1}} {B^{I}}_{j_{1}}$ +

$\overset{n_{2}}{\sum_{j_{2} = 1}}$ (

${B^{I}}_{j_{2}}$ )^T

${K^{I}}_{j_{2}} {B^{I}}_{j_{2}}$

(30)

式中: ${K^{I}}_{j_{1}}$ 和 ${K^{I}}_{j_{2}}$ 分别为结构I中j₁和j₂的刚度矩阵,可由式(24)和(26)推得.

根据式(23),定义j₁和 j₂的广义弹性力向量 ${Q^{I}}_{j_{1}}$ 和 ${Q^{I}}_{j_{2}}$ ,经单元组装后得到结构I的广义弹性力向量:

Q_I=

$\overset{n_{1}}{\sum_{j_{1} = 1}}$ (

${B^{I}}_{j_{1}}$ )^T

${Q^{I}}_{j_{1}}$ +

$\overset{n_{2}}{\sum_{j_{2} = 1}}$ (

${B^{I}}_{j_{2}}$ )^T

${Q^{I}}_{j_{2}}$

(31)

2.3 非规则大型薄壁结构装配变形协调模型

在装配非规则大型薄壁结构时,一般先利用外加载荷将具有偏差零件的装配区域校正至标准位置后再连接结构,装配结束后释放外加载荷.大型结构具有柔性特征,装配体结构因弹性应变能释放而发生回弹变形.忽略由连接工艺引入的应力影响,假定装配过程中产生的应力均由外加校形载荷作用产生且为弹性力,则取消外力加载后,应力能够完全释放,导致处于标准位置下的装配体发生变形.因此,非规则大型薄壁结构的装配过程可等效为校形回弹过程.计算结构校形到标准零件位置所产生的外力,在装配结束后释放该外力,即可得到非规则大型薄壁结构装配后的结构偏差.

含初始偏差的结构I和结构II的装配过程如图5所示.AD和B'C'为装配过程中的约束边,校形力F_I和F_II分别将BC和A'D'校形至标准位置.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 非规则大型薄壁结构装配过程

Fig.5 Assembly of irregular large thin-walled structures

由校形回弹产生的广义弹性力Q_I和Q_II与相对应的F_I和F_II大小相等、方向相反.校形后,标准装配体结构I-II的刚度矩阵为

K^I-II=

$B_{I}^{T}$ K_IB_I+

$B_{I I}^{T}$ K_IIB_II

(32)

式中:B_I和B_II分别为结构I和结构II相对于结构I-II的坐标变换矩阵,参照式(28)定义.在总回弹力(F_S)作用下,标准装配体发生回弹变形,F_S与总校形力大小相等、方向相反.因此,装配过程中的装配体变形协调静力学平衡方程为

F_S=

$B_{I}^{T}$ Q_I+

$B_{I I}^{T}$ Q_II=-(

$B_{I}^{T}$ F_I+

$B_{I I}^{T}$ F_II)FS₌KI_-IIeI_-II

(33)

式中:e_I-II为装配体结构由标准位置发生回弹后的节点坐标.由式(23)、(24)和(32)可知,K_I-II为e_I-II的函数,因此式(33)为非线性方程组.利用数值迭代法求解e_I-II,根据式(3)并结合单元形函数计算得到非规则大型薄壁结构任意一点的装配后变形.

3 不同厚度分布下瓜瓣结构装配仿真分析

3.1 瓜瓣结构厚度分布及初始偏差定义

基于所提非规则高阶曲面单元及其复合单元构建重型火箭箱底瓜瓣模型,并利用非规则大型薄壁结构装配协调变形计算方法计算装配偏差,分析边缘加厚区域尺寸对瓜瓣装配变形的影响.

重型火箭箱底瓜瓣结构如图6所示.其上下底边分别为直径为900 mm和 9500 mm的1/12标准圆弧,假设结构最小厚度为18 mm,边缘区域厚度为25 mm,而边缘区域的宽度未知.将上底边加厚区域的宽度定义为从上边线沿母线方向的覆盖长度(l_U),将加厚区域底边所在圆弧对应的圆心角(φ_U)定义为其沿圆周方向的宽度尺寸,同理定义下底边加厚区域尺寸l_D和φ_D.侧边加厚区域的宽度即为母线长度(h),而沿圆周方向的宽度尺寸为其底边圆弧所对应的圆心角(φ_LR).

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 瓜瓣结构尺寸定义

Fig.6 Size of irregular thin-walled structure

瓜瓣结构的初始偏差主要由加工成形过程的回弹产生.偏差受多因素影响,在一定范围内波动.根据瓜瓣几何型面特征选取关键控制点,按照零件的偏差波动范围得到各点的随机偏差值.利用关键控制点拟合偏差型面得到含有初始偏差的瓜瓣结构,如图7所示.曲面最大偏差为50.2 mm.装配过程参照图5,其中瓜瓣I存在初始偏差,瓜瓣II则为标准结构,两者的厚度分布及尺寸参数相同,可计算装配体的装配变形.

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 带初始偏差瓜瓣结构

Fig.7 Irregular thin-walled structure with deviation

3.2 加厚区域宽度对瓜瓣装配变形影响

利用MATLAB仿真软件实现所提装配偏差计算模型的数值化,具体过程如图8所示.定义材料为2219铝合金,设E=70 GPa,设定密度和泊松比.

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 MATLAB仿真计算流程

Fig.8 Calculation procedure with MATLAB simulation

分别利用所提计算模型和Abaqus仿真软件对相同加厚区域尺寸参数的瓜瓣结构进行装配偏差仿真计算,其中,φ_LR=3°,l_U=l_D=250 mm,最大厚度为25 mm.两种方法计算得到的下底边AB装配后偏差曲线如图9所示.可知,所提装配偏差计算模型结果与Abaqus仿真结果相符,最大误差为3.2%,验证了装配偏差计算模型的准确性.而Abaqus仿真中采用三维实体单元离散结构,为保证计算准确性,模型沿厚度方向至少划分为3层单元.因此,相较于Abaqus的三维实体单元建模,所提方法能以更少单元数量进行建模仿真计算,降低了模型复杂度,同时保证了数值计算准确性.

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 本文模型与Abaqus仿真结果对比

Fig.9 Comparison of model proposed and Abaqus simulation

仿真计算不同侧边加厚区域尺寸对装配变形的影响,侧边加厚区域尺寸变化由φ_LR值表示.图10为不同φ_LR值对瓜瓣I各边装配偏差的影响.其中,ε_max为偏差最大值.可知,随着φ_LR值变化,瓜瓣结构整体仍近似于大型薄壁结构,装配变形形式受加厚区域参数变化的影响较小;受边缘加厚区域影响,底边装配后最大偏差点出现在靠近装配边端点B和C处,而装配边装配后的最大偏差点出现在靠近底边端点B处.

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 不同φ_LR值对底边和装配边装配偏差的影响

Fig.10 Assembly deviation of sides AB, CD, and BC at different values of φ_LR

图11为不同φ_LR值对瓜瓣I各边ε_max值的影响.可知,随着φ_LR值增大,各边ε_max值均呈非线性减小,结构整体刚度增加.而在相同宽度尺寸下,与上底边相比,下底边加厚区域的增加量更大,因此下底边区域的装配变形减少和刚度增大更显著,即径厚比大的区域结构,其刚度受加厚区域宽度变化影响更显著.

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 不同φ_LR值对底边和装配边装配偏差最大值的影响

Fig.11 Maximal error of sides AB, CD, and BC at different values of φ_LR

图12为l_U和l_D对瓜瓣I各边ε_max值的影响.可知,随着l_U值增大,上底边的ε_max值呈非线性减小,而下底边和装配边的ε_max值基本保持不变;随着l_D值增大,上底边的ε_max值保持不变,而下底边和装配边的ε_max值呈非线性减小.底边区域的结构刚度随对应加厚区域宽度增加呈非线性增加,而侧边区域结构刚度随l_D值增大而非线性增加.由于大型薄壁结构的尺寸较大,局部尺寸变化对对应位置的影响显著,而下底边区域比上底边区域的覆盖面积更大,所以对整体的影响更显著.

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 l_U和l_D值对底边和装配边装配偏差最大值的影响

Fig.12 Maximal error of sides AB, CD, and BC at different values of l_U and l_D

3.3 加厚区域厚度对瓜瓣装配变形影响

瓜瓣结构的焊接区域厚度对焊后结构的强度影响显著,因此根据焊接区域厚度定义边缘加厚区域厚度.图13为不同最大厚度(t_max)对瓜瓣I各边ε_max值的影响.可知,随着t_max值增大,上底边和装配边的ε_max值呈非线性减小,上底边区域刚度显著增加,厚度增加对局部(径厚比小)区域的刚度影响较大.

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13 不同t_max值对底边和装配边装配偏差最大值的影响

Fig.13 Maximal error of sides AB, CD, and BC at different values of t_max

综上可知,在不同加厚区域尺寸参数变化影响下,ε_max值均呈非线性变化,且总体趋势为随着加厚区域尺寸增大,ε_max值减小,说明装配体整体刚度增加.考虑实际瓜瓣的轻量化设计,保证结构刚度并减少加厚区域体积,选取各结果曲线拐点,即ε_max值减小量最小处的尺寸作为结构设计参数.当l_U=250 mm、l_D=500 mm、t_max=25 mm时,随着尺寸增加,刚度变化不显著.随着φ_LR值增大,整体刚度增加量呈非线性增加,考虑实际轻量化设计及满足焊接区域宽度需求(200 mm),选取φ_LR=3°,能够减小装配过程中因整体柔性特征而产生的变形,增加整体结构刚度,减轻后续装配校形难度和工作量,同时满足实际结构设计轻量化需求.

4 结语

基于绝对节点坐标法,提出非规则高阶曲面单元及其复合单元模型;构建不同厚度单元耦合方法,用以离散存在厚度差异且对边不等的非规则大型薄壁曲面结构,有效解决了由单层单元离散曲面结构产生的节点连接问题,以此描述非规则大型薄壁结构的大变形;构建非规则大型薄壁结构装配协调变形模型并解决了装配偏差计算问题.

以重型火箭瓜瓣结构为例,利用构建的变形计算模型对瓜瓣装配偏差进行仿真分析,研究瓜瓣边缘加厚区域尺寸对装配偏差和结构刚度的影响.不同区域的刚度对结构参数变化敏感度不同.径厚比大的区域,其结构刚度受加厚区域宽度变化影响显著;而径厚比小的区域(局部区域),其结构刚度受加厚区域厚度变化影响显著.根据分析结果提出合理的瓜瓣加厚区域尺寸参数组合,为重型火箭贮箱箱底设计提供指导.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

王国庆, 李曙光, 吴会强.

重型火箭贮箱大型结构制造技术现状及发展分析

[J]. 宇航材料工艺, 2014(Sup.1): 1-6.

[本文引用: 1]

WANG

Guoqing

, LI

Shuguang

, WU

Huiqiang

Status and development analyses on manufacturing technologies for large scale structures of heavy-lift launch vehicle propellant tanks

[J]. Aerospace Materials & Technology, 2014(Sup.1): 1-6.

[本文引用: 1]

[2]

CHEN

H J

, TSAI

S W

Analysis and optimum design of composite grid structures

[J]. Journal of Composite Materials, 1996, 30(4): 503-534.

DOI:10.1177/002199839603000405 URL [本文引用: 1]

[3]

OLSON

M D

, HAZELL

C R

Vibration studies on some integral rib-stiffened plates

[J]. Journal of Sound and Vibration, 1977, 50(1): 43-61.

DOI:10.1016/0022-460X(77)90550-8 URL [本文引用: 1]

[4]

王博, 田阔, 郑岩冰, 等.

超大直径网格加筋筒壳快速屈曲分析方法

[J]. 航空学报, 2017, 38(2): 220387.

[本文引用: 1]

WANG

, TIAN

Kuo

, ZHENG

Yanbing

, et al.

A rapid buckling analysis method for large-scale grid-stiffened cylindrical shells

[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(2): 220387.

[本文引用: 1]

[5]

CARRERA

Historical review of Zig-Zag theories for multilayered plates and shells

[J]. Applied Mechanics Reviews, 2003, 56(3): 287-308.

DOI:10.1115/1.1557614 URL [本文引用: 1]

[6]

D H

Layerwise theories of laminated composite structures and their applications: A review

[J]. Archives of Computational Methods in Engineering, 2021, 28(2): 577-600.

DOI:10.1007/s11831-019-09392-2 URL [本文引用: 1]

[7]

SHABANA

A A

. An absolute nodal coordinate formulation for the large rotation and deformation analysis of flexible bodies[R]. Chicago, USA: University of Illinois at Chicago, 1996.

[本文引用: 1]

[8]

MIKKOLA

A M

, SHABANA

A A

A non-incremental finite element procedure for the analysis of large deformation of plates and shells in mechanical system applications

[J]. Multibody System Dynamics, 2003, 9(3): 283-309.

DOI:10.1023/A:1022950912782 URL [本文引用: 1]

[9]

MATIKAINEN

M K

, VALKEAPÄÄ

A I

, MIKKOLA

A M

, et al.

A study of moderately thick quadrilateral plate elements based on the absolute nodal coordinate formulation

[J]. Multibody System Dynamics, 2014, 31(3): 309-338.

DOI:10.1007/s11044-013-9383-6 URL [本文引用: 1]

[10]

EBEL

, MATIKAINEN

M K

, HURSKAINEN

V V

, et al.

Analysis of high-order quadrilateral plate elements based on the absolute nodal coordinate formulation for three-dimensional elasticity

[J]. Advances in Mechanical Engineering, 2017, 9(6): 168781401770506.

[本文引用: 1]

[11]

ABBAS

L K

, RUI

X T

, MARZOCCA

Aerothermoelastic analysis of panel flutter based on the absolute nodal coordinate formulation

[J]. Multibody System Dynamics, 2015, 33(2): 163-178.

DOI:10.1007/s11044-014-9410-2 URL [本文引用: 1]

[12]

H D

, ZHAO

Z J

, YANG

, et al.

A new composite plate/plate element for stiffened plate structures via absolute nodal coordinate formulation

[J]. Composite Structures, 2020, 247: 112431.

DOI:10.1016/j.compstruct.2020.112431 URL [本文引用: 1]

重型火箭贮箱大型结构制造技术现状及发展分析

2014

... 厚度非均匀分布的非规则大型薄壁结构件具有重量轻、强度高等优点,广泛运用于航空航天、船舶等领域^[1].非规则大型薄壁结构件的几何特征造成其结构刚度的非均匀分布,在加工和装配过程中易产生变形以及转动耦合的非线性翘曲变形,导致结构存在较大偏差.结构的复杂性使得利用解析法求解其非线性变形十分困难.而在数值仿真中,为满足精度要求,传统有限元法受单元线性假设限制,离散结构需要大量单元而导致计算效率低.因此,针对非规则大型薄壁结构,构建合理的高阶单元模型以及采用有效的结构离散方法是描述其翘曲变形并提高计算效率的关键. ...

重型火箭贮箱大型结构制造技术现状及发展分析

2014

Analysis and optimum design of composite grid structures

1996

... 利用等效单元法^[2]和复合单元法^[3]可以解决壳单元离散的局限性并减少单元数量.等效单元法通过内力等效将附加厚度结构的作用等效至主体结构,但无法描述两者之间的变形耦合^[4];复合单元法将主体及附加厚度结构等效为单元叠加,通过单元间变形协调及连续性条件构建复合单元模型.针对板类结构,Carrera^[5]总结了两种常用的板壳叠层结构建模方法, 其中假定各层单元均以单层板形式变形的单一板模型忽略了柔性结构的非线性大变形特征,Layer-wise模型(一种预训练语言模型)通过在各单层板接触界面处插入高阶边界条件提高计算精度,但结构自由度的增加导致计算效率低,且该模型常用于复合材料结构建模及变形计算^[6].目前的复合单元主要基于小变形和小转动假设,在针对大变形和大旋转的非线性问题方面,Shabana等^[7-8]提出绝对节点坐标法(ANCF),用斜率坐标替代传统的转角来精确描述单元截面变形,并提出非线性弹性力和大变形的计算方法;Matikainen等^[9]通过在厚度方向引入二阶位移插值构建高阶ANCF板单元,克服了弯曲产生的泊松锁定;Ebel等^[10]通过增加单元节点消除前者产生的附加剪切锁定现象,提高了计算精度.Abbas等^[11]通过修正板单元的形函数构建变厚度的板单元模型,消除了ANCF对规则单元模型的限制,并利用模态分析以及ANSYS软件的对比结果验证了单元模型的正确性.基于复合单元法,Yu等^[12]考虑板单元间接触界面的变形协调条件,构建全参数的板/板单元模型来离散含筋板结构,消除了筋板穿透现象且具有较高的计算精度,但该单元在描述弯曲单元时会产生泊松锁闭问题.现有的ANCF单元主要用于离散厚度非均匀的规则薄壁结构或单一厚度的非规则薄壁结构,而针对厚度非均匀且形状非规则薄壁结构离散的高阶单元模型研究较少. ...

Vibration studies on some integral rib-stiffened plates

1977

超大直径网格加筋筒壳快速屈曲分析方法

2017

超大直径网格加筋筒壳快速屈曲分析方法

2017

Historical review of Zig-Zag theories for multilayered plates and shells

2003

Layerwise theories of laminated composite structures and their applications: A review

2021

1996

A non-incremental finite element procedure for the analysis of large deformation of plates and shells in mechanical system applications

2003

A study of moderately thick quadrilateral plate elements based on the absolute nodal coordinate formulation

2014

Analysis of high-order quadrilateral plate elements based on the absolute nodal coordinate formulation for three-dimensional elasticity

2017

Aerothermoelastic analysis of panel flutter based on the absolute nodal coordinate formulation

2015

A new composite plate/plate element for stiffened plate structures via absolute nodal coordinate formulation

2020

〈

〉