文章快速检索     高级检索
  重庆大学学报  2016, Vol. 39 Issue (5): 82-89  DOI: 10.11835/j.issn.1000-582X.2016.05.011 RIS(文献管理工具)
0

引用本文 

韩志林, 程长征, 盛宏玉. 基于状态空间法的叠层材料热力分析[J]. 重庆大学学报, 2016, 39(5): 82-89. DOI: 10.11835/j.issn.1000-582X.2016.05.011.
HAN Zhilin, CHENG Changzheng, SHENG Hongyu. Analysis of thermal stress in laminated material based on state space method[J]. Journal of Chongqing University, 2016, 39(5): 82-89. DOI: 10.11835/j.issn.1000-582X.2016.05.011. .

基金项目

国家自然科学基金资助项目(11372094)

通信作者

程长征(联系人), 男, 合肥工业大学教授, 博士生导师, (E-mail)changzheng.cheng@hfut.edu.cn

作者简介

韩志林(1990-), 男, 合肥工业大学博士研究生, 主要从事计算力学研究, (E-mail)hanzhilin1990@163.com

文章历史

收稿日期: 2016-04-10
基于状态空间法的叠层材料热力分析
韩志林, 程长征, 盛宏玉     
合肥工业大学 土木与水利工程学院, 合肥 230009
摘要: 通过引入热力问题的应力应变关系,建立了求解叠层材料温度应力的半解析状态空间法。该方法沿叠层材料水平方向用有限元法离散位移、应力分量,在竖向采用状态空间法求解,得到叠层材料各交界面上的位移、应力分量。为验证该方法的有效性,针对某一给定温度进行热力学分析,并与有限元结果进行对比。算例表明,该方法计算精度高、计算量小。
关键词: 叠层材料    状态空间法    热弹性力学    有限元法    
Analysis of thermal stress in laminated material based on state space method
HAN Zhilin , CHENG Changzheng , SHENG Hongyu     
School of Civil Engineering, Hefei University of Technology, Hefei 230009, P.R.China
Supported by National Natural Science Foundation of China(11372094)
Abstract: A state space finite element method for 2D thermal elasticity of laminated materials was deduced by introducing stress-strain relations of thermal elastic problems. Finite element approximation was introduced for the variations of displacements and stresses on the horizontal interfaces while the through thickness distributions of the variations of displacements and stresses were obtained by utilizing state space method. The displacement and stress at the interface of laminated materials deduced by the present method were illustrated. It is shown that the present method has high calculation accuracy and low computational cost by comparing with the finite element method.
Key Words: laminated material    state space method    thermal elasticity    finite element method    

在复合材料的力学计算中,通常把复合材料简化成叠层材料进行分析。状态空间法在分析叠层材料时具有计算精度高、计算量小的特点。范家让[1]首次提出了状态空间法并将其应用在板壳问题中。状态空间法主要包括输入变量、传递矩阵、输出变量三部分,其中输入变量为叠层材料各交界面上已知的位移、应力组成的向量,输出变量为各界面未知的位移、应力分量,传递矩阵由材料属性决定。在以后的研究过程中,受到有限元思想的启发,Sheng[2]及Ye[3]等提出了与有限元方法相结合的状态空间法,该方法沿叠层材料水平方向用有限元法离散位移、应力分量,而在竖向采用状态空间法求解。当沿叠层材料水平向采用样条插值法[4],在竖向采用状态空间法时,同样可以准确计算叠层材料的弹性力学问题。由于状态空间法中输入变量与输出变量间关系简洁明了,因此该方法可以作为传递方法,与其他计算方法结合,例如,Benedetti等[5]将状态空间法与快速边界元法结合,求得了压电粘结材料的应力场。

