Analysis on the stress intensity factor(SIF)for mixed-mode crack atbottom base layer with weight functions expressed with matrix
-
摘要: 由于基层材料的变异性、基层干缩(尤其是水泥基基层)、温度应力、土基压实不够、施工和养护不当,以及混合型荷载等都会在基层形成混合型裂缝,这些裂缝很容易形成反射裂缝,因而分析基层混合型裂缝应力强度因子就非常有必要。根据能量准则、叠加原理、贝蒂互换定理等推导出用矩阵权函数计算混合型应力强度因子的方法。并分析了待定权函数系数的求法,即结合有限元计算混合型裂缝应力强度因子方法求出待定的权函数系数,继而得到了矩阵型权函数。并利用有限元检验得到的权函数,两种不同方法在不同荷载作用下的计算结果吻合得很好,表明该方法可行。并编制了相应的程序来实现上述理论,对研究反射裂缝扩展和工程应用非常有益。Abstract: Due the variations of materials, shrinkage of base layer(especially the cementitious binder is used), temperature stresses, poorly compacted soil, poor constraction and poor maintenance and the mixed-mode loads subjected to the pavement, mixed-mode cracks are easy to occur at the bottom of the base layer of pavements, so it is necessary to analyze the stress intensity factors for the mixed-mode cracks at the bottom of base layer.The method to calculate the mixed-mode stress intensity factors is derived with the matrix weight functions according to the energy rule, superpostition theory and Betti's reciprocity theory.The method to get the constant coefficients of weight functions is also analyzed, which is based on the SIFs by finite element method.The weight functions are checked by the finite element method, and the agreement appears to be fairly good for the whole results by the two methods, which reveals the weight function method is very efficient.Thecorresponding programs arealso made, and it will be useful for analysis on the reflecting cracks in pavements and the engineering applications.
-
Key words:
- weight function /
- matrix /
- base layer /
- stress intensity factor(SIF) /
- finite element
-
由于基层材料的变异性、基层干缩(尤其是水泥基基层)、温度应力、土基压实不够而导致承载力下降和施工接缝,以及施工和养护不当等都会在基层形成裂缝。裂缝在车荷载和环境荷载作用下扩展,继而反射到面层,形成反射裂缝。通常,裂缝方向大致是竖直的,但由于材料的变异性和混合荷载作用下,裂缝也有可能是倾斜的;并且,车荷载与裂缝的相对位置不同,都将车荷载看为Ⅰ型裂缝荷载是很不恰当的,必须对基层混合型裂缝进行深入研究。分析裂缝应力强度因子的方法很多,如有限元、边界元、能量法等,但对于道路特殊的层状结构,权函数方法[1,2]有独特之处,它不仅精度高,计算简便,而且用权函数分析裂缝的扩展寿命非常简洁,可以减少大量的计算工作。笔者也利用权函数方法对沥青路面Ⅰ型裂缝应力强度因子进行了相应的研究[3~5]。
1. 基本理论
对于如图 1所示裂缝,在不同荷载作用下,裂缝扩展Δa释放的能量可表示为
G1Δa=1H[(KⅠ,1)2+(KⅡ,2)2]Δa=12∫Γ→t1→Δu1dΓ (1) G2Δa=1H[(KⅠ,1)2+(KⅡ,2)2]Δa=12∫Γ→t1→Δu2dΓ (2) H={E/(1−v2)平面应变E平面应力 (3) 式中:下标“1”、“2”分别为荷载状态;G为能量释放率;E和v分别为模量和泊松比;$\vec t$和$\overrightarrow {\Delta u} $分别为应力矢量和位移增量矢量。将两种状态叠加得
1H[(KⅠ,1+KⅠ,2)2+(KⅡ,1+KⅡ,2)2]Δa=12∫Γ[→t1(→Δu1+→Δu2)+→t2(→Δu2+→Δu1)]dΓ (4) 利用贝蒂互换定理(Betti’s reciprocal theory),即
∫Γ→t1→u2dΓ=∫Γ→t2→u1dΓ (5) 当裂缝扩展Δa时
∫Γ→t1(→u2+→Δu2)dΓ=∫Γ→t2(→u1+→Δu1)dΓ (6) 所以
∫Γ→t1→Δu2dΓ=∫Γ→t2→Δu1dΓ (7) KⅠ,1KⅠ,2+KⅡ,1KⅡ,2=H2∫Γ→t1→∂u2∂adΓ (8) 由式(8)可知,要得到任意荷载状态下的混合型裂缝应力强度因子,就应该找到两种参考荷载状态,根据上式建立方程组联立求解。下面用两种参考荷载具体表示为
KⅠ(a)KⅠ,1(a)+KⅡ(a)KⅡ,1(a)=H2∫a0[σ(x)∂v1(x,a)∂a+τ(x)∂u1(x,a)∂a]dx (9) KⅠ(a)KⅠ,2(a)+KⅡ(a)KⅡ,2(a)=H2∫a0[σ(x)⋅∂v2(x,a)∂a+τ(x)∂u2(x,a)∂a]dx (10) KⅠ(a)=H2D(a){∫a0[KⅡ,2(a)∂v1(x,a)∂a−K(1)Ⅱ,1(a)∂v2(x,a)∂a]σ(x)dx+∫a0[KⅡ,2(a)⋅∂u1(x,a)∂a−KⅡ,1(a)∂u2(x,a)∂a]τ(x)dx} (11) KⅡ(a)=H2D(a){∫a0[KⅠ,2(a)∂v2(x,a)∂a−KⅠ,2(a)∂v1(x,a)∂a]σ(x)dx+∫a0[KⅠ,1(a)⋅∂u2(x,a)∂a−KⅠ,2(a)∂u1(x,a)∂a]τ(x)dx} (12) D(a)=KⅠ,1(a)KⅡ,2(a)−KⅡ,1(a)KⅠ,2(a) (13) 只要D(a)≠0,上述方程组才有唯一解。
同样可利用权函数的基本理论[1~6],将上述公式用矩阵表示如下
(KⅠ(a)KⅡ(a))=∫a0(hⅠσ(x,a)hⅠτ(x,a)hⅡσ(x,a)hⅡτ(x,a))(σ(x)τ(x))dx (14) 权函数用矩阵表示为
(hⅠσ(x,a)hⅠτ(x,a)hⅡσ(x,a)hⅡτ(x,a))=H2D(a)⋅(KⅡ,2(a)∂v1(x,a)∂a−KⅡ,1(a)∂v2(x,a)∂aKⅡ,2(a)∂u1(x,a)∂a−KⅡ,1(a)∂u2(x,a)∂aKⅠ,1(a)∂v2(x,a)∂a−KⅠ,2(a)∂v1(x,a)∂aKⅠ,1(a)∂u2(x,a)∂a−KⅠ,2(a)∂u1(x,a)∂a) (15) 若仅考虑Ⅰ型裂缝,由式(9)或式(10)便可得到Ⅰ型应力强度因子公式
K=∫a0σ(x)h(a,x)dx (16) h(a,x)=HKⅠ,i∂vi(a,x)∂a (17) Petroski H J、Achenbach J D在分析Ⅰ型裂缝应力强度因子时,即将William关于裂缝尖端位移的序列级数[7]近似表达为
u(a,x)=σ0H√2[4F(a)a1/2(a−x)1/2+G(a)a−1/2(a−x)3/2] (18) 所以
∂u(a,x)∂a=σ0√2H{2√aF(a)(a−x)−1/2+[4∂F(a)∂a√a+2F(a)√a+3G(a)2√a].(a−x)1/2+[∂G(a)∂a1√a−G(a)2a3/2]⋅(a−x)3/2} (19) 再对裂缝长度a微分后,代入式(17),上述权函数实际上近似表达为
h(a,x)=√2πa[D0(1−ρ)−1/2+D1(1−ρ)1/2+D2(1−ρ)3/2] (20) 式中:ρ=x/a,x轴以裂缝嘴为原点,方向沿裂缝线方向;系数D0、D1和D2是关于裂缝长度a的常数。其实,Petroski H J和Achenbach J D只取了William裂缝尖端位移的序列系数的前两项来近似表示裂缝的实际位移。由于只取了两项,估计会在一些特殊情况产生误差,所以将权函数用无穷级数表示会更精确,如式(21)所示
h(a,x)=√2πa∞∑n=0Dn(1−ρ)(n−1)/2 (21) 同理,也可将混合型权函数用无穷级数表示为
hⅠσ=√2πa∞∑n=0DⅠσ,n(1−ρ)(n−1)/2 (22) hⅡσ=√2πa∞∑n=0DⅡσ,n(1−ρ)(n−1)/2 (23) hⅠτ=√2πa∞∑n=0DⅠτ,n(1−ρ)(n−1)/2 (24) hⅡτ=√2πa∞∑n=0DⅡτ,n(1−ρ)(n−1)/2 (25) 当然,实际应用不可能取无数项,为此,本文取前面四项进行分析,并且可以证明
{DⅠσ=DⅡτ=1DⅡσ=DⅠτ=0,(n=0) (26) 如图 1裂缝,在距裂缝尖点ε(ε/a→0)作用有一对垂直裂缝的单位集中力P,因此应力强度因子可应用无限长裂缝计算公式计算,并且可不考虑几何不对称对应力强度因子的影响,所以这对法向集中力产生的应力强度因子为
KⅠ=√2πεP,KⅡ=0 (27) σ(x)=Pδ(x−x0),τ(x)=0,δ(z)={1(z=0)0(z≠0) (28) 所以
KⅠσ=∫a0σ(x)hⅠσ(x,a)dx=PⅠσ(x0,a)=hⅠσ(x0,a)KⅡσ=hⅡσ(x0,a)=0} (29) 将式(29)代入式(22)、(23),由于(ε/a→0),所以
√2πaDⅠσ,0(εa)−1/2=√2πε√2πaDⅡσ,0(εa)−1/2=0} (30) 所以
DⅠσ=1,DⅡσ=0 (31) 同理,在靠近裂缝尖端作用一对裂缝切向集中力Q,即
τ(x)=Qδ(x−x0) (32) 可以得到DⅡτ,0=1,DⅠτ,0=0。
为了求出其它未知权函数系数,可在裂缝面上不同位置作用单位法向集中力和单位切向集中力,将法向集中力和切向集中力按照式(28)和(32)用δ函数表示,再代入式(14)和(22~25),可得到
KⅠσ=hⅠσ(x,a)=√2πa∞∑n=0DⅠσ,n(1−ρ)(n−1)/2KⅡσ=hⅡσ(x,a)=√2πa∞∑n=0DⅡσ,n(1−ρ)(n−1)/2KⅠτ=hⅠτ(x,a)=√2πa∞∑n=0DⅠτ,n(1−ρ)(n−1)/2KⅡτ=hⅡτ(x,a)=√2πa∞∑n=0DⅡτ,n(1−ρ)(n−1)/2} (33) 由于未知权函数系数为三个,总共12个未知系数,所以只需要三对单位法向集中力和三对单位切向集中力。只要能得到这些集中力的应力强度因子,便可联立方程组求出权函数系数。
2. 有限元方法
从公式(33)可得,只要能计算出任意位置的集中力作用下的应力强度因子,便可以联立方程组得到矩阵权函数余下的系数。为此,下面应用有限元方法来进行辅助分析。
同前文[5]一样,裂缝尖端单元采用单参八节点变异单元,如图 2、3所示,有限元计算混合型应力强度因子公式为[8~12]
KⅠ=2μκ+1√2πl(4Uθ|B−Uθ|C)KⅡ=2μκ+1√2πl(4Ur|B−Ur|C)} (34) 式中:Uθ和Ur为垂直裂缝面和沿裂缝面的位移。这里所需要特别强调一点是,在计算Uθ和Ur时,一定要扣除裂缝尖端节点的位移。
三对集中力位置选在ρ=0、0.4和0.8。利用有限元计算出各单元的节点位移,再转化为裂缝法向和切向位移,代入式(34)可得各集中力产生的混合型应力强度因子KⅠ和KⅡ。为了减少误差,分别计算裂缝上下面的应力强度因子,再取平均值,即
KⅠ=(KⅠ,upper+KⅠ,down)/2KⅡ=(KⅡ,upper+KⅡ,down)/2} (35) 应力强度因子结果如表 1所示。沥青和基层厚度分别为15 cm和30 cm,模量分别为2000、15000 MPa,土基模量为35 MPa。三层的泊松比分别为0.25、0.25和0.35。裂缝长度为12 cm,角度为75°。整个路面简化为平面应力计算。
表 1 有限元计算的混合型应力强度因子 (MPa·cm1/2 )Table 1. Mixed-mode SIFs by finite element method将计算的应力强度因子代入公式,联立求解方程组(33)可得到权函数系数矩阵,如表 2所示。权函数如图 4所示。
表 2 混合型权函数系数Table 2. Coefficients of mixed┐mode weight function3. 检验权函数系数
尽管得到了权函数,但是否可行,还需要进一步检验,下面利用公式(33)对所得到的权函数进行检验。为此,本文对三种荷载条件进行了分析,即在裂缝面上作用均匀的法向单位应力、切向单位应力和法向单位应力与切向单位应力组合作用,分别简称Load1、Load2和Load3。利用权函数计算应力强度因子可参照文献[1]的方法。应力强度因子计算结果如表 3所示。
表 3 有限元和权函数方法的应力强度因子计算结果 (MPa·cm1/2)Table 3. SIF results by finite element method and weight function method4. 结语
从裂缝扩展的能量释放原理、叠加原理和贝蒂互换定理以及William关于裂缝尖端的位移级数序列推导出了混合型裂缝的权函数计算方法,并利用简单的荷载形式,即单位集中力作用下的应力强度因子建立联立方程组,计算出了混合型裂缝的矩阵权函数。并分析了长度为12 cm,角度为75°的基层底裂缝的权函数和应力强度因子。与有限元方法对比分析,该方法得到的权函数计算的混合型应力强度因子计算结果和有限元方法吻合得非常好,说明该方法可行。由于该方法简单,对于分析反射裂缝的扩展很有好处,非常适合工程应用。所有过程都编制了相应的程序来实现,可直接应用于工程实践。
-
表 1 有限元计算的混合型应力强度因子 (MPa·cm1/2 )
Table 1. Mixed-mode SIFs by finite element method
表 2 混合型权函数系数
Table 2. Coefficients of mixed┐mode weight function
表 3 有限元和权函数方法的应力强度因子计算结果 (MPa·cm1/2)
Table 3. SIF results by finite element method and weight function method
-
[1] BÜECKNER H F. A novel principle for the computation of stress intensity factors[J]. ZAMM, 1970, 50 (9): 529-546. [2] RICE J R. Some remarks on elastic crack-tip stress fields[J]. Int. J. Solids Structure, 1972, 8 (10): 751- 758. [3] LUORui, HUANG Xiao-ming, LU Yuan, ZHANG Wang. Research on the influence by base layer to the stress intensity factor of pavement with continual layers[J]. Journal of Southeast University, 2001, 31 (5): 61-64. [4] LUO Rui, HUANG Xiao-ming. Research on the method to calculate the stress intensity factor of asphalt pavement's edge crack[J]. Journal of Highway and Transportation Research and Development, 2002, 19(1): 12-15. [5] LUO Rui, HUANG Xiao-ming. Calculation on the stress intensity factor for the bottom crack of asphalt pavement's asphalt layer considering partial constraint by weight function[J]. Chinese Journal of Geotechnical Engineering, 2001, 23 (5): 610-614. [6] PETROSKI H J, ACHENBACH J D. Computation of the weight function from a stress intensity factor[J]. Engng. Frac. Mech., 1978, 10(2): 257-266. doi: 10.1016/0013-7944(78)90009-7 [7] WILLIAMS M L. Stress singularities resulting from various boundary conditions in angular corners of plates in extension[J]. J. Appl. Mech., 1952, 19 (4): 526-528. doi: 10.1115/1.4010553 [8] MANU C. Quadratic isoparametric elements as tran- sition elements[J]. Engng. Frac. Mech., 1986, 24 (4): 509-512. doi: 10.1016/0013-7944(86)90224-9 [9] CHONG Rhee H. Stress intensity factor evaluation from displacements along arbitrary crack tip radial lines for warped surface flaws[J]. Engng. Frac. Mech., 1989, 32(5): 723-732. doi: 10.1016/0013-7944(89)90168-9 [10] LESLIEBanks-sills, SHERMAN D. Comparison of methods for calculating stress intensity factors with quarter-point elements[J]. Int. J. Frac. 1986, 32 (2): 127-140. doi: 10.1007/BF00019788 [11] ZHANG S. FEM analysis on mixed-mode fracture of CSM-GRP[J]. Engng. Frac. Mech., 1986, 23 (3): 521-535. doi: 10.1016/0013-7944(86)90160-8 [12] CHONG Rhee H, MAMDOUH M Salama. Mixedmode stress intensity factor solutions of a warped surface flaw by three-dimensional finite element analysis[J]. Engng. Frac. Mech., 1987, 28 (2): 203- 209. doi: 10.1016/0013-7944(87)90214-1 -