第27卷 第4期 岩 土 工 程 学 报 Vol.27 No.4 2005年 4月 Chinese Journal of Geotechnical Engineering Apr.,2005
利用刚度矩阵法求解多层弹性半空间体的温度应力
Thermo-stress in multi-layered elastic half space solved with stiffness matrix
method
钟 阳，赵晓雷
（1.大连理工大学 土木学院，辽宁 大连116024；2.哈尔滨工业大学，黑龙江 哈尔滨 150006）
1
2
摘 要：从热弹性力学的基本方程出发，采用Hankel 积分变换和Laplace 积分变换等数学手段，首先推导出了单层弹性半空间轴对称体的温度应力问题的刚度矩阵，然后按传统有限元的方法组成总体刚度矩阵。通过求解由总体刚度矩阵所构成的代数方程组，再对其进行Hankel 和Laplace 积分逆变换就可解出在外荷载和温度联合作用下多层弹性半空间轴对称问题的解析解。由于刚度矩阵的元素中不含有正指数项，计算时不会出现溢出的现象，从而克服了传递矩阵法的缺点。在推导过程中，因不用事先人为的选择应力函数，使得问题的求解更加合理，同时也为进一步研究这类问题如湿度场、动力学等奠定了理论基础。最后，文中还给出了计算实例来证明推导结果的正确性。
关键词：多层弹性体半空间； 刚度矩阵； 积分变换；温度应力
中图分类号：TU 431 文献标识码：A 文章编号：1000–4548(2005)04–0374–04 作者简介：钟 阳(1955– )，男。大连理工大学，教授，博士生导师。
ZHONG Yang1, ZHAO Xiao-lei2
（1.Dep.of Civil. Engrg. Dalian Technology University,Dalian116024,China; 2.Harbin Institute of Technology, Harbin 150006, China）
Abstract: In the paper, thermo-stress in multilayered elastic half space is presented. Firstly, the stiffness matrix for a layer is derived based on the fundamental elasticity equations and some mathematic methods such as Hankel integral transformation. Then the global stiffness matrix is established for multilayered elastic half space using the finite element concepts in which layers are completely contacted. Therefore, explicit solution for axisymmetrical problems in multilayered elastic half space is obtained from the solution of the algebraic equation formed by global stiffness matrix and the inverse Hankel integral transformation. Because positive exponential function is not included in the element of matrix, the calculation is not overflowed. Therefore, the shortages of transfer matrix method are overcome. This method is clear in concept, and the corresponding formulas given in the paper are not only simple but also convenient for application. More important is that this method can be used to solve other problems of multilayered elastic half space such as thermo field and dynamics. An example of road surface deflection is presented to prove the calculated results.
Key words: multilayered elastic half space ; stiffness matrix; integral transformation；thermo-stress
0 引 言
在路面结构的计算及层状地基的分析中，多是以多层弹性半空间的理论解为依椐。虽然，国内外学者在这方面做了大量的研究工作，并且提出了许多求解方法，例如应力函数法[6]、传递矩阵法[1, 2]，但是，这些方法都不完善，主要缺点是在解的公式中含有正指数项，使得在计算时经常会出现溢出或病态矩阵现象，从而导致计算上的失败。刚度矩阵法[3, 4, 7]可以克服这一缺点，是一种较好的分析方法。
路面结构是裸露于大自然的，除了承受行驶车辆的静力荷载外，还要承受温度荷载。但是在现行的路面设计规范和计算方法中，还没有考虑温度荷载的影响。本文首先从热弹性力学的基本方程出发，利用 Hankel积分变换和Laplace 积分变换等数学手段，先推导出了单层弹性半空间轴对称体的温度应力问题的刚度矩
阵，然后按传统的有限元方法组成总体刚度矩阵。通过求解由总体刚度矩阵所构成的代数方程并对其进行Hankel 积分变换和Laplace 积分逆变换，可解出外荷载和温度联合作用下多层弹性半空间轴对称问题的解析解。由于在刚度矩阵的元素中不含有正指数项，所以计算时不会出现溢出的现象。从而克服了传递矩阵法的缺点。在本文的推导过程中，因不用事先人为的选择应力函数，使得问题的求解更加合理，同时也为进一步研究这类问题如湿度场、动力学等奠定了理论基础。最后，文中还给出了计算实例来证明推导结果的正确性。
1 刚度矩阵的推导
不计体力时，轴对称问题用位移表示的热应力平衡
───────
收稿日期:2004–05–19
第4期 钟 阳，等．利用刚度矩阵法求解多层弹性半空间体的温度应力 375
方程[1]为
1∂e u 2(1+µ) ∂T
+∇2u −=α ， (1)
1−2µ∂r r 1−2µ∂r 1∂e 2(1+µ) ∂T
+∇2w =α 。 (2)
1−2µ∂z 1−2µ∂z
物理方程式可表示成如下形式
ˆ(r , z , s ) e st d s 。 f (r , z , t ) =∫f
∞
又根据Hankel 变换的定义，对函数f (r , z , s ) 的n 阶
Hankel 正变换和逆变换分别为
~
f (ξ, z , s ) =
∞
∞
σr =2G ⎢
⎡µ∂u ⎤αET
， (3) e +⎥−
∂r ⎦1−2µ⎣1−2µ
(3)
−∞
∫rf (r , z , s ) J
n
(r ξ) d r ;
~
f (r , z , s ) =∫ξf (ξ, z , s ) J n (r ξ) d ξ 。
−∞
⎡µu ⎤αET
， (4) e +⎥− σθ=2G ⎢
r ⎦1−2µ⎣1−2µ
σz =2G ⎢
⎛∂u ∂w ⎞~ˆ~~ˆ 2 ˆ2(1+µ) ∂T + (6) d 2 ⎟ 。1d e w ~ˆ=0 ，(11) −ξw +−⎝∂z ∂r ⎠
d z 21−2µd z 1−2µ∂z
热传导方程为
~ˆd 2e ~ˆ2~∂T ˆ−α(1+µ) s T 2=0 ， (12) ξ −e λ∇T = ， (7) 2d z (1−µ) λ∂t
ˆ2~∂u u ∂w ∂21∂∂22T d ˆ2~式中 e =++++2；；∇= −q T =0 ， (13) 22∂r r ∂z r ∂r ∂r ∂r d z
E s 22G =；E 与G 分别为弹性模量和剪切弹性其中 q =ξ+ 。显然，经过Laplace 变换 和 2(1+µ) λ模量；µ与α分别为材料的泊松比和热膨胀系数；T Hankel 变换后，控制方程（9）转化为由式（10）~（13）
τzr =G ⎜
和λ分别为温度与热传导系数。式（1）、（2）和（7）构成弹性半空间轴对称体温度应力问题的控制方程。为了求解方便，对式（1）施加微分算子（2）对z 求偏导数并将两式相加可得到
⎡µ∂w ⎤αET 2~ˆξ ~d ˆ ， (5) u 2 ~ 2 ( 1 + µ ) ~−e +ˆ − ˆ +⎥ ξ=0 ， (10) −ξu e T 12µ1−2µ−∂z 2⎦⎣d z 1−2µ1−2µ
其中 J n (r ξ) 是n 阶贝塞尔函数。再对式（9）的第
一式施加一阶Hankel 变换，对式（9）的其它各式施加零阶Hankel 变换可得到
∂1
+，式∂r r
组成的常微分方程。其解为
~ˆqz −qz
T =m A 1e +B 1e ， (14)
~ˆ=s A e qz +B e −qz +4(1−2µ) A e ξz +B e −ξz ，e 1122
()
(λ((
)()
∇2e =α
1+µ2
∇T 。 (8) 1−µ
这样，弹性半空间轴对称体的温度应力问题的控制方程可以表示为如下一组方程：
⎧1∂e 2u −u =2(1+µ)α∂T ，+∇⎪1−2µ∂r ∂r r 1−2µ⎪
⎪1∂e +∇2w =2(1+µ)α∂T ，⎪∂z 1−2µ ⎨1−2µ∂z (9)
+µ1⎪∇2e =∇2T ，
⎪1−µ⎪2∂T
。⎪λ∇T =∂t ⎩
(15)
~ˆ=−ξA e qz +B e −qz +d A e ξz −d B e −ξz +A e ξz +B e −ξz ， u 11122233
(16)
~ˆ=q A e qz −B e −qz −d A e ξz −d B e −ξz +A e ξz +B e −ξz ， w 11122244
(17)
)
)
式中 m =α
2z ξ+12z ξ−1(1+µ) λ
。 ; d 1=; d 2=ξ(1−µ) s ξ
在式（14）~（17）中有8个积分常数A i 和B i (i =1, 2, 3, 4) ，但实际只有4个是独立的。把式（16）
首先对方程（9）实施Laplace 变换，其正变
换和逆变换分别记为
f ˆ(r , z , s ) =
∞
∫
f (r , z , t ) e −st d t ；
~ˆ
~~ˆ=ξu ˆ+∂w 可和（17）代入体积应变的积分变换式 e
∂z
以得到
ξA 4=−ξA 3+2⋅(3−4µ) A 2； ξB 4=ξB 3−2⋅(3−4µ) B 2 。
将A 4和B 4代入式（17）得 ~ˆ=q A e qz −B e −qz +d A e ξz −d B e −ξz −A e ξz +B e −ξz ，w 11324233
(18)
式中 ξd 3=7−8µ−2z ξ；ξd 4=7−8µ+2z ξ
()
376 岩 土 工 程 学 报 2005年
再对式（6）的第一式施加一阶Hankel 变换，对式（5）施加零阶Hankel 变换，并将式（16）、（17）和（18）代入可得到
qz −qz ξz −ξz ~ˆ=−2qG τξA −d 5A 2e ξz +d 6B 2e −ξz +ξA zr 1e −B 1e 3e −B 3e
， (19)
Q =λ
()()
~ˆ=2G ξ2(A e qz +B e −qz )+d A e ξz +d B e −ξz −ξ(A e ξz +B e −ξz )σz 11728233
， (20)
其中 d 5=a 1−2z ξ；d 6=a 1+2z ξ；
d 7=a 2−2z ξ；d 8=a 2+2z ξ；a 1=3−4µ；a 2=5−4µ 。
根据热力学中的热量的定义有
∂T
。 （21） ∂z
如图1所示，对于弹性层状半空间轴对称体中的任意一层的解可表示为
{U }=[A ]{C }；{σ}=[B ]{C } 。 （22）
式中
T
~~ˆˆ⎡~⎤~~~ˆˆˆˆ {U }=⎢u (ξ, h , s ) w (ξ, h , s ) u (ξ, 0, s ) w (ξ, 0, s ) T (ξ, h , s ) T (ξ, 0, s ) ⎥，
⎣⎦
T
~~ˆˆ⎡⎤~~~~ˆˆˆˆ {σ}=⎢τzr (ξ, h , s ) σz (ξ, h , s ) τzr (ξ, 0, s ) σz (ξ, 0, s ) Q (ξ, h , s ) /λQ (ξ, 0, s ) /λ⎥，
⎣⎦
{C }=[A 1
B 1A 2B 2A 3
B 3] ，
T
图 1 单一层的应力与位移及温度与热量关系的示意图
Fig1. Relationship of stress and strain with temperature and heat conduction for a layer
⎡−ξe qh −ξe −qh −d 1e ξh −d 2e −ξh e ξh e −ξh ⎤
⎢⎥
ξξqh −qh h −ξh h −ξh
⎢qe −qe −d 4e −e e ⎥d 3e ⎢⎥
⎢−ξ −ξ−−11⎥ ，⎢⎥[A ]=2G
⎢q −q −a −11⎥a ⎢⎥⎢me qh /2G me −qh /2G 0000⎥⎢⎥⎢m /2G m /2G 0000⎥⎣⎦⎡−ξqe qh ξqe −qh d 5e ξh d 6e −ξh ξe ξh −ξe −ξh ⎤
⎥⎢
ξh ξh −ξh −ξh 2−qh
⎢ξ2e qh d 8e ξe −d 7e −ξe −ξe ⎥
⎥⎢
⎢−ξq a 1ξq ξ−a 1−ξ⎥ ，
⎥⎢[B ]=2
⎢ξa 2a 2ξ2−ξ−ξ⎥
⎥⎢
⎢mqe qh −mqe −qh
0000⎥
⎥⎢
⎢mq 0000⎥−mq ⎦⎣
T
~~ˆˆ⎡⎤~~~~ˆˆˆˆ{σ}=⎢τzr (ξ, h , s ) σz (ξ, h , s ) τzr (ξ, 0, s ) σz (ξ, 0, s ) Q (ξ, h , s ) /λQ (ξ, 0, s ) /λ⎥,
⎣⎦
T
~~ˆˆ~~~~ˆ(ξ, h , s ) w ˆ(ξ, h , s ) −u ˆ(ξ, 0, s ) −w ˆ(ξ, 0, s ) T u {U }e =⎡(ξ, h , s ) −T (ξ, 0, s ) ⎤。
⎢⎥⎣⎦刚度矩阵[K ]e 中各元素k ij 的计算，可利用e
Mathematica [10] 完成，具体表达式见附录。
2 刚度矩阵在多层弹性体系中的应用
通过上面的推导，可建立单一层的状态方程为
~ˆ(ξ, h , s ) ⎫⎧τ⎡k 11k 12−k 13−k 14k 15
⎪⎪~⎢k ˆ
⎪σ(ξ, h , s ) ⎪⎢12k 22k 14−k 24k 25~ˆ(ξ, 0, s ) ⎪⎪⎪2G ⎢−k 13k 14k 11−k 12−k 16⎪τ⎨~ˆ(ξ, 0, s ) ⎬=∆⎢−k −k −k k 22k 26σ2412⎢14⎪⎪~ˆ⎢0⎪Q 000k 55(ξ, h , s ) ⎪
⎢⎪⎪~ˆ000−k 56⎢⎪⎣0⎭⎩Q (ξ, 0, s ) ⎪
e
−k 16⎤
−k 26⎥⎥k 15⎥⎥−k 25⎥−k 56⎥
⎥k 55⎥⎦
e
~ˆ(ξ, h , s ) ⎫⎧u
⎪⎪~ˆ
⎪w (ξ, h , s ) ⎪~ˆ(ξ, 0, s ) ⎪⎪⎪⎪−u ⎨~ˆ(ξ, 0, s ) ⎬−w ⎪⎪~ˆ⎪T (ξ, h , s ) ⎪
⎪⎪~ˆ⎪⎭⎩−T (ξ, 0, s ) ⎪
e
其中 a ξ=7−8µ 由式（21）和（22）可得到 {σ}e =[K ]e {U }e 。 (23) 上角标e 表示单元，[K ]表示弹性层状半空间轴对称
体中的任意一层的刚度矩阵。其中
k 12−k 13−k 14k 15−k 16⎤⎡k 11
⎢k ⎥−−k k k k k [1**********]6⎢⎥⎢−k k k 11−k 12−k 16k 15⎥
[K ]e =2G ⎢1314⎥
k 26−k 25⎥∆⎢−k 14−k 24−k 12k 22
⎢0k 55−k 56⎥000⎢⎥
−k k 0000⎢5655⎥⎣⎦
e
，(24)
对于层间完全接触的条件，以及其Laplace 变换和Hankel 变换可分别表示为
~~ˆ(ξ, z , s ) ⎤⎡u ˆ(ξ, z , s ) ⎤⎡u (r , z i , t ) ⎤⎡u (r , z i , t ) ⎤⎡u i i
⎢⎥⎢w (r , z , t ) ⎥⎢w (r , z , t ) ⎥⎢~~ˆ(ξ, z , s ) ˆ(ξ, z , s ) ⎥w w i i ⎥⎥⎢i i ⎥⎢⎥⎢⎢
~~ˆ(ξ, z , s ) ⎥⎢σˆ(ξ, z , s ) ⎥⎢σz (r , z i , t , ) ⎥⎢σz (r , z i , t ) ⎥⎢σz i z i
⎥ ⎢⎥⎢τ(r , z , t ) ⎥=⎢τ(r , z , t ) ⎥⎢τ~~=ˆˆ⎢τzr (ξ, z i , s ) ⎥zr (ξ, z i , s ) ⎥⎢zr i ⎥⎢zr i ⎥⎢
⎥⎥⎢~~ˆˆ⎢T (r , z i , t ) ⎥⎢T (r , z i , t ) ⎥⎢⎢T (ξ, z i , s ) ⎥⎢T (ξ, z i , s ) ⎥
⎢Q (r , z , t ) ⎥⎢Q (r , z , t ) ⎥⎢⎥⎥⎢~i i ⎦i +1~⎦i ⎣⎣ˆˆ
⎢⎦i +1⎣Q (ξ, z i , s ) ⎥⎦i ⎢⎣Q (ξ, z i , s ) ⎥ ，（25）
第4期 钟 阳，等．利用刚度矩阵法求解多层弹性半空间体的温度应力 377
其中下脚标i 表示第i 层。
利用有限元的概念和方法，根据式（24）和（25）可以建立多层体系的总刚度矩阵为
⎢ ⎢⎢⎢⎢⎣
1
2
N −1
~(ξ, h , s ) ⎫~ˆ(ξ, h , s ) ⎫⎧τ⎧u rz
⎤⎪~⎪⎪⎪~ˆ(ξ, h , s ) ⎥⎪w ⎪⎪σz (ξ, h , s ) ⎪⎥⎪⎪⎪⎪⋅. ⎪⎪⎥⎨⎬=⎨⋅⎬ 。⎥⎪. ⎪⎪~⎪ˆ~ˆ⎥⎪T (ξ, h h , s ) ⎪⎪Q (ξ, h n , s ) ⎪N ⎥⎪⎪~⎪⎦⎪~ˆˆQ (ξ, 0, s ) ⎪⎪n ⎭⎩T (ξ, 0, s ) ⎭⎩
图 4 路面温度应力的计算结果
Fig. 4 Calculated results of thermal stress of pavement
(26)
求解方程式（26）可求出任一深度处的状态向量，对状态向量进行Laplace 逆变换和Hankel 逆变换就得出真实问题的解析解。逆变换的具体算法见文献[8]。
4 结 论
(1)利用刚度矩阵法推导出了多层弹性半空间轴对称问题在层间完全接触情况下温度应力的解析解。
(2)由于刚度矩阵的元素中只含有负指数项，计算时不会出现溢出或病态矩阵，所以利用这种方法的计算非常稳定。
(3)在路面上表面层温度应力出现最大值。
(4)温度应力的大小和方向随时间和温度而变化，如图4所示，在不同的温度下路面的温度应力会出现正应力和负应力。这正是路面出现低温裂纹的原因。 参考文献：
[1] 钟 阳,王哲人,等.求解多层弹性半空间轴对称问题的传
3 计算实例
计算示例采用沥青路面结构的三层体系，见图2。
图2 沥青路面计算示意图
Fig.2 The calculation model of a practical flexible pavement
递矩阵法[J].土木工程学报,1992, 25(6): 37–43. [2] 钟 阳,王哲人.求解多层弹性半空间非轴对称问题的传递
矩阵法[J].土木工程学报,1995, 28 (1): 66–72.
[3] 钟 阳，陈静云,等.求解动荷载作用下多层粘弹性半空间
轴对称问题的精确刚度矩阵法[J].计算力学学报，2003，20(6): 749–755.
[4] 钟 阳，张永山.求解多层弹性半空间轴对称问题的精确刚
度矩阵法[J].力学季刊,2003, 24(3): 395–400.
[5] 王 凯.层状弹性连续体系在圆形均布垂直荷载作用下的
力学计算[J].土木工程学报,1982, 15 (2): 89–95. [6] 朱照宏,王秉纲,郭大智.路面力学计算[M].北京:人民交通
出版社,1985.
[7] Pestel E C, Lekie F A. Matrix Methods in Elasto-Mechanics
[M].1963.
[8] 陈光敬.成层横观各向同性地基的解析法分析及其工程应
用[D].上海:同济大学,1998.
[9] AL-Khoury R. Poroelastic spectral element for wave
propagation and parameter identification in multi-layer system[J]. International Journal of Solid and Structure, 2002, (39): 4073-4071.
[10] 嘉禾工作室.Mathematica 应用实例教程[M].北京:机械工
业出版社,2002.
图 3 路面温度的计算结果
Fig.3 Calculated results of pavement temperature
材料参数为：E 3=2500 MPa，h 3=0.15 m，α3=0.0022，
µ3=0.25，E 2=1500 MPa，h 2=0.25 m，α2=0.0028，λ3=1.0，
λ2=1.2，µ2=0.25，E 1=50 MPa，α1=0.0030;λ1=1.0，µ1=0.35。分别计算出路表面以下5 ，10，20 cm处的温度和温度应力随时间的变化，如图3和图4所示。由计算结果知，在路面上表面层温度应力出现最大值。温度应力的大小和方向随时间和温度而变化，见图4，在不同温度下路面的温度应力会出现正应力和负应力。这正是路面出现低温裂纹的原因。
第27卷 第4期 岩 土 工 程 学 报 Vol.27 No.4 2005年 4月 Chinese Journal of Geotechnical Engineering Apr.,2005
利用刚度矩阵法求解多层弹性半空间体的温度应力
Thermo-stress in multi-layered elastic half space solved with stiffness matrix
method
钟 阳，赵晓雷
（1.大连理工大学 土木学院，辽宁 大连116024；2.哈尔滨工业大学，黑龙江 哈尔滨 150006）
1
2
摘 要：从热弹性力学的基本方程出发，采用Hankel 积分变换和Laplace 积分变换等数学手段，首先推导出了单层弹性半空间轴对称体的温度应力问题的刚度矩阵，然后按传统有限元的方法组成总体刚度矩阵。通过求解由总体刚度矩阵所构成的代数方程组，再对其进行Hankel 和Laplace 积分逆变换就可解出在外荷载和温度联合作用下多层弹性半空间轴对称问题的解析解。由于刚度矩阵的元素中不含有正指数项，计算时不会出现溢出的现象，从而克服了传递矩阵法的缺点。在推导过程中，因不用事先人为的选择应力函数，使得问题的求解更加合理，同时也为进一步研究这类问题如湿度场、动力学等奠定了理论基础。最后，文中还给出了计算实例来证明推导结果的正确性。
关键词：多层弹性体半空间； 刚度矩阵； 积分变换；温度应力
中图分类号：TU 431 文献标识码：A 文章编号：1000–4548(2005)04–0374–04 作者简介：钟 阳(1955– )，男。大连理工大学，教授，博士生导师。
ZHONG Yang1, ZHAO Xiao-lei2
（1.Dep.of Civil. Engrg. Dalian Technology University,Dalian116024,China; 2.Harbin Institute of Technology, Harbin 150006, China）
Abstract: In the paper, thermo-stress in multilayered elastic half space is presented. Firstly, the stiffness matrix for a layer is derived based on the fundamental elasticity equations and some mathematic methods such as Hankel integral transformation. Then the global stiffness matrix is established for multilayered elastic half space using the finite element concepts in which layers are completely contacted. Therefore, explicit solution for axisymmetrical problems in multilayered elastic half space is obtained from the solution of the algebraic equation formed by global stiffness matrix and the inverse Hankel integral transformation. Because positive exponential function is not included in the element of matrix, the calculation is not overflowed. Therefore, the shortages of transfer matrix method are overcome. This method is clear in concept, and the corresponding formulas given in the paper are not only simple but also convenient for application. More important is that this method can be used to solve other problems of multilayered elastic half space such as thermo field and dynamics. An example of road surface deflection is presented to prove the calculated results.
Key words: multilayered elastic half space ; stiffness matrix; integral transformation；thermo-stress
0 引 言
在路面结构的计算及层状地基的分析中，多是以多层弹性半空间的理论解为依椐。虽然，国内外学者在这方面做了大量的研究工作，并且提出了许多求解方法，例如应力函数法[6]、传递矩阵法[1, 2]，但是，这些方法都不完善，主要缺点是在解的公式中含有正指数项，使得在计算时经常会出现溢出或病态矩阵现象，从而导致计算上的失败。刚度矩阵法[3, 4, 7]可以克服这一缺点，是一种较好的分析方法。
路面结构是裸露于大自然的，除了承受行驶车辆的静力荷载外，还要承受温度荷载。但是在现行的路面设计规范和计算方法中，还没有考虑温度荷载的影响。本文首先从热弹性力学的基本方程出发，利用 Hankel积分变换和Laplace 积分变换等数学手段，先推导出了单层弹性半空间轴对称体的温度应力问题的刚度矩
阵，然后按传统的有限元方法组成总体刚度矩阵。通过求解由总体刚度矩阵所构成的代数方程并对其进行Hankel 积分变换和Laplace 积分逆变换，可解出外荷载和温度联合作用下多层弹性半空间轴对称问题的解析解。由于在刚度矩阵的元素中不含有正指数项，所以计算时不会出现溢出的现象。从而克服了传递矩阵法的缺点。在本文的推导过程中，因不用事先人为的选择应力函数，使得问题的求解更加合理，同时也为进一步研究这类问题如湿度场、动力学等奠定了理论基础。最后，文中还给出了计算实例来证明推导结果的正确性。
1 刚度矩阵的推导
不计体力时，轴对称问题用位移表示的热应力平衡
───────
收稿日期:2004–05–19
第4期 钟 阳，等．利用刚度矩阵法求解多层弹性半空间体的温度应力 375
方程[1]为
1∂e u 2(1+µ) ∂T
+∇2u −=α ， (1)
1−2µ∂r r 1−2µ∂r 1∂e 2(1+µ) ∂T
+∇2w =α 。 (2)
1−2µ∂z 1−2µ∂z
物理方程式可表示成如下形式
ˆ(r , z , s ) e st d s 。 f (r , z , t ) =∫f
∞
又根据Hankel 变换的定义，对函数f (r , z , s ) 的n 阶
Hankel 正变换和逆变换分别为
~
f (ξ, z , s ) =
∞
∞
σr =2G ⎢
⎡µ∂u ⎤αET
， (3) e +⎥−
∂r ⎦1−2µ⎣1−2µ
(3)
−∞
∫rf (r , z , s ) J
n
(r ξ) d r ;
~
f (r , z , s ) =∫ξf (ξ, z , s ) J n (r ξ) d ξ 。
−∞
⎡µu ⎤αET
， (4) e +⎥− σθ=2G ⎢
r ⎦1−2µ⎣1−2µ
σz =2G ⎢
⎛∂u ∂w ⎞~ˆ~~ˆ 2 ˆ2(1+µ) ∂T + (6) d 2 ⎟ 。1d e w ~ˆ=0 ，(11) −ξw +−⎝∂z ∂r ⎠
d z 21−2µd z 1−2µ∂z
热传导方程为
~ˆd 2e ~ˆ2~∂T ˆ−α(1+µ) s T 2=0 ， (12) ξ −e λ∇T = ， (7) 2d z (1−µ) λ∂t
ˆ2~∂u u ∂w ∂21∂∂22T d ˆ2~式中 e =++++2；；∇= −q T =0 ， (13) 22∂r r ∂z r ∂r ∂r ∂r d z
E s 22G =；E 与G 分别为弹性模量和剪切弹性其中 q =ξ+ 。显然，经过Laplace 变换 和 2(1+µ) λ模量；µ与α分别为材料的泊松比和热膨胀系数；T Hankel 变换后，控制方程（9）转化为由式（10）~（13）
τzr =G ⎜
和λ分别为温度与热传导系数。式（1）、（2）和（7）构成弹性半空间轴对称体温度应力问题的控制方程。为了求解方便，对式（1）施加微分算子（2）对z 求偏导数并将两式相加可得到
⎡µ∂w ⎤αET 2~ˆξ ~d ˆ ， (5) u 2 ~ 2 ( 1 + µ ) ~−e +ˆ − ˆ +⎥ ξ=0 ， (10) −ξu e T 12µ1−2µ−∂z 2⎦⎣d z 1−2µ1−2µ
其中 J n (r ξ) 是n 阶贝塞尔函数。再对式（9）的第
一式施加一阶Hankel 变换，对式（9）的其它各式施加零阶Hankel 变换可得到
∂1
+，式∂r r
组成的常微分方程。其解为
~ˆqz −qz
T =m A 1e +B 1e ， (14)
~ˆ=s A e qz +B e −qz +4(1−2µ) A e ξz +B e −ξz ，e 1122
()
(λ((
)()
∇2e =α
1+µ2
∇T 。 (8) 1−µ
这样，弹性半空间轴对称体的温度应力问题的控制方程可以表示为如下一组方程：
⎧1∂e 2u −u =2(1+µ)α∂T ，+∇⎪1−2µ∂r ∂r r 1−2µ⎪
⎪1∂e +∇2w =2(1+µ)α∂T ，⎪∂z 1−2µ ⎨1−2µ∂z (9)
+µ1⎪∇2e =∇2T ，
⎪1−µ⎪2∂T
。⎪λ∇T =∂t ⎩
(15)
~ˆ=−ξA e qz +B e −qz +d A e ξz −d B e −ξz +A e ξz +B e −ξz ， u 11122233
(16)
~ˆ=q A e qz −B e −qz −d A e ξz −d B e −ξz +A e ξz +B e −ξz ， w 11122244
(17)
)
)
式中 m =α
2z ξ+12z ξ−1(1+µ) λ
。 ; d 1=; d 2=ξ(1−µ) s ξ
在式（14）~（17）中有8个积分常数A i 和B i (i =1, 2, 3, 4) ，但实际只有4个是独立的。把式（16）
首先对方程（9）实施Laplace 变换，其正变
换和逆变换分别记为
f ˆ(r , z , s ) =
∞
∫
f (r , z , t ) e −st d t ；
~ˆ
~~ˆ=ξu ˆ+∂w 可和（17）代入体积应变的积分变换式 e
∂z
以得到
ξA 4=−ξA 3+2⋅(3−4µ) A 2； ξB 4=ξB 3−2⋅(3−4µ) B 2 。
将A 4和B 4代入式（17）得 ~ˆ=q A e qz −B e −qz +d A e ξz −d B e −ξz −A e ξz +B e −ξz ，w 11324233
(18)
式中 ξd 3=7−8µ−2z ξ；ξd 4=7−8µ+2z ξ
()
376 岩 土 工 程 学 报 2005年
再对式（6）的第一式施加一阶Hankel 变换，对式（5）施加零阶Hankel 变换，并将式（16）、（17）和（18）代入可得到
qz −qz ξz −ξz ~ˆ=−2qG τξA −d 5A 2e ξz +d 6B 2e −ξz +ξA zr 1e −B 1e 3e −B 3e
， (19)
Q =λ
()()
~ˆ=2G ξ2(A e qz +B e −qz )+d A e ξz +d B e −ξz −ξ(A e ξz +B e −ξz )σz 11728233
， (20)
其中 d 5=a 1−2z ξ；d 6=a 1+2z ξ；
d 7=a 2−2z ξ；d 8=a 2+2z ξ；a 1=3−4µ；a 2=5−4µ 。
根据热力学中的热量的定义有
∂T
。 （21） ∂z
如图1所示，对于弹性层状半空间轴对称体中的任意一层的解可表示为
{U }=[A ]{C }；{σ}=[B ]{C } 。 （22）
式中
T
~~ˆˆ⎡~⎤~~~ˆˆˆˆ {U }=⎢u (ξ, h , s ) w (ξ, h , s ) u (ξ, 0, s ) w (ξ, 0, s ) T (ξ, h , s ) T (ξ, 0, s ) ⎥，
⎣⎦
T
~~ˆˆ⎡⎤~~~~ˆˆˆˆ {σ}=⎢τzr (ξ, h , s ) σz (ξ, h , s ) τzr (ξ, 0, s ) σz (ξ, 0, s ) Q (ξ, h , s ) /λQ (ξ, 0, s ) /λ⎥，
⎣⎦
{C }=[A 1
B 1A 2B 2A 3
B 3] ，
T
图 1 单一层的应力与位移及温度与热量关系的示意图
Fig1. Relationship of stress and strain with temperature and heat conduction for a layer
⎡−ξe qh −ξe −qh −d 1e ξh −d 2e −ξh e ξh e −ξh ⎤
⎢⎥
ξξqh −qh h −ξh h −ξh
⎢qe −qe −d 4e −e e ⎥d 3e ⎢⎥
⎢−ξ −ξ−−11⎥ ，⎢⎥[A ]=2G
⎢q −q −a −11⎥a ⎢⎥⎢me qh /2G me −qh /2G 0000⎥⎢⎥⎢m /2G m /2G 0000⎥⎣⎦⎡−ξqe qh ξqe −qh d 5e ξh d 6e −ξh ξe ξh −ξe −ξh ⎤
⎥⎢
ξh ξh −ξh −ξh 2−qh
⎢ξ2e qh d 8e ξe −d 7e −ξe −ξe ⎥
⎥⎢
⎢−ξq a 1ξq ξ−a 1−ξ⎥ ，
⎥⎢[B ]=2
⎢ξa 2a 2ξ2−ξ−ξ⎥
⎥⎢
⎢mqe qh −mqe −qh
0000⎥
⎥⎢
⎢mq 0000⎥−mq ⎦⎣
T
~~ˆˆ⎡⎤~~~~ˆˆˆˆ{σ}=⎢τzr (ξ, h , s ) σz (ξ, h , s ) τzr (ξ, 0, s ) σz (ξ, 0, s ) Q (ξ, h , s ) /λQ (ξ, 0, s ) /λ⎥,
⎣⎦
T
~~ˆˆ~~~~ˆ(ξ, h , s ) w ˆ(ξ, h , s ) −u ˆ(ξ, 0, s ) −w ˆ(ξ, 0, s ) T u {U }e =⎡(ξ, h , s ) −T (ξ, 0, s ) ⎤。
⎢⎥⎣⎦刚度矩阵[K ]e 中各元素k ij 的计算，可利用e
Mathematica [10] 完成，具体表达式见附录。
2 刚度矩阵在多层弹性体系中的应用
通过上面的推导，可建立单一层的状态方程为
~ˆ(ξ, h , s ) ⎫⎧τ⎡k 11k 12−k 13−k 14k 15
⎪⎪~⎢k ˆ
⎪σ(ξ, h , s ) ⎪⎢12k 22k 14−k 24k 25~ˆ(ξ, 0, s ) ⎪⎪⎪2G ⎢−k 13k 14k 11−k 12−k 16⎪τ⎨~ˆ(ξ, 0, s ) ⎬=∆⎢−k −k −k k 22k 26σ2412⎢14⎪⎪~ˆ⎢0⎪Q 000k 55(ξ, h , s ) ⎪
⎢⎪⎪~ˆ000−k 56⎢⎪⎣0⎭⎩Q (ξ, 0, s ) ⎪
e
−k 16⎤
−k 26⎥⎥k 15⎥⎥−k 25⎥−k 56⎥
⎥k 55⎥⎦
e
~ˆ(ξ, h , s ) ⎫⎧u
⎪⎪~ˆ
⎪w (ξ, h , s ) ⎪~ˆ(ξ, 0, s ) ⎪⎪⎪⎪−u ⎨~ˆ(ξ, 0, s ) ⎬−w ⎪⎪~ˆ⎪T (ξ, h , s ) ⎪
⎪⎪~ˆ⎪⎭⎩−T (ξ, 0, s ) ⎪
e
其中 a ξ=7−8µ 由式（21）和（22）可得到 {σ}e =[K ]e {U }e 。 (23) 上角标e 表示单元，[K ]表示弹性层状半空间轴对称
体中的任意一层的刚度矩阵。其中
k 12−k 13−k 14k 15−k 16⎤⎡k 11
⎢k ⎥−−k k k k k [1**********]6⎢⎥⎢−k k k 11−k 12−k 16k 15⎥
[K ]e =2G ⎢1314⎥
k 26−k 25⎥∆⎢−k 14−k 24−k 12k 22
⎢0k 55−k 56⎥000⎢⎥
−k k 0000⎢5655⎥⎣⎦
e
，(24)
对于层间完全接触的条件，以及其Laplace 变换和Hankel 变换可分别表示为
~~ˆ(ξ, z , s ) ⎤⎡u ˆ(ξ, z , s ) ⎤⎡u (r , z i , t ) ⎤⎡u (r , z i , t ) ⎤⎡u i i
⎢⎥⎢w (r , z , t ) ⎥⎢w (r , z , t ) ⎥⎢~~ˆ(ξ, z , s ) ˆ(ξ, z , s ) ⎥w w i i ⎥⎥⎢i i ⎥⎢⎥⎢⎢
~~ˆ(ξ, z , s ) ⎥⎢σˆ(ξ, z , s ) ⎥⎢σz (r , z i , t , ) ⎥⎢σz (r , z i , t ) ⎥⎢σz i z i
⎥ ⎢⎥⎢τ(r , z , t ) ⎥=⎢τ(r , z , t ) ⎥⎢τ~~=ˆˆ⎢τzr (ξ, z i , s ) ⎥zr (ξ, z i , s ) ⎥⎢zr i ⎥⎢zr i ⎥⎢
⎥⎥⎢~~ˆˆ⎢T (r , z i , t ) ⎥⎢T (r , z i , t ) ⎥⎢⎢T (ξ, z i , s ) ⎥⎢T (ξ, z i , s ) ⎥
⎢Q (r , z , t ) ⎥⎢Q (r , z , t ) ⎥⎢⎥⎥⎢~i i ⎦i +1~⎦i ⎣⎣ˆˆ
⎢⎦i +1⎣Q (ξ, z i , s ) ⎥⎦i ⎢⎣Q (ξ, z i , s ) ⎥ ，（25）
第4期 钟 阳，等．利用刚度矩阵法求解多层弹性半空间体的温度应力 377
其中下脚标i 表示第i 层。
利用有限元的概念和方法，根据式（24）和（25）可以建立多层体系的总刚度矩阵为
⎢ ⎢⎢⎢⎢⎣
1
2
N −1
~(ξ, h , s ) ⎫~ˆ(ξ, h , s ) ⎫⎧τ⎧u rz
⎤⎪~⎪⎪⎪~ˆ(ξ, h , s ) ⎥⎪w ⎪⎪σz (ξ, h , s ) ⎪⎥⎪⎪⎪⎪⋅. ⎪⎪⎥⎨⎬=⎨⋅⎬ 。⎥⎪. ⎪⎪~⎪ˆ~ˆ⎥⎪T (ξ, h h , s ) ⎪⎪Q (ξ, h n , s ) ⎪N ⎥⎪⎪~⎪⎦⎪~ˆˆQ (ξ, 0, s ) ⎪⎪n ⎭⎩T (ξ, 0, s ) ⎭⎩
图 4 路面温度应力的计算结果
Fig. 4 Calculated results of thermal stress of pavement
(26)
求解方程式（26）可求出任一深度处的状态向量，对状态向量进行Laplace 逆变换和Hankel 逆变换就得出真实问题的解析解。逆变换的具体算法见文献[8]。
4 结 论
(1)利用刚度矩阵法推导出了多层弹性半空间轴对称问题在层间完全接触情况下温度应力的解析解。
(2)由于刚度矩阵的元素中只含有负指数项，计算时不会出现溢出或病态矩阵，所以利用这种方法的计算非常稳定。
(3)在路面上表面层温度应力出现最大值。
(4)温度应力的大小和方向随时间和温度而变化，如图4所示，在不同的温度下路面的温度应力会出现正应力和负应力。这正是路面出现低温裂纹的原因。 参考文献：
[1] 钟 阳,王哲人,等.求解多层弹性半空间轴对称问题的传
3 计算实例
计算示例采用沥青路面结构的三层体系，见图2。
图2 沥青路面计算示意图
Fig.2 The calculation model of a practical flexible pavement
递矩阵法[J].土木工程学报,1992, 25(6): 37–43. [2] 钟 阳,王哲人.求解多层弹性半空间非轴对称问题的传递
矩阵法[J].土木工程学报,1995, 28 (1): 66–72.
[3] 钟 阳，陈静云,等.求解动荷载作用下多层粘弹性半空间
轴对称问题的精确刚度矩阵法[J].计算力学学报，2003，20(6): 749–755.
[4] 钟 阳，张永山.求解多层弹性半空间轴对称问题的精确刚
度矩阵法[J].力学季刊,2003, 24(3): 395–400.
[5] 王 凯.层状弹性连续体系在圆形均布垂直荷载作用下的
力学计算[J].土木工程学报,1982, 15 (2): 89–95. [6] 朱照宏,王秉纲,郭大智.路面力学计算[M].北京:人民交通
出版社,1985.
[7] Pestel E C, Lekie F A. Matrix Methods in Elasto-Mechanics
[M].1963.
[8] 陈光敬.成层横观各向同性地基的解析法分析及其工程应
用[D].上海:同济大学,1998.
[9] AL-Khoury R. Poroelastic spectral element for wave
propagation and parameter identification in multi-layer system[J]. International Journal of Solid and Structure, 2002, (39): 4073-4071.
[10] 嘉禾工作室.Mathematica 应用实例教程[M].北京:机械工
业出版社,2002.
图 3 路面温度的计算结果
Fig.3 Calculated results of pavement temperature
材料参数为：E 3=2500 MPa，h 3=0.15 m，α3=0.0022，
µ3=0.25，E 2=1500 MPa，h 2=0.25 m，α2=0.0028，λ3=1.0，
λ2=1.2，µ2=0.25，E 1=50 MPa，α1=0.0030;λ1=1.0，µ1=0.35。分别计算出路表面以下5 ，10，20 cm处的温度和温度应力随时间的变化，如图3和图4所示。由计算结果知，在路面上表面层温度应力出现最大值。温度应力的大小和方向随时间和温度而变化，见图4，在不同温度下路面的温度应力会出现正应力和负应力。这正是路面出现低温裂纹的原因。