第五章  沥青路面应力分析 一.古典设计方法 1.麻省公式  图5-1 古典公式示意图     1901年,美国麻省道路委员会第八次年会上发表了世界上第一个路面设计的公式。它假定汽车是一个集中荷载P,荷载以45(角通过碎石基层分布于边长为碎石层厚2倍的正方形面积的土基上,所以:  (5-1)  2.Downs公式     1933年,Downs对麻省公式进行修正,认为荷载在路面层内的传布与垂直方向成某一分布角(的圆锥上,所以传到路面的顶面时,压力分布于一个圆形的面积上而不是正方形,但他仍假定汽车荷载为集中荷载。据此:  图5-2 古典公式改进 (5-2)  3.Gray公式     1934年、Gray认为由于汽车荷载轮胎接触路面由一个面积,所以不应当假定汽车荷载为集中荷载,而应当假定汽车荷载为圆形均布荷载,并设轮载接地圆形面积的半径为a,即:  (5-3)  4.评述   古典理论公式是假定路面只要起分布荷载的作用,采用简单的分布角的概念,这个朴素思想的路面力学理论应予解决的问题; 从各公式得知,路面厚度主要取决于土基承载力得大小,这就是土基强度得问题。但初期没有提出土基参数的测定问题;   古典公式以轮载作为交通荷载,它不能反映交通量的因素,这在当时轻交通时代可能矛盾不突出,但随着交通得发展,不考虑交通量是无法使用的解决的办法就是在土基承载力取值上应根据交通量的大小采取不同的安全系数。 二.弹性半空间体     1.解答过程     1887~1885 布辛尼斯克得到完整的解答,方法是采用半逆解法。     1925年 A.E.Love势能法得到了解答。     采用路面力学中的方法,同样可以得到解答。     2.A.E.Love解     轮隙弯沉的计算及应用采用以上公式   三.多层体系    1.解答过程    1945年,D.M.Burmister得到理论解.    1945-1955 研究层状体系的工程应用    1955,R.L.希夫曼得到非轴对称的解    2.计算方法     采用查诺模图法     采用程序计算法 四.计算程序   沥青路面通常是多层体系。自从本世纪四十年代以来无论在理论分析,还是在数值计算方面,都取得很大进展,特别是计算机科学的发展及其在工程技术中广泛应用,使层状体系理论的研究的日趋完善,其中有波米斯特(D.M.Burmister)(1945年)及英因福克斯(L.Fox)、阿堪姆(W.E.Acum)、苏联科岗(Korah)及英国琼斯(A.Jones)等所作的贡献。在荷载形式方面,包括轴对称均布荷载与非轴对称单向水平荷载,都可直接进行数值计算,在层次结构方面,由双层体系、三层体系发展到多层体系。在计算机程序方面,有壳牌公司编制的Bisar程序,雪弗隆公司编制的Chevron程序,美国地沥青学会所采用的DAMA程序。   1.基本图式与基本假定 多层体系在圆形均布垂直荷载作用下的计算图式如图1所示。   层状体系基本假定:   (1)各层都是由均质、各向同性的弹性材料组成,这种材料的力学性能服从虎克定律;   (2)假定土基在水平方向和向下的深度方向均为无限,其上的各层厚度均为有限,但水平方向仍为无限;   (3)上层表示作用着轴对称圆形均布垂直荷载,同时在上层无限深度处及水平无限远处应力和应变都是零;   (4)层间接触面假定完全连续。  图1计算图式    2.基本原理   根据弹性理论,对于轴对称空间体,其几何方程为:    其物理方程为:       式中:G为剪数模量,      μ为弹性体的泊桑比   轴对称空间课题微分单元的平衡微分方程为:  从式1~式3看出,三式中共有十个变量,并且已有十个方程式,结合边界条件即可解出未知量值。但这种解法相当困难,甚至不可能得到应力分量。因此一般采用应力函数求解。研究物体的变形一般是针对物体内部割出的一块微分单元体,显然各相邻单元体的变形应是谐调的。所以物体在变形前是一个连续体,在变形后也应是一个连续体。消去位移分量,可得变形连续方程为:  Θ——第一应力不变量, 变形连续方程又称相容条件,是由圣维南(B.desaint-Venant)于1864年提出的。实际上式4应有四个相容条件,但确是等效的。采用应力函数法求解轴对称课题主要有Love函数法及Southwell函数法,这里介绍Love函数法。 设应力函数φ=φ(r,z)并给定     将式5代入平衡微分方程式3和变形连续方程式4,除平衡微分方程中第一个恒等于零外,其余全部转化为重调和方程,即    这就是说,如果应力函数φ是重调和方程的解,则能满足平衡微分方程的变形连续方程。并可由式5求得应力分量,再由物理方程求得应变分量。位移分量可由下式求得。    重调和方程的求解可采用分离变量法。对于多层体系中某一层j,可以给定应力函数为:  代入重调和方程可以得出:  令 ,解以上方程式,可得应力函数为:  式中:ξ——参数;    J0(ξρ)——第一类零阶贝塞尔函数;    ——无量纲系数;    Aj,Bj,Cj,Dj为积分常数,可由每一层的边界条件和层间结合条件等确定。 下标j从1到n,表示同该层次相应的计算参数。将应力函数式(10)代入洛夫应力函数与应力关系式。  以上表达的各项分量并非因荷载p(ρ)所引起,而是由σz=-ξJ0(ξr)所引起。通过Hankel变换等推演过程,可以求得:  (12) 假设:  (13) 假设式(11)的各项分量为Y*,而由实际荷载p(ρ)产生的各项应力、位移分量为Y,这两者有以下关系:  (14) 将式(11)代入:  由σz=q(ρ)这个实际荷载引起的各项分量。  3.积分常数计算   对n层体系具有4n个积分常数。由以上算式可以看出,多层体系的应力应变计算的关键是要确定对应于各个层次的积分常数,然后通过贝塞尔函数的无穷积分计算,便可完成全部计算分析工作。确定积分常数,可以根据相应的边界条件与层间结合条件来进行。在多层体系顶面(j=1,λ=0)具有以下边界条件: (σ(z)1=-ξJ0(ξρ) (τ(zr)1=0 (16)   在第j层与第j+1层之间的结合面上(λ=λj),若这两层是完全连续的,则具有以下连续条件:   (σ(z)j=(σ(z)j+1   (τ(zr)j=(τ(zr)j+1    (u()j=(u()j+1    (w()j=(w()j+1 (17)   此外,在地基的无限深处,应力与位移皆满足    (σz,σ(,τzr,u,w)r→∞=0  (18)   则要求应力函数φn|(((=0,因e((|(((=∞,即: An=Cn=0(19) 则对于n层体系,还有4n-2个待定积分系数,而根据边界条件可以建立4n-2个方程式,因此全部积分系数均可以求解。确定待定积分系数,用矩阵法非常简单,便于使用计算机分析计算。为此可将应力和位移中包含有Aj, Bj,Cj,Dj的系数写成矩阵形式:  式中hj=λj-λj-1 ; M——4×4的矩阵   根据连续条件,可以写成:  由式21可以看出,第j层积分常数可由第j+1层的积分常数求得。 通过逐层计算,可以将第一层的积分常数与第n层的积分常数联系起来。并利用下式可得:  由多层体系顶面的边界条件代入(11)得:   则:     因此,在计算积分常数时,可按以下步骤进行计算: 1)形成矩阵[C] 2)形成矩阵[F] 3)计算Bn,Dn 4)由下而上逐层计算各层的积分常数。   在积分常数确定之后,通过贝塞尔函数及无穷积分数值可计算应力分量及位移分量。   4.数值积分   在进行应力分量及位移分量计算时,可以归纳为以下的形式:  5.贝塞尔函数的计算 贝塞尔方程  (25) 贝塞尔函数的解  n=0  n=1  (1)函数特性     当x=0时,J0(x)=1,其它各阶均为零     是衰减函数,趣近于正弦函数   (2)贝塞尔函数的计算    当x小于4时  i ai bi i ai bi  0 1. 2. 4 0.4442584 0.1777583  1 -3.999999 -4. 5 -0.0709253 -0.0236617  2 3.9999973 2.6666661 6 0.0076772 0.0022069  3 -1.777756 -0.888884 7 -0.0005014 -0.0001290  当x大于4时   i  P0i   Q0i   P1i   Q1i  0 0.3989423 -0.0124669 0.3989423 0.0374008  1 -0.0017531 0.0004564 0.0029218 -0.000639  2 0.0001734 -0.000087 -0.0002232 0.0001065  3 -0.0000488 0.0000342 0.0000581 -0.0000399  4 0.0000174 -0.0000142 -0.0000201 0.0000162  5 -0.0000037 0.0000032 0.0000042 -0.0000037   6.贝塞尔函数的无穷积分   (1)高斯积分   (2)高斯系数 高斯型点P的选择,可以根据计算精度确定。 p t Ak  2 0.5773502692 1.  3 0.7745966692 0.55555555556   0. 0.88888888889  4 0.8611363116 0.3478548451   0.3399810436 0.65211451549  5 0.9061798459 0.2369268851   0.5384693101 0.4786286705   0. 0.5888888889  (3)贝塞尔函数的积分 为避免因贝塞尔函数的波动性而引起积分误差,在积分时采用贝塞尔函数零点积分的办法,即首先找到贝塞尔函数的零点,然后在该范围内进行高斯积分,则:   4.计算程序 五.层状体系计算程序   层状体系应力,应变及位移的计算时主要涉及的是贝塞尔函数的无穷积分,因此无穷积分的精度将直接影响计算精度。因为贝塞尔函数是波动衰减函数,如果采用第二章所述的高斯积分,必须合理选取高斯积分段。在一般的数值积分法中主要采用等值增长的办法,使计算相对误差达到规定的精度。但由于贝塞尔函数是波动衰减函数,如果高斯积分区段一端函数值为正,另一端函数值为负,高斯积分点或为正,或为负,那么计算结果误差则比较大,为了防止以上这种情况的发生,在程序中采用零点分段的办法。由于贝塞尔函数的零点为已知,那么零点与其它任何数相乘均为零,那么积分时选用的分段区间为相邻两零点之间,则积分结果精度较高。 6.求解N层体系积分常数Ai,Bi,Ci,Di子程序(SOL) 积分常数计算顺序由下而上进行,即由第n层的积分常数An,Bn,Cn,Dn计算An-1,Bn-1,Cn-1,Dn-1,然后逐层向上,直到希望计算的某一层。程序中,积分常数计算的主要任务是确定系数阵[F]及[C]。其执行程序为(SOL)。 程序中符号说明 TM——相当于ξhj; H(NH)——每层的结构厚度(取总厚度的相对值); Z(NH)——各层界面的竖向坐标(取总厚度的相对值); NH——结构层数(不包括第N层土基,NH=N-1); PR(N)——各层次泊桑比; EE(NH)——相当于 VV(4)——存放每个层次的积分常数Aj,Bj,Cj,Dj; LPT(NH)——层次结构的顺序; A(4,4)——系数矩阵[Nj]; C(4,4)——系数矩阵[C]; F(2,4)——系数矩阵[F]。 7.积分计算子程序(同前) 8.余项计算 在某一区段内无穷积分积的执行程序可见多层体系计算程序。在程序中余项值的计算由控制变量INTT控制执行。当INTT值为零时,表明计算点在层状体系顶面(Z=0),则要求计算余顶值,当INTT为1时,表明计算点在层状体系顶面以下点(Z≠0),则不要计算余项。余项计算值代表符VS11,RS11,TS11及W11分别表示σz,σr,τzr及ω在有限积分段的余项值。 9.多层体系应力、应变及位移计算程序说明 本程序适用于多层体系结构,对双圆或多圆均布荷载可采用应力迭加原理得到。执行程序如下: 程序中主要变量说明 INTT——计算点位置符(INTT=0,表示要求余项); IC——积分次数; VS(…),RS(…),TS(…),w(·)分别存放每次积分后σz,σr,(zr,w的对应值  DEL——积分的相对精度 WT(·),D(·)分别为高斯积分宽度和高斯积分点。 其它变量名同子程序SOL 六、双圆或多圆荷载应力的计算 现行的路面设计规范多用双圆均匀荷载,利用本程序进行应力计算时,必须对源程序进行修正。修正的基本方法是根据应力迭加原理。要求荷轮隙、轮印中心及轮印中心两侧的应力或位移,可利用计算两点位移或应力迭加的方法。 七.弹性多层体系应力、位移分析程序(AP01)    1.程序功能   弹性多层体系应力,位移分析程序适用于N层组成的多层结构体系,具有如下功能: 1)适用于多层弹性体系,层数不限,在此最大值定为L=6; 2)每个层次的弹性模量和泊桑比不受限制; 3)适用于计算各层体系任意一点的应力,位移计算,可同时算出多个点的应力及位移,计算点最大值为24点; 4)荷载为单园垂直均布荷载,作用于上层顶面; 5)对双圆荷载,则利用单圆荷载进行应力迭加。   2.程序的输入输出说明 输入变量 NL——层状体系的层数(NL≤6); NS——应力、应变及 位移计算点数(NS≤24); NIC——积分最多次数(NIC≤40); INTT——计算点状态参数(INTT=0,表面点,INTT≠0,内部点); DEL——近似积分的精度,常用0.0001; CR——荷载圆的半径(cm); CP——荷载的单位接触压力Kg/cm^2; R(NS)——每个计算点离荷载中心的径向坐标值(cm); Z(NS)——每个计算点离表面的垂直坐标值(cm); E(NL)——每一层的弹性模量,kg/cm^2; PR(NL)——每一层的泊桑比; HA(NH)——每一结构层的厚度(NH=NL-1)cm; 输出变量 STRESS Z——Z方向正应力 kg/cm^2; STRESS R——R方向正应力 Kg/cm^2; STRESS T——T方向正应力 Kg/cm^2; STRESS ZR——ZR方向剪应力 Kg/cm^2; STRAIN Z——Z方向应变; STRAIN R——R方向应变; STRIIN T——T方向应变; DISP Z——Z方向位移,cm。 3. 计算实例 例1、某一弹性三层体系,有关参数如下: 3,10,18 8*0,1,1 0.0001,7.07,10.65 0.,10.,20.,30.,50.,80.,150.,250.,0.,0. 8*0.,10.,50. 0.25,0.25,0.35 20000.,10000.,500. 10.,40. 4. 计算结果 STRESSES, STRAINS & VERTICAL DISPLACEMENTS IN AN ELASTIC ULTILAYER SYSTEM UNDER A SURFACE CIRCLE UNIFORM LOAD NUMBER OF LAYERS ------------- -----= 3 NUMBER OF P OINTS ------------- = 10 CONTACT RADIUS OF LOAD --------- = 7.07000 CM CONTACT PRESSURE OF LOAD ----- = 10.65000 KSC TOLERANCE FOR INTEGRATION --- = .00010 RADIUS COORDINATES AND VERTICAL COORDINATES FOR EACH CONPUTING POINT NO. R(CM) Z(CM) 1 .000 .000 2 10.000 .000 3 20.000 .000 4 30.000 .000 5 50.000 .000 6 80.000 .000 7 150.000 .000 8 250.000 .000 9 .000 10.000 10 .000 50.000 PARAMETERS FOR EACH LAYER NO. PR E(KSC) H(CM) 1 .250 20000. 10.000 2 .250 10000. 40.000 3 .350 500. NUMBER OF ITERATIONS IC : 5 RETURN= 1 STRESSE,STRAIN& VERTICAL DISPLACEMENT FOR EACH COMPUTING POINT NO. STRESS Z STRESS R STRESS T STRESS ZR 1 -10.64922 -9.76979 -9.76979 .00000 2 .00000 .40168 -2.61501 .00000 3 .00000 .01439 -1.08695 .00000 4 .00000 -.11002 -.70087 .00000 5 .00000 -.12058 -.43806 .00000 6 .00000 -.00761 -.25757 .00000 7 .00000 .06889 -.08904 .00000 8 .00000 .05590 -.02340 .00000 9 -4.02456 1.53290 1.53290 .00000 10 -.06498 .61827 .61827 .00000 NO. STRAIN Z STRAIN R STRAIN T DISP. Z 1 -.28822E-03 -.23325E-03 -.23325E-03 .21370E-01 2 .27667E-04 .52772E-04 -.13577E-03 .16424E-01 3 .13407E-04 .14306E-04 -.54527E-04 .14085E-01 4 .10136E-04 .32601E-05 -.33668E-04 .12942E-01 5 .69829E-05 -.55306E-06 -.20396E-04 .11447E-01 6 .33147E-05 .28392E-05 -.12783E-04 .96647E-02 7 .25193E-06 .45575E-05 -.53133E-05 .65413E-02 8 -.40623E-06 .30877E-05 -.18690E-05 .39669E-02 9 -.23955E-03 .10779E-03 .10779E-03 .18071E-01 10 -.37411E-04 .47994E-04 .47994E-04 .13806E-01 八.计算结果分析 (1)路基应力  图5-6路基竖向应力系数随深度、(E1/E0) 变化情况 铺设路面结构层的主要目的是扩散车轮荷载,以减少传给路基的应力值,因为过大的应力值使路基出现剪切破坏或出现塑性变形,从而使路面结构破坏。图5-6是相对刚度不同的双层体系,沿荷载截面中轴上路基竖向应力系数随深度而变化的情况。图中可以明显地看出,在路面厚度不变的情况下,随路面材料刚度的增长(E1/E0),路基的应力急剧减少,特别是路基顶面处的应力值降得更快。例如,在两层分界面处,按均质半无限体(E1/E0=1)计算所得的约为竖向应力的68%,而设置模量增大9倍的面层后,约为竖向应力的30%。 利用三层体系的数值解,可以分析基层或面层的厚度和刚度对路基顶面竖向应力的影响。面层和路基的刚度不变时,竖向应力系数随基层刚度和厚度而变化的情况下(如图5-7a)所示。可以看出,随层刚度和厚度的增加而减少。面层刚度的影响如5-7b,路基应力随面层刚度的增加而减少。面层刚度很大时,基层厚度对路基应力的影响很小。  基层的厚度和刚度 面层的厚度和刚度 图5-7基层或面层的厚度和刚度对路基顶面竖向应力的影响 由此可见,为把路基应力降到某一容许值,可以采用增加面层或基层厚度或刚度办法,其中增加刚度比增加厚度效果大。这个规律对于设计沥青路面的基层有重要的意义。采用粒料基层时,由于本身的模量值很低,只能通过增加厚度来减少路基应力;而采用刚度较大的稳定类基层,则可明显减少路基应力,并且在相同的路基类型和容许应力(弯沉)条件下,其厚度可比粒料基层减少很多。  图5-8 a)厚度影响(E2/E1= E0/E2=0.2) b)路基和面层刚度的影响(h1/(=0.5 h2/(=1.0) (2)路面弯沉 路面弯沉是路基和路面结构不同深度处竖向应变的总和。对于等级不太高的路面来说,其中70-95%由路基提供,各点的应变是三向应力状态的函数,因此,影响路基应力的诸因数也会影响路面弯沉。  图5-9 基层底面拉应力系数(E2/E1=0.2, h1/(=0.5) 图5-8给出了三层体系荷载面中轴处的表面弯沉系数随层厚和模量而变化的情况。增加面层或基层的厚度都可减少路面弯沉;但在面层或基层厚度较薄时,增加厚度对降低弯沉量的影响比层后大时显著得多。也可通过增加路基、基层或面层的刚度使路面弯沉量降低。对比图中曲线变化可以看出,在路基刚度低时,路基刚度对弯沉量的影响要比基层和面层的影响明显得多。 (3)基层底面的拉应力 采用刚度较大基层将提高荷载扩散能力,使路基的应力和弯沉量减少。但是随着基层相对刚度的增大,基层底面的拉应力增大。此拉应力如果超过材料的抗拉强度,基层将会断裂,并导致路面破坏。图5-9给出了在面层相对刚度和厚度不变时,基层底面拉应力系数随基层相对刚度和厚度而变化的情况。可以看出,增加基层的相对刚度,将导致增大,而在基层较薄时,刚度对的影响要比厚基层严重得多。因此,为降低路基的应力或路面弯沉值而选用相对刚度较大的基层时,应验算基层底面的拉应力,使材料的抗拉强度与之相适应。  图5-10面层底面拉应力系数(E2/E0=10., h2/(=2) 基层底面最大拉应力位置一般在荷载作用面中轴处;在双圆荷载作用下,则出现在其中一个荷载作用面的中轴处。 (4)面层的底面的拉应力 垂直荷载作用下,面层底面的径向应力并非都是拉应力(如图5-10)。 面层较薄而相对刚度较小时(E2/E1>0.35, h1/(<0.5,或E2/E1>0.1, h1/(<0.25),可能出现拉应力。面层较厚和刚度较大时,面层底面便出现拉应力。它随面层相对刚度的增大而增大,特别是面层的相对刚度很大时(E2/E1<0.3),拉应力随刚度的增大而急剧增大。底面最大拉应力的位置,一般在荷载面的中轴处;双圆荷载时,最大拉应力一般出现在某一荷载面的中轴处,但在面层很厚时,随层厚增增大而移向双圆荷载面的对称轴处。 在圆形均布的单向水平荷载作用下,面层内会出现较大的径向拉应力。特别在路面荷载作用面的边缘处,其数值很大(图5-11)。面层较薄时,其底面也会出现较大的径向拉应力(5-12)。  图5-11水平荷载作用的情况基层底面拉应力系数(E1/E2=5 E2/E0=10.,h2/h1=2)  图5-12 水平单向荷载作用下,面层底面的径向应力系数 (5)剪应力 增加上层的刚度,还将导致层内剪应力的增加。垂直荷载作用下,面层内任一水平面上的最大剪应力一般出现在通过荷载作用面边缘的垂直线上(图5-13)。在面层相对刚度增大时,最大剪应力随深度而变化的情况:最大应力出现在面层中部,并随面层刚度的增大而增大。但在面-基的分界面上,最大剪应力最大剪应力随面层刚度的增加而减少。面层的厚度对最大剪应力也由很大的影响。在面层和基层的相对刚度不变的情况下,随面层厚度的减少,最大剪应力增加。其最大值出现的位置逐渐上移。由二分点上升到三分点附近。因此,面层相对刚度很大而厚度较薄时,垂直荷载将产生较大的剪应力。应采取措施以防止面层出现出现较大的应力。  图5-14 水平荷载作用下径向剪应力系数 面层受到圆形均布的单向水平荷载作用时,面层内各水平面上受到的最大剪应力随深度的增加而衰减的很快(5-14)。面-基曾分界处,最大值已下降大水平力的不到10%。而在基层底面,最大剪应力已经可以忽略不计。 垂直荷载和水平荷载共同作用下,面层内最大剪应力也随深度的增加而减少,并随面层相对刚度和厚度增加而增加。 路面的最大剪应力出现在荷载面边缘处。其值主要受水平大小的影响,同时也受面层厚度和刚度的影响,当面层相对刚度较小时,面层刚度和厚度对最大剪应力的影响很小。  图5-13 a)随深度的变化 b)面层厚度的影响  图 沥青混合料压缩蠕变试验(1)σ1=30Kpa,(2)σ2=480Kpa,温度60℃,侧应力=0 九.粘弹性层状体系 1.路面材料的基本性能 弹性(线性)关系 非线性关系 粘弹性关系 塑性关系  图  图8-2 沥青劲度随时间的关系 2.概述 1955年,E.L.Lee首先发表了粘弹性体的应力分析论文; 六十年代,Monismith,Y.H.Huang等提出双层、三层粘弹性体的位移解; 七十年代,提出了路面设计体系,即设计、使用性能预估、维修对策等应一体化考虑; 八十年代,提出VESYS-IV,逐步从理论走向实际。 3.粘弹性模型 弹性体: 粘性体: 粘弹性模型: Maxwell模型  Kelvin模型:  Van de Poel模型:Kelvin串联弹簧  Burgers模型:Maxwell与Kelvin的串联  Lechersich模型:Kelvin与粘壶的串联 Jeffreys模型: Maxwell与粘壶的并联 4.弹性半空间体 弹性半空间体的一般解:  5.对应法则 要得到粘弹性问题的解,首先要求出对应的弹性问题的解,然后对解中的模量及荷载求Laplace变换,将粘弹性算子代入,再求逆变换得到粘弹性问题的解。 6.弹性半空间体粘弹性解 ①变换   ②粘弹性算子 Maxwell模型: 粘弹性算子为: Kelvin: 粘弹性算子为: ③代入①得  ④Laplace逆变换 Kelvin  t=0 w=0 t=( w=弹性体 Maxwell  t=0 w=弹性体 t=( w=无限大 7.广义本构方程 Jeffreys模型: Maxwell与粘壶的并联  ;  ;   ;  则:  同理:;   8.非线性弹性体系 (1)概述 结构计算经常采用;路面结构沥青的非线性 (2)计算方法 增量法;迭代法;混合法 (3)增量法 荷载:;位移:;应力: 设第i-1步的应力已知,根据及应力-应变关系,确定弹性矩阵,得到刚度矩阵,然后第i步的刚度矩阵为,得I步的位移增量: i=1,2,3,…..m