板状结构自发大变形问题的三维数值分析1)

张默涵 李录贤

(西安交通大学航天航空学院机械结构强度与振动国家重点实验室,西安 710049)

(西安交通大学飞行器环境与控制陕西省重点实验室,西安 710049)

引言

自然界中的物体有着丰富多彩的形貌,其中又以树叶、花朵、海带等板状结构最具代表[1-4].板状结构是指完全相同的面状结构在厚度方向堆砌而形成的厚度尺寸比面内尺寸相比较小的一类特殊三维结构[5].与一般三维结构相同的是,板状结构因内部生长或外部环境等因素也存在内部应力,并通过文献[6-7]沿叶片主脉方向等间距切割树叶的实验得到了证实.但与一般三维结构不同的是,板状结构更易通过变形模式转变释放内部应力而具有复杂的平衡构型[8-9],例如大豆的豆荚在张开后具有螺旋扭曲的形状[10-11],这正是板状三维结构的特殊之处.

从能量角度考察,板状结构的变形遵守系统总能量最小原理;从形貌角度考察,因总是选取刚度较小、易于变形的模式,板状结构常常形成更复杂的形状,这正是大自然中花朵和树叶具有迷人而复杂图案的物理机制[12].该机制可为工业设计和制造提供新的灵感:文献[13-14]通过释放柔性板状结构的预应变分别得到形状复杂的器官芯片和各类特定形状的飞行器;文献[15-16]通过外界刺激使水凝胶制成的板状结构内部发生不均匀变形,进而驱动水凝胶运动;文献[17]通过释放粘接在一起的双层板状结构因变形不同而产生的内部应力制备了纳米管.这种机制还可解释生物生命活动中的许多现象,例如团藻内细胞片的内陷[18]、人类双眼皮的形成[19]等.

自然界和人工系统中许多物体的变形特征是,在无外界载荷作用或几何约束时,结构受内部存在的不协调几何因素激发而发生自发大变形[5,20],同时产生残余应力.该问题是一个经典的非协调弹性力学问题,为此发展了三种基本理论[21]:第一种理论是将残余应力处理为一个物理场、并表征其特性;第二种理论是将变形梯度分解为生长或外部环境引起的不协调变形梯度与弹性松弛变形梯度之积;第三种理论是将应变分解为生长或外部环境引起的不协调应变与弹性应变之和.由于大变形问题的非线性和三维问题的复杂性,借助于有限元方法优势的数值分析可望成为该类问题求解的一种有效途径.

实际上,板状结构因厚度尺寸比面内尺寸明显较小,一般将其转化为薄板问题加以研究.例如,Efrati 等[20]采用曲线坐标系,以目标度量张量描述薄板局部无应力时的期望距离,度量张量g描述变形后结构的实际距离,推导了“不协调弹性理论”,并经退化建立了非欧板理论.但是,由于曲线坐标系的引入及内部应力场分布的复杂性,运用板理论可以求解的板状结构自发大变形问题仅限目标度量张量为简单圆锥曲线函数[5,20]的几个实例.

本文借助于有限元分析,对板状结构的自发大变形问题进行研究.

1 基本理论

1.1 自发大变形问题描述

对于因生长、外部环境等因素而具有内部应力的板状结构 Ω,建立如图1 所示的笛卡尔坐标系统,其中o-x-y为板状结构的中面,z为板状结构的离面方向.结构的上下面均自由;侧面边界 Γ 分为给定位移(本质) 边界条件u(x,y,z)=u0(x,y) 的边界 Γu和零外力(自然) 边界条件 σij(x,y,z)nj(x,y)=0 的边界ΓS,其中 σij为柯西应力,nj为板状结构侧面的外法向余弦.另外,板状结构内部还存在一个沿厚度不变的初始不协调应变场(x,y),它是结构发生自发大变形而产生内部应力的动因,也是本文所研究自发大变形问题的特色.对于边界 ΓS上沿厚度方向作用非零外力 σij(x,y,z)nj(x,y)=ti(x,y) 的情形,可等效为相应内部应力场产生的自发大变形问题加以研究.

图1 板状结构及坐标系统Fig.1 A plate-like structure and the coordinated system

其中Vi j为左柯西-格林变形张量.

结合柯西应力场 σij(x,y,z),结构变形能E可表示为