无论在建筑工业还是机械工业中,状态空间法在分析热传导问题时都得到了广泛应用[6]。在分析瞬态热传导问题时,时域分析是状态空间法求解的关键。目前,主要有两种时域分析方法:差分格式法[7-8]以及拉普拉斯变换法[9]。差分格式法较为直观,其精度取决于时间步长的选择。以上两种方法结合状态空间法均能准确计算瞬态热传导问题,但在建筑工业范围内,求解墙体热传导的方法除了状态空间法,还有直接求根法及频域回归法。经学者研究发现,状态空间法相较于频域回归法,其准确性略差[10]

由于受到边界的约束,在热传导过程中会产生温度应力,当该应力过大时会导致结构失效。状态空间法在温度应力计算中同样得到广泛应用[11],例如,李家宇等[12]就利用状态空间法求得了四边简支压电热弹性层合正交双曲壳的精确解。另外,Youssef等[13]提出了一种广义热弹性力学理论并将其与状态空间法结合求解了一维热弹性力学问题。后来又有其他学者扩大了该理论的应用范围,例如,Ezzat等[14]进一步求得了磁热弹问题的数值解。而以上的热应力分析大多是针对无限大结构,对于有一定约束条件的有限尺寸结构,则可以延续状态空间有限元思想[2],与热弹性力学的基本方程结合在一起,利用状态空间法分析温度应力。

1 二维模型变分原理

考虑如图 1所示由M层各向同性材料组成的叠层材料,每层材料的弹性模量E、泊松比μ以及热膨胀系数α可以不同,图中x向和y向分别表示横向和纵向,横向和纵向的位移分量为uvLH分别为长度、高度,ΔT为温度荷载,两边虚线表示结构只受y向约束。叠层材料几何形状及荷载均为对称。截取左边一半用状态空间法分析,产生边界条件为:x=0边的位移v=0,应力σxx=0;x=L/2边的位移u=0,应力σxy=0,其中σxxσxy分别表示x向正应力和切应力。

图 1 直角坐标系中的叠层材料 Figure 1 Laminated material in the rectangular coordinate system

根据Hellinger-Reissner变分原理[15],对于叠层材料中的任意一层可以建立变分方程

$\mathop {\smallint \smallint }\limits_\varOmega \delta {\boldsymbol{\varsigma} ^{\rm{T}}}\left[{{\boldsymbol{E}^{\rm{T}}}\left( \nabla \right)\boldsymbol{u}-\boldsymbol{\varepsilon} } \right]{\rm{d}}\varOmega - \mathop {\smallint \smallint }\limits_\varOmega \delta {\boldsymbol{u}^{\rm{T}}}\left[{\boldsymbol{E}\left( \nabla \right)\boldsymbol{\sigma} + f} \right]{\rm{d}}\varOmega \\ - \int_{{B_u}} {\delta {\boldsymbol{P}^{\rm{T}}}\left( {\boldsymbol{u} -\boldsymbol{\bar u}} \right){\rm{d}}S{\rm{ + }}} \int_{{B_\varsigma }} {\delta {\boldsymbol{u}^{\rm{T}}}\left( {{\boldsymbol{P}_{\rm{S}}} -{{\boldsymbol{\bar P}}_{\rm{S}}}} \right){\rm{d}}S = 0, } $ (1)

式中:u为位移分量uv组成的向量; σε分别表示应力、应变向量; BuBσ分别表示已知的位移边界和应力边界;$\boldsymbol{\bar u}, {{\boldsymbol{\bar P}}_{\rm{S}}}$分别表示边界上已知的位移和应力条件,${\boldsymbol{E}^{\rm{T}}}\left( \nabla \right)$为微分算子[2]。若位移、应力场均满足边界条件,则式(1)中的$\boldsymbol{u} = \boldsymbol{\bar u}, {\boldsymbol{P}_{\rm{S}}} = {{\boldsymbol{\bar P}}_{\rm{S}}}$。因此,式(1)可以等价表示成