本文研究的自发大变形问题,是一个典型的大挠度、小应变问题,根据文献[23],材料仍服从胡克定律,本构关系可表示为

其中G=Y/(1+2ν)和λ=νY/[(1+ν)(1-2ν)] 为拉梅常数,Y和ν 为材料的杨氏模量和泊松比.

由于无外力做功,结构的总势能等同于系统的变形能,于是,结构将在变形能最小原理支配下发生自发大变形.对于这样一个非线性三维大变形问题,本文运用有限元方法加以分析.

1.2 板状结构的特性分析

1.2.1 板状结构的变形模式分析

为了说明板状结构可能的变形模式,考察面内尺寸为 2L、厚度为h的方形板状结构.假定从中切出边长为L的方形块,这样就形成一个U 型结构,如图2(a).接下来,在方形槽中插入一个材料相同,三边长为L、第四边(外边)长为L0(＞L) 的梯形块,使边长为L的三边与U 型结构的相应边粘在一起.这样,U 型结构会因L0＞L而张开,同时梯形块将受到压缩,如图2(b)所示;但是,如果板状结构足够薄,梯形块将不再保持面内的压缩变形模式,而是弯曲成如图2(c)所示的形状[20].

图2 板状结构及其变形模式Fig.2 A plate-like structure and its deformation modes

上述分析表明,存在内部应力的板状结构,具有两种可能的变形模式:一种是面内的伸缩模式,其特征是沿板厚度方向的每一个平面都相同,在厚度方向保持上下对称;另一种是离面的弯曲模式,其特征是上下表面之其一发生伸长、而另一则发生收缩,在厚度方向不再上下对称.可以看出,该结构的变形由于梯形块与正方形空槽间的几何不协调而产生,而变形模式则与结构的厚度密切相关.

1.2.2 板状结构的变形能分离

当厚度h较小时,通过引入Kirchhoff-Love 假设,板状结构可转化为板问题加以研究,例如Föppl-Von Kármán 板理论[24]、Koiter 板理论[23]等.这些板理论均表明,大变形时板状结构的变形能E可分离为伸缩变形能Es和弯曲变形能Eb之和,即

薄板弯曲时,伸缩刚度随厚度h线性变化,伸缩能Es是中面应变的二次函数,其表达式为[25]

式(6)表明,随着面内应变的增加,结构的挠度将经历一个从无到有、由小到大并逐渐占主导的过程1理论上已证明[21],只有含压缩应变问题才可诱发横向位移.,这是薄板问题之所以发生失稳的力学机制.

薄板弯曲时,弯曲刚度随厚度h三次变化[5,26],弯曲能Eb是中面弯曲应变的二次函数,表达式为[25]

需要说明的是,由于自发大变形过程中同时存在的面内伸缩变形和离面弯曲变形也与厚度h有关,将伸缩能表述为与厚度h成线性变化和将弯曲能表述为与厚度h成三次变化的说法(例如文献[5,26])是不恰当的,与文献[20]中图3 所描述的规律也不相符.实际上,当目标度量张量对应的高斯曲率为正时,板状结构的伸缩能和弯曲能均正比于厚度的2.5 次方;高斯曲率为负时,伸缩能则正比于厚度的4 次方[27].

图3 正方形板状结构的几何尺寸Fig.3 Dimensions of the square plate-like structure

仿照薄板结构的变形能分离方法,将一般厚度板状结构的总变形能E分离为

其中Em称为类伸缩变形能(简称类伸缩能),指结构随中面一起伸缩的变形能,Er则为除去Em的结构剩余变形能(简称剩余能).

利用有限元分析时,板状结构的总变形能E为每个单元的变形能之和,即

根据定义,类伸缩能Em经结构中面的变形能计算得到,即

根据式(8)～ 式(10),剩余能为

随着结构厚度h变小,板状结构的三维变形因逐渐符合Kirchhoff-Love 假设而趋近于薄板结构的变形,此时,三维角度的Em和Er相应地就与薄板角度的Es和Eb逐渐接近.

2 自发大变形的三维分析方法

板状结构的自发大变形行为虽然只在薄板情形下得以凸显,但考虑到薄板结构大变形理论求解的复杂性,以及薄板结构可视为厚度h逐渐变薄的中等厚度三维板状结构的客观事实,本文借助于三维大变形有限元分析,提出板状结构自发大变形问题的求解方法.

第1 节的分析已表明,板状结构的变形能由类伸缩能Em和剩余能Er两部分组成.本节将通过受外力作用板状结构的三维大变形数值计算,对板状结构的变形能构成予以定量分析,从能量角度提出屈曲失稳条件,进而揭示板状结构失稳的力学机制.

2.1 板状结构的三维大变形分析

如图3 所示,考虑边长为600 mm 的正方形板状结构,假定前后两边自由、左右两边简支,左端沿x方向作用沿厚度均匀分布的压缩载荷q.选用具有三维大变形分析功能的二次缩减积分单元C3D20R,经收敛性验证后确定沿厚度方向均匀划分11 层,面内划分成20×20 格,共计4400 个单元.分析时取q=228 N/mm,材料的杨氏模量Y=2.55 GPa、泊松比ν=0,计算得到不同厚度时板状结构中心位置o点处(参考图3)的横向位移uz,如图4所示.

由图4可以看出,当厚度h＞ 8.08 mm 时,横向位移uz很小,结构处于单一的面内伸缩变形模式;当厚度h略小于8.08 mm 时,横向位移陡然增大,并因厚度的不同或正或负,出现分岔现象,表明结构的变形中此时增添了新的离面弯曲模式.

图4 方板中心位置 o 处横向位移 uz随板厚度 h 的变化Fig.4 Variation of uzat opoint with thicknessh

2.2 变形能计算及屈曲失稳的能量条件

根据2.1 节的大变形有限元计算,利用式(9)可得到结构的总变形能E,利用式(10)可得到板状结构的类伸缩能Em,进而利用式(11)得到剩余能Er.这样,图3 所示问题中类伸缩能Em和剩余能Er随厚度的变化如图5 所示.

图5 方形板状结构能量构成随厚度的变化Fig.5 Variation of the strain energy with thickness for a square platelike structure

从图5可以看出,在8.01 mm 至8.10 mm 较小的厚度变化范围内,结构的总变形能大小及构成却发生了大幅变化.当厚度h＞ 8.08 mm 时,板状结构的剩余能Er为零,类伸缩能Em在总变形能中占主导地位,称8.08 mm 为该结构的临界厚度hcr,也就是说,厚度h＞hcr时,板状结构内部不会存在弯曲变形能.

分析图5 还发现,板状结构存在一个类伸缩能Em与剩余能Er相同的交叉点,称所对应厚度为变形模式转变厚度htr,在此例中其值为8.076 mm.从三维角度,模式转变厚度htr表征了剩余能Er和类伸缩能Em相等时的厚度,也就是说,在厚度由大变小过程中板状结构具有类伸缩能与剩余能大小关系发生反转的特点[21,28-29],在变形模式上表现为以面内伸缩模式为主导变成为以离面弯曲模式为主导.

与文献[20-21]类似,基于模式转变厚度时所表现出的能量关系,本文提出板状结构屈曲失稳的能量条件为

对照附录中图A1 基于薄板经典稳定性理论中qcr=228 N/mm 对应的厚度h=8.076 mm,本文基于三维大变形数值分析得到的转变厚度htr与之完全吻合,侧证了本文分析方法及式(12)能量条件的正确性.

本节工作表明,三维大变形有限元分析是定量确定板状结构模式转变厚度htr的有效途径.

3 自发大变形实例分析

运用温度场作用下结构的热胀冷缩效应,可有效模拟结构由于生长或外部环境引起的内部不协调变形(x,y)[30].本节运用不同种类温度场变化,通过三维大变形有限元计算,分析板状结构的不同自发大变形行为.

3.1 正方形板状结构的自发大变形

如图6 所示,考察600 mm×600 mm 正方形板状结构.取材料的杨氏模量为2.55 GPa、泊松比为0.从中划出一个深度300 mm、宽度为 2a的小矩形区域,设该区域作用一沿y方向线性非均匀变化的温度场 ΔT(参考图6),并假定材料仅在x方向具有非零的热胀系数 α=7.5 ×10-4K-1.这样,正方形板状结构将在外部边界自由、内部不协调变形作用下产生内部应力场,并因之诱发自发变形.此时,x方向正应力σx的典型分布如图7 所示,可以看出,结构内部既有受压(蓝色)区域,也有受拉(红色)区域.

同样选取二次缩减积分单元C3D20R,它也具有分析温度场作用下三维大变形问题能力.结构沿厚度方向均匀划分11 层;将温度作用的区域面划分为15×15 格,共2475 个单元;其它区域在中面内划分320 格,共3520 个单元.

首先分析厚度h对该正方形板状结构屈曲失稳的影响.为此,假定温度场作用的半宽度a=150 mm,令厚度h在[2.5 mm,10 mm]间每隔0.5 mm 共16 种变化,分别进行三维大变形有限元计算.

运用与2.2 节相同的技术计算结构的变形能及其构成,它们随厚度的变化如图8 所示.

图8 局部温度场作用时正方形板状结构能量随厚度的变化Fig.8 Variation of the strain energy with thickness due to partial change in temperature

由图8可以看出,当板的厚度h＞ 8.5 mm 时,结构的总能量几乎都是类伸缩能,剩余能可忽略不计,此时受温度作用区域的变形模式为纯粹的面内伸缩,此时,图6 所示问题的临界厚度hcr=8.5 mm.当厚度h在8.5 mm 至6.89 mm 之间变化时,总能量中剩余能的占比逐渐增加,板状结构处于面内伸缩模式和离面弯曲模式的混合变形阶段,此时,弯曲特征还不十分明显(参考图9(a)中的h=7.5 mm 情形).当厚度h＜ 6.89 mm 时,剩余能较类伸缩能已占据主导地位,离面弯曲特征已较为明显(参考图9(b)中的h=6 mm 情形),此时,图6 所示问题的模式转变厚度htr=6.89 mm.

图9 局部温度场作用下板状结构沿厚度的变形与应力分布Fig.9 Deformed shape and stress distribution in thickness direction due to change in temperature

由图8可进一步看出,由于与图5 外部压力q作用下结构内仅具有均匀单向压应力的情形不同(参考图7),内部不协调变形引起的结构变形能及构成变化更为复杂,临界厚度hcr与模式转变厚度htr不再接近,两者相差18.9%,这也是本文采用式(12)作为屈曲失稳条件的主要原因.

其次分析温度场作用的半宽度a对该结构屈曲失稳的影响.为此,假定板的厚度h=7 mm,令半宽度a在[100 mm,200 mm]间每隔10 mm 共11 种变化,分别进行三维大变形有限元计算.

运用与2.2 节相同的技术计算结构的变形能及其构成,它们随局部温度场半宽度a的变化如图10所示.从图可以看出,温度作用范围a所代表的不协调变形因素,导致结构内变形能构成的变化,也会引起结构的屈曲失稳,与改变结构厚度h的情形类似.

图10 正方形板状结构能量构成随局部温度场半宽度 a 的变化Fig.10 Variation of the strain energy with the half-width of local temperature in a square plate-like structure

3.2 圆形板状结构的自发大变形

图11 所示为半径50 mm 的圆形板状结构;材料为仿树叶类PA 水凝胶弹性材料[30-31],杨氏模量为0.35 MPa、泊松比为0.25.圆形板状结构处于沿径向r幂律变化的温度场 ΔT=β(r/R)n中,并假定结构仅具有沿周向的非零热胀系数 α=1 ×10-3K-1.这样,温差幅值β和指数n两个参数的变化在圆形板状结构内将产生不同大小和分布的内部应力场,升温和降温时周向正应力 σθ的典型分布如图12 所示,可以看出,整个结构中既有受压区域(蓝色)、也有受拉区域(红色).该结构在边界自由情形下将发生自发变形.

图11 处于非均匀变化温度场中的圆形板状结构Fig.11 A circular plate-like structure subjected to a power-law change in temperature

图12 圆形板状结构的典型周向正应力分布Fig.12 Typical distributions of circumferentially normal stresses in the circular plate-like structure

取结构的1/4 进行建模.仍选用二次缩减积分单元C3D20R,整体结构沿厚度方向均匀划分11 层,面内划分442 格,共4862 个单元.运用与2.2 节相同的技术计算结构的变形能及其构成,并仿图8 得到模式转变厚度htr,进而研究温差幅值β和指数n两个参数对htr的影响.

3.2.1 温差幅值β对htr的影响

给定指数n=0.5,1,2和5,在[-5°C,5°C]间每隔1°C 共11 组β值进行计算,其中β的正负分别表示对应点处为升温或降温.图13 为给定n时圆形板状结构模式转变厚度htr随β的变化关系.

由图13可以看出,模式转变厚度htr随β绝对值的增加而增大,但由于升降温时内应力分布的差异(参考图12),其效果并不对等.以指数n=2 为例,β=5°C 时的htr是β=1°C 时的2.35 倍,而β=-5°C 时的htr只是β=-1°C 时的2.25 倍.总之,温差幅值β的绝对值对该问题模式转变厚度htr的影响显著.

图13 幅值β对模式转变厚度 htr 的影响关系Fig.13 Influence of amplitude β onhtr

3.2.2 指数n对htr的影响

给定β=-5°C,-3°C,3°C和5°C,在[0.5,5]间每隔0.5 共取10 组非零n值进行计算,图14 所示为给定β时 圆形板状结构的模式转变厚度htr随n的变化关系.

图14 指数 n 对模式转变厚度 htr 的影响关系Fig.14 Influence of index number n onhtr

由图14可以看出,温差幅值β正负的不同,模式转变厚度htr随n的变化规律明显不同.对于同等温度幅值,升温时结构更容易发生屈曲失稳.温差幅值β为正时,模式转变厚度htr随n先略微增大后又略微减小,n=1 时htr为最大;温差幅值β为负时,模式转变厚度htr随n单调减小;但总体影响幅度并不十分显著1进一步计算表明,指数h 对屈曲波数具有显著影响..例如,对于β=3°C 的情形,n=5 时的htr比n=0.5 时减小了23.38%;对于β=-3°C,n=5 时的htr比n=0.5 时减小了35.29%.总之,与幅值β的影响相比,指数n对模式转变厚度htr的影响相对较小.

本节所研究问题中转变厚度htr的变化规律表明,板状结构因内部不协调变形诱发的屈曲失稳是一个复杂的自发大变形过程,与产生内部不协调变形的各个因素密切相关.

4 结论

本文在板状结构中引入沿厚度不变的不协调初始应变,通过三维大变形有限元分析,将结构的变形能分离为类伸缩能和剩余能,提出板状结构屈曲失稳的能量条件,建立了板状结构自发大变形问题的研究方法,从三维大变形角度,揭示了板状结构屈曲失稳的物理机制.根据本文工作,可得到以下结论.

(1)在板状结构大变形过程中,当剩余能从零增加到与类伸缩能相等时,板状结构将发生由面内伸缩模式为主导到离面弯曲模式为主导的屈曲失稳,经典的载荷屈曲条件是本文能量屈曲条件在外力作用下的特殊情形.

(2)对于受外力作用的板状结构,由于结构内部压应力场分布和变化相对简单,大变形时剩余能将从零快速增加至超越类伸缩能,使得结构发生骤然屈曲失稳现象.

(3)在内部不协调变形因素作用下,由于板状结构内部压应力分布的不均匀性、甚至拉应力的存在等多种原因,大变形时剩余能经历一个从小到大量变、再到超越类伸缩能质变的逐步变化过程,因而,结构一般不会发生骤然的屈曲失稳.

(4)不协调因素是板状结构在无外部约束下仍存在复杂内部应力分布的主要原因,所诱发的自发大变形(包括屈曲失稳、后屈曲等)是树叶、花朵等板状生物结构形成丰富形貌的物理机制.

本文借助三维有限元方法,数值分析了板状结构在多种因素作用下的大变形过程,重点考察了该类结构的屈曲失稳行为,揭示了不协调因素对转变厚度的影响规律.实际上,基于本文的三维大变形有限元方法,还可进一步研究板状结构在不协调因素作用下的后屈曲行为及复杂结构形貌,这些内容的特点是运用微分几何的曲面论知识分析板状结构中面在经历自发大变形后的形貌特征,是本文作者正在开展的工作.

附录

对于前后两边自由、左右两边简支、左端沿x方向施加沿厚度均匀分布压缩载荷q的正方形板状结构,薄板理论给出的屈曲载荷计算公式为[32]

式中,qcr为屈曲载荷,L为正方形板的边长.

对于图3 所示的问题,L=600 mm,Y=2.55 GPa,v=0,屈曲载荷值qcr随厚度h的变化如图A1 所示.特别地,经式(A1) 计算,qcr=228 N/mm 对应的厚度为h=8.076 mm.