$\mathop {\smallint \smallint }\limits_\varOmega \delta {\boldsymbol{u}^{\rm{T}}}\left[{\boldsymbol{E}\left( \nabla \right)\boldsymbol{\sigma} + \boldsymbol{f}} \right]{\rm{d}}\varOmega = 0, $ (2a)
$\mathop {\smallint \smallint }\limits_\varOmega \delta {\boldsymbol{\sigma} ^{\rm{T}}}\left[{\boldsymbol{\varepsilon}-{\boldsymbol{E}^{\rm{T}}}\left( \nabla \right)\boldsymbol{u}} \right]{\rm{d}}\varOmega = 0, $ (2b)

以上2个方程分别表示平衡方程和几何方程。

文中方法仅在叠层材料的上边划分单元,即在第一层材料的上边划分N-1个线性单元,产生N个节点。e1ekeN-1分别表示第1个、第k个、第N-1个单元,如图 1所示。引入沿x轴方向有等参元性质的单元,节点的位移、应力参量仅为y的函数,现以第k个单元为例进行说明。第k个单元的位移和应力可以采用如下插值形式

${u^k} = \sum\limits_{i = 1}^n {N_i^k\left( \xi \right)u_i^k\left( y \right), {\upsilon ^k} = \sum\limits_{i = 1}^n {N_i^k\left( \xi \right)\upsilon _i^k\left( y \right), } } $ (3a)
$\sigma _{xy}^k = \sum\limits_{i = 1}^n {N_i^k\left( \xi \right)\sigma _{xyi}^k\left( y \right), \sigma _{yy}^k = \sum\limits_{i = 1}^n {N_i^k\left( \xi \right)\sigma _{yyi}^k\left( y \right), \sigma _{xx}^k = \sum\limits_{i = 1}^n {N_i^k\left( \xi \right)\sigma _{xxi}^k\left( y \right)} } } $ (3b)

式中:ξ为局部坐标值;Nik(ξ)表示第i个形函数。uik(y)、σxyik(y)分别表示节点的位移、应力参量,且仅为y的函数;n为插值节点数,本文中n=2,表示采用线性单元。

另外,温度应力问题中的物理方程为

$\boldsymbol{\varepsilon} = \boldsymbol{S\sigma} + \boldsymbol{J}$ (4)

式中:S为各向同性材料的柔度矩阵,J=[αT  αT  0]T。将式(4)代入式(2b),并将其拆分为平面问题表达式(5)与反平面问题表达式(6):

$\mathop {\smallint \smallint }\limits_\Omega {\left\{ \begin{array}{l} \delta {\sigma _{xy}}\\ \delta {\sigma _{yy}} \end{array} \right\}^{\rm{T}}}\left( {\frac{\partial }{{\partial y}}\left\{ \begin{array}{l} u\\ v \end{array} \right\}- \left[{\begin{array}{*{20}{c}} 0&{-\frac{\partial }{{\partial x}}}&{{S_3}}&0\\ 0&0&0&{{S_1}} \end{array}} \right]\left\{ \begin{array}{l} u\\ v\\ {\sigma _{xy}}\\ {\sigma _{yy}} \end{array} \right\} + \left[\begin{array}{l} 0\\ -{S_2} \end{array} \right]\left\{ {{\sigma _{xx}}} \right\} + \left[\begin{array}{l} 0\\ -\alpha T \end{array} \right]} \right){\rm{d}}\varOmega = 0, $ (5)
$\mathop {\smallint \smallint }\limits_\Omega {\left\{ {\delta {\sigma _{xx}}} \right\}^{\rm{T}}}\left( {{S_1}{\sigma _{xx}}- \left[{\begin{array}{*{20}{c}} {\partial /\partial x}&0&0&{-{S_2}} \end{array}} \right]\left\{ \begin{array}{l} u\\ v\\ {\sigma _{xy}}\\ {\sigma _{yy}} \end{array} \right\} + \alpha T} \right){\rm{d}}\varOmega = 0.$ (6)

将平面问题表达式(5)与(2a)合并

$\begin{array}{l} \mathop {\smallint \smallint }\limits_\Omega {\left\{ \begin{array}{l} \delta u\\ \delta v\\ \delta {\sigma _{xy}}\\ \delta {\sigma _{yy}} \end{array} \right\}^{\rm{T}}}\left( {\left[ {\begin{array}{*{20}{c}} 0&0&1&0\\ 0&0&0&1\\ 1&0&0&0\\ 0&1&0&0 \end{array}} \right]\frac{\partial }{{\partial y}}\left\{ \begin{array}{l} u\\ v\\ {\sigma _{xy}}\\ {\sigma _{yy}} \end{array} \right\} - } \right.\\ \left. {\left[ {\begin{array}{*{20}{c}} 0&0&0&0\\ 0&0&{ - \frac{\partial }{{\partial x}}}&0\\ 0&{ - \frac{\partial }{{\partial x}}}&{{S_3}}&0\\ 0&0&0&{{S_1}} \end{array}} \right]\left\{ \begin{array}{l} u\\ v\\ {\sigma _{xy}}\\ {\sigma _{yy}} \end{array} \right\} + \left[ {\begin{array}{*{20}{c}} {\frac{\partial }{{\partial x}}}\\ 0\\ 0\\ { - {S_2}} \end{array}} \right]\left\{ {{\sigma _{xx}}} \right\} + \left[ {\begin{array}{*{20}{c}} 0\\ 0\\ 0\\ { - \alpha } \end{array}} \right]T} \right){\rm{d}}\Omega = 0, \end{array}$ (7)

式(5)~(7)中的T表示每一层材料的温度;S1S2S3为柔度矩阵S中的元素。将每一层材料上边用线性单元离散,将式(3)代入到式(6)(7)中得到

$\int_y {{{\left\{ {\boldsymbol{\delta s}} \right\}}^{\rm{T}}}\left( {\boldsymbol{Ds}-\boldsymbol{E}\left\{ \begin{array}{l} \boldsymbol{p}\\ \boldsymbol{q} \end{array} \right\} + \boldsymbol{GT}} \right){\rm{d}}y{\rm{ = 0, }}} $ (8)
$\int_y {{{\left\{ \begin{array}{l} \delta \boldsymbol{p}\\ \delta \boldsymbol{q} \end{array} \right\}}^{\rm{T}}}\left( {\boldsymbol{A}\frac{{\rm{d}}}{{{\rm{d}}y}}\left\{ \begin{array}{l} \boldsymbol{p}\\ \boldsymbol{q} \end{array} \right\}-B\left\{ \begin{array}{l} \boldsymbol{p}\\ \boldsymbol{q} \end{array} \right\} + \boldsymbol{Cs} + \boldsymbol{HT}} \right){\rm{d}}y = 0, } $ (9)

式(8)(9)中pT=[uT(y)vT(y)],qT=[σxyT(y)σyyT(y)],sT=[σxxT(y)],T则表示离散单元后每个节点的温度。向量的排列形式按照节点号升序排列,例如,u(y)=[u1(y), u2(y), …, uN(y)]TABCDEGH表示常系数矩阵。

图 1所示,将已知边界条件代入式(8)(9)可等价写为

${\boldsymbol{D}_f}{\boldsymbol{s}_f} = {\boldsymbol{E}_f}\boldsymbol{R}-{\boldsymbol{G}_f}\boldsymbol{T}, $ (10)
${\boldsymbol{A}_f}\frac{{\rm{d}}}{{{\rm{d}}y}}\boldsymbol{R} = {\boldsymbol{B}_f}\boldsymbol{R}-{\boldsymbol{C}_f}{\boldsymbol{s}_f}-{\boldsymbol{H}_f}\boldsymbol{T}, $ (11)

式中:R=[pfT qfT]Tpfqfsf表示每一层材料上、下边的未知状态量。常数矩阵AfBfCfDfEfGfHf表示经过代入已知边界状态量后推导出的新的常数矩阵。T表示划分单元后所有节点的温度。联立式(10)、(11)并消去平面内应力分量sf得到了如下方程

${\boldsymbol{A}_f}\frac{{\rm{d}}}{{{\rm{d}}y}}\boldsymbol{R} = {\boldsymbol{K}_f}\boldsymbol{R} + {\boldsymbol{B}_f}$ (12)

式中:Kf=(Bf-CfDf-1Ef),Bf=(CfDf-1Gf-Hf)T

在代入边界条件时应该注意到,由于节点1处的σyy未知,状态未知量中应力未知量的个数比位移未知量的个数多1个。为了消去节点1的σyy,引入温度应力的通式

${\sigma _{yy}} = \frac{{\mu E}}{{1-{\mu ^2}}}{\varepsilon _{xx}} + \frac{E}{{1-{\mu ^2}}}{\varepsilon _{yy}}-\frac{{E\alpha T}}{{1 - \mu }}.$ (13)

在如图 1所示的问题中,由于x=0边仅受竖向约束,横向自由,因此该边界上的应变εyyεxx均为零,代入式(13)得

${\sigma _{yy}} =-\frac{{E\alpha {T_1}}}{{1-\mu }}, $ (14)

将式(14)代入式(12)中,形成任意一层的非齐次状态方程

$\frac{{\rm{d}}}{{{\rm{d}}y}}\boldsymbol{\tilde R} = {{\boldsymbol{\tilde N}}_f}\boldsymbol{\tilde R} + {{\boldsymbol{\tilde M}}_f}, $ (15)

式中:${{\boldsymbol{\tilde N}}_f} = \boldsymbol{\tilde A}_f^{-1}-{{\boldsymbol{\tilde K}}_f}$,其中的${\boldsymbol{\tilde R}}$$\tilde A_f^{-1}、{{\tilde K}_f}$分别表示消去多余未知量σyy后推导出的未知状态量和常数矩阵;Mf不仅与式(12)中的Bf有关,还与节点1处的温度有关。

2 状态方程的求解

对于叠层材料中的第i层,求解温度应力问题得出的非齐次状态方程为

$\frac{{\rm{d}}}{{{\rm{d}}y}}{\boldsymbol{R}_i} = {\boldsymbol{N}_i}{\boldsymbol{R}_i} + {\boldsymbol{M}_i}, $ (16)

由于在分析过程中温度荷载为已知量,则Mi为常数列阵。根据式(16)的解可得

$\left\{ \begin{array}{l} {\boldsymbol{R}_i}\left( y \right) = {\boldsymbol{D}_i}\left( y \right){\boldsymbol{R}_i}\left( {{y_{i- 1}}} \right) + {\boldsymbol{H}_i}\left( y \right), y \in \left[{{y_{i-1}}, {y_i}} \right], \\ {\boldsymbol{D}_i}\left( y \right) = {{\rm{e}}^{{\boldsymbol{N}_i}\left( {y -{y_{i -1}}} \right)}}, \\ {\boldsymbol{H}_i}\left( y \right) = \int_{{y_{i -1}}}^y {{{\rm{e}}^{{\boldsymbol{N}_i}\left( {y - \tau } \right)}}{\boldsymbol{M}_i}{\rm{d}}\tau {\rm{.}}} \end{array} \right.$ (17)

y=yi时,有

$\left\{ \begin{array}{l} {\boldsymbol{R}_i}\left( {{y_i}} \right) = {\boldsymbol{D}_i}\left( {{y_i}} \right){\boldsymbol{R}_i}\left( {{y_{i-1}}} \right) + {\boldsymbol{H}_i}\left( {{y_i}} \right), \\ {\boldsymbol{D}_i}\left( {{y_i}} \right) = {{\rm{e}}^{{N_i}{h_i}}}, \\ {\boldsymbol{H}_i}\left( {{y_i}} \right) = \int_{{y_{i-1}}}^{{y_i}} {{{\rm{e}}^{{N_i}\left( {{y_{i-\tau }}} \right)}}{\boldsymbol{M}_i}{\rm{d}}\tau {\rm{, }}} \end{array} \right.$ (18)

式中:hi表示第i层材料的厚度。对于式(18)中第2式,采用泰勒级数计算其近似值,截取10阶

${{\rm{e}}^{\boldsymbol{A}t}} = \boldsymbol{I} + \boldsymbol{A}t + \frac{{{\boldsymbol{A}^2}}}{{2!}}{t^2} + \frac{{{\boldsymbol{A}^3}}}{{3!}}{t^3} + \cdots + \frac{{{\boldsymbol{A}^{10}}}}{{10!}}{t^{10}}, $ (19)

对于式(18)中的第3式,采用六点高斯积分计算其数值解。将式(18)中的第2、第3式代入第1式,经整理得到线性方程组

${\boldsymbol{A}_i}{\boldsymbol{x}_i} = {\boldsymbol{b}_i}, $ (20)

式中:Ai为第i层材料的系数矩阵; bi为式(18)第1式中已知参量组成的列阵; xi为第i层材料上、下界面的未知位移及应力分量。

对第i+1层材料进行如上所述的计算后,同样可以形成线性方程组

${\boldsymbol{A}_{i + 1}}{\boldsymbol{x}_{i + 1}} = {\boldsymbol{b}_{i + 1}}.$ (21)

利用各交界面上的位移关系以及应力关系可得

$\left\{ {\boldsymbol{p}_f^{i + 1}\left( {{y_i}} \right)} \right\} = \left\{ {\boldsymbol{p}_f^i\left( {{y_i}} \right)} \right\}, \left\{ {\boldsymbol{q}_f^{i + 1}\left( {{y_i}} \right)} \right\} = \left\{ {\boldsymbol{q}_f^i\left( {{y_i}} \right)} \right\}, $ (22)

式中:{pfi+1(yi)}及{pfi(yi)}分别表示第i+1层材料上端界面及第i层材料下端界面的位移分量; {qfi+1(yi)}及{qfi(yi)}分别表示第i+1层材料上端界面及第i层材料下端界面的应力分量,即,Ri(yi)=Ri+1(yi)。将各层材料形成的线性方程组组合形成叠层材料的整体线性方程组

$\boldsymbol{A}x = \boldsymbol{b}, $ (23)

解得的x即叠层材料各界面上所有未知的位移及应力分量。

3 算例

图 2为由3层各向同性材料组成的叠层材料,长L=2 m,高H=0.3 m,每一层材料高h=0.1 m。x=0、x=L边仅受沿y轴的位移约束,如图 2中虚线所示。沿y轴正方向从上至下的材料参数分别为弹性模量E1=2.1×1011 Pa,E2=1.05×1011 Pa,E3=2.1×1011 Pa;泊松比均为μ=0.3;热膨胀系数α1=1×10-5/℃,α2=2×10-5/℃,α3=1×10-5/℃。图 2中的1+,1-分别表示第1层材料的上下界面,其余的2+,3-等表达方式类似。整个叠层材料受到的温度荷载表达式为ΔT=-10x2+20x

图 2 叠层材料几何形状以及温度荷载 Figure 2 Geometry of laminated material and thermal load

由于整个材料形状以及所受温度载荷对称,仅分析左边一半结构。首先进行有限元法计算的收敛性分析:分别沿x轴离散为10、20、30、40等份,沿y轴离散为12等份,每一层均划分4等份。有限元法均采用八节点二次单元。在每一层交界面上取几个特殊位置的点,并将该点处的位移及应力分量列入表 1

表 1 有限元法的收敛性分析 Table 1 The convergence of finite element method

表 1中的1+、1-分别表示第1层材料的上、下界面,3+、3-分别表示第3层材料的上、下界面。由表 1可以看出,当采用10×12的单元划分时,特定点的位移以及应力值已经收敛。由于本算例中最左侧竖边仅受竖向的位移约束,横向自由,其应变εyyεxx为零,且本算例的温度载荷在此处为零,代入式(13)可以解得左侧竖边上各点的应力σyy均为零。对于有限元法,可以发现,随着单元逐渐加密,由于左侧竖边的应变εxx逐渐逼近真实值零,使得应力σyy虽不准确等于边界条件,却逐渐变小,趋近于边界条件,说明单元加密能够提高计算精度。

然后进行状态空间法与有限元法的对比分析。对于本文方法,仅在上边划分20个线性单元。同时采用有限元软件ANSYS对结构进行分析,为提高计算精度,采用40×12的单元划分方式,在叠层材料左侧一半的矩形域共划分480个8节点实体二次单元,如图 3所示,其中SSM表示本文的状态空间法结果,FEM表示有限元法计算结果。

图 3 有限元单元网格划分 Figure 3 Mesh of laminated material for ANSYS

图 4图 5分别绘出了1、2层材料交界面处的位移uv。从中可以看出1~2层界面上位移的本文计算结果与有限元法计算结果吻合得很好,且均准确满足给定的位移边界条件。

图 4 1~2层界面横向位移 Figure 4 Horizontal interfacial displacement between the first and second layer
图 5 1~2层界面竖向位移 Figure 5 Vertical interfacial displacement between the first and second layer

图 6图 7分别给出了1、2层材料交界面处的应力场σxyσyy,其中σxy的本文解与有限元解能很好地吻合,σyy的本文解与有限元解在边界上略微有所偏差。由于将式(13)处理后得到的σyy值直接作为应力边界条件代入本文算法,因此本文在边界处的σyy严格等于边界条件。有限元法的计算过程为,首先利用变分原理求出各节点的位移量,再由位移量求解各节点的应力分量。导致有限元法在边界上只能准确满足位移边界条件,而应力边界条件只能近似满足。因此,本文方法与有限元方法结果在边界处略有差异,如图 7所示。且本文方法根据混合变分原理可以同时求出位移以及应力分量,减少了计算误差,在交界面上得到的应力值更加准确。

图 6 1~2层界面应力σxy Figure 6 Interfacial stress σxy between the first and second layer
图 7 1~2层界面应力σyy Figure 7 Interfacial stress σyy between the first and second layer

由于材料属性的不同,导致不同的材料有不同的应力值σxx图 8图 9绘出了1、2层材料交界面对应不同材料的应力σxx。由图 8图 9可以发现,本文解与有限元解吻合较好,且有限元解在x=0处近似满足边界条件σxx=0。

图 8 对应第1层材料的界面应力σxx Figure 8 Stress σxx of interface for the first layer
图 9 对应第2层材料的界面应力σxx Figure 9 Stress σxx of interface for the second layer
4 结语

在半解析状态空间有限元法的基础上,引入温度应力问题的物理方程,得出用于计算热应力问题的状态空间方程。计算发现,本文方法相对于传统的有限元法具有划分单元少的优点,且采用少量线性单元即可达到传统有限元法用二次单元加密划分所得结果的精度,节省了计算成本。本文方法边界处σyy为精确的应力边界条件,且本文方法采用迭代算法计算,单元数量与叠层材料的层数没有关系,使得计算过程更加简洁明了。在以后的研究中,可以在时域内应用差分格式,用状态空间法研究叠层材料瞬态温度应力问题。

参考文献
[1] 范家让. 强厚度叠层板壳的精确解[M]. 北京: 科学出版社, 1996 .
FAN Jiarang. Strong laminated plate and shell thickness of exact solutions[M]. Beijing: Science Press, 1996 . (in Chinese)
[2] Sheng H Y, Ye J Q. A state space finite element for laminated composite plates[J]. Computer Methods in Applied Mechanics and Engineering , 2002, 191 (37/38) : 4259–4276.
[3] Ye J Q, Sheng H Y, Qin Q H. A state space finite element for laminated composites with free edges and subjected to transverse and in-plane loads[J]. Computers & Structures , 2004, 82 (15/16) : 1131–1141.
[4] Wu C P, Jiang R Y. A state space differential reproducing kernel method for the 3D analysis of FGM sandwich circular hollow cylinders with combinations of simply-supported and clamped edges[J]. Composite Structures , 2012, 94 (11) : 3401–3420. DOI:10.1016/j.compstruct.2012.05.005
[5] Benedetti I, Aliabadi M H, Milazzo A. A fast BEM for the analysis of damaged structures with bonded piezoelectric sensors[J]. Computer Methods in Applied Mechanics and Engineering , 2010, 199 (9/10/11/12) : 490–501.
[6] Luo C, Moghtaderi B, Page A. Modelling of wall heat transfer using modified conduction transfer function, finite volume and complex Fourier analysis methods[J]. Energy and Buildings , 2010, 42 (5) : 605–617. DOI:10.1016/j.enbuild.2009.10.031
[7] 盛宏玉, 李和平, 叶建乔, 等. 层状介质热传导瞬态分析的一种新半数值解法[J]. 合肥工业大学学报(自然科学版) , 2010, 33 (5) : 709–712.
SHENG Hongyu, LI Heping, YE Jianqiao, et al. A new semi-numerical approach for the transient heat conduction analysis of laminated medium[J]. Journal of Hefei University of Technology(Natural Science) , 2010, 33 (5) : 709–712. (in Chinese)
[8] 盛宏玉, 李和平, 叶建乔, 等. 组合实心圆柱体热传导瞬态分析的状态变量法[J]. 上海交通大学学报 , 2011, 45 (2) : 247–251.
SHENG Hongyu, LI Heping, YE Jianqiao, et al. State variable method for transient heat conduction analysis of composite sold cylinder[J]. Journal of Shanghai Jiaotong University , 2011, 45 (2) : 247–251. (in Chinese)
[9] Maestre I R, Cubillas P R, Pérez-Lombard L. Transient heat conduction in multi-layer walls: An efficient strategy for Laplace's method[J]. Energy and Buildings , 2010, 42 (4) : 541–546. DOI:10.1016/j.enbuild.2009.10.023
[10] Li X Q, Chen Y, Spitler J D, et al. Applicability of calculation methods for conduction transfer function of building constructions[J]. International Journal of Thermal Sciences , 2009, 48 (7) : 1441–1451. DOI:10.1016/j.ijthermalsci.2008.11.006
[11] 田晓耕, 沈亚鹏. 广义热弹性问题研究进展[J]. 力学进展 , 2012, 42 (1) : 18–28.
TIAN Xiaogeng, SHEN Yapeng. Research progress in generalized thermoelastic problems[J]. Advances in Mechanics , 2012, 42 (1) : 18–28. (in Chinese)
[12] 李家宇, 程秀全, 卿光辉. 四边简支压电热弹性层合正交双曲壳的精确解[J]. 工程力学 , 2010, 27 (8) : 83–89.
LI Jiayu, CHENG Xiuquan, QING Guanghui. Exact solution for laminated piezothermoelastic generic shells with four simply supported edges[J]. Engineering Mechanics , 2010, 27 (8) : 83–89. (in Chinese)
[13] Youssef H M, El-Bary A A. Two-temperature generalized thermoelasticity with variable thermal conductivity[J]. Journal of Thermal Stresses , 2010, 33 (3) : 187–201. DOI:10.1080/01495730903454793
[14] Ezzat M A, Youssef H M. State space approach for conducting magneto-thermoelastic medium with variable electrical and thermal conductivity subjected to ramp-type heating[J]. Journal of Thermal Stresses , 2009, 32 (4) : 414–427. DOI:10.1080/01495730802637233
[15] Necas J, Hlavácek I. Mathematical theory of elastic and elasto-plastic bodies: an introduction[M]. New York: Elsevier Scientific Publishing Co. Amsterdam, 1981 .