弹粘塑性孔隙介质在冲击荷载作用下的一种本构关系 ――第一部分：状态方程

岩石力学与工程学报 22(9)：1405～1410

2003年9月 Chinese Journal of Rock Mechanics and Engineering Sep.，2003

弹粘塑性孔隙介质在冲击荷载作用下的一种本构关系

——第一部分：状态方程

戚承志 王明洋 钱七虎

(北京建工学院土木工程系 北京 100044) (解放军理工大学 南京 210007) (中国工程院土木水利与建筑学部 100038)

*

摘要 利用弹塑性介质的不可逆热力学的普遍关系，对弹粘塑性孔隙介质在冲击静水压力作用下的力学行为进行了研究。在本构关系中，使用了对于固体来讲更为合理的德拜形式的自由能。利用等效应变法，考虑了孔隙率对孔隙介质在冲击静水压力作用下的力学行为的影响，并得到了孔隙率演化方程。最终得到的本构关系符合不可逆热力学定律，形式相对简单，物理意义明确。数值计算结果与已有结果吻合很好。

关键词 粘弹性、粘塑性力学，本构关系，弹粘塑性孔隙介质，不可逆热力学，等效应变法 分类号

TU 435 文献标识码 A 文章编号 1000-6915(2003)09-1405-06

CONSTITUTIVE RELATION OF ELASTOVISCOPLASTIC POROUS

MEDIA SUBJECTED TO SHOCK

——PART I：STATE EQUATION

Qi Chengzhi1，Wang Mingyang2，Qian Qihu3

(1Department of Civil Engineering，Beijing Institute of Civil Engineering and Architecture， Beijing 100044 China)

(2PLA University of Technology， Nanjing 210007 China) (3Chinese Academy of Engineering， Beijing 100038 China)

Abstract The process of deformation and fracture of deformable solid media subjected to shock is accompanied by many irreversible dissipative mechanical，physical and chemical processes，and the constitutive equations of solid subjected to shock should express these irreversible features. The mechanical behaviors of elastoviscoplastic porous media subjected to hydrostatic shock are studied here in the framework of irreversible thermodynamics of elastoviscoplastic media. The corresponding constitutive equations are expressed by Helmholtz free energy of Debye form，as it is more rational for solid. The influence of porosity on the mechanical behaviors of elastoviscoplastic porous media subjected to hydrostatic shock is remarkable，and here it is considered by means of effective strain method. The evolution equation of porosity is also obtained on this basis. The constitutive equations presented here are relatively simple，and their physical meaning is clear，corresponding to irreversible thermodynamics. In order to validate the proposed constitutive relations，a numerical example is given，and the agreement is good with the existing results.

Key words viscoelastic and viscoplastic mechanics，constitutive relation，elastoviscoplastic porous media，irreversible thermodynamics，effective strain method

2002年1月8日收到初稿，2002年3月20日收到修改稿。

* 国家自然科学基金(A50179038)及教育部回国留学人员启动基金(教外司留[1999]747号)资助项目。

作者 戚承志 简介：男，1965年生，1986年毕业于工程兵工程学院，1992年于前苏联圣彼得堡建筑工程大学获工学博士学位，现任副教授、解放军理工大学在站博士后，主要从事结构抗震、随机振动、破坏力学方面的教学和研究工作。

• 1406 • 岩石力学与工程学报 2003年

1 引 言

在诸如核爆及化爆的较强作用下，变形介质中会发生相应的力学、热学及构造变化。这种构造变化一般包括硬化程度、相变程度、化学反应程度、晶格缺陷程度、位错密度、损伤积累程度等等[1

～3]

，

它们在一般情况下是相互作用的。介质的这种变化要求在介质的本构关系中引入内变量以描述这种变化[4]。介质的动力破坏是一个复杂的、多阶段的过程，包括微损伤的产生、发展及连接、胚胎裂纹的形成，连接及发展成宏观裂纹。宏观裂纹的最终发展导致材料的最终破坏分裂。通常破坏可分为延性破坏、脆性破坏及绝热剪切带的形成[5]。

介质与结构在爆炸冲击荷载作用下的动力响 应是一个令人感兴趣的问题。在爆炸源附近，应力偏张量通常远小于应力球张量，因此，通常材料强度的细节可以忽略。但当冲击波球面向外传播而接近自由表面时，强度效应变得十分明显。

岩石等介质的非线性变形与破坏通常与微损伤的形成及晶粒的流动有关。岩石介质通常为孔隙介质。岩石的孔隙度本质上与岩石构成的颗粒性有关，它不仅用于度量微孔洞，也用于度量微裂纹。当压应力足够高时，介质会发生压缩，孔隙度减少，同时，由于微结构的破坏，介质内产生损伤，使介质的抗压及抗拉强度减少[6]。

对孔隙材料的本构关系，有些模型是基于连续介质力学建立的，如文[7～9]，另外一些是基于混合物理论而建立的，如文[10]。目前，随着损伤破坏力学的发展，固体材料的损伤等耗散因素都被考虑进去，这样，材料行为的耗散性及不可逆性都被表现出来。这种耗散性及不可逆性必须符合不可逆热力学的基本定律。因此，用热力学法研究此类问题就显得更合理。基于热力学理论框架的孔隙介质模型，需要确定介质的Helmholtz自由能函数。一些研究者基于最终的微分结果，唯象地来反推Helmholtz自由能，这种构造方法带有人为的和牵强的感觉，而不是基于物质内在的微观构造机制。在文[11，12]中，利用介质的不可逆热力学理论及统计物理的流体爱因斯坦(Einstein)形式的自由能函数，给出了介质的本构关系。为了修正与实际介质之间的差异，又引入了一个孔隙率变化函数。这一模型的优点是利用介质的统计物理学的基本理论，并符合不可逆热力学基本定律。缺点是引入了太多的参数来描述介质的孔隙度及塑性变形的演化方程，参数共有60多个，应用起来不够简便。

根据固体物理理论，固体的德拜(Debye)模型比爱因斯坦模型能更好地反映固体的真实情况。尽管德拜模型是基于简单物质推导而来，但计算表明，适当选取参数，德拜模型能成功地应用于较复杂介质的行为描述。因而，本文的目的是基于热力学定律和德拜形式的自由能函数，利用有效应变法，给出孔隙介质在爆炸冲击作用下弹塑性区域的状态方程，形式相对简洁而描述又足够精确。本文不涉及爆炸源附近的蒸发及相变区域，因而，所采用的应变为工程应变εij。

2 弹粘塑性介质的不可逆热力学

介质的不可逆过程可以用热力学第二定律表述。热力学第二定律可用Clausius-Duhem不等式表述如下[13]：

ρ0TS

&(n)≥0 (1) 式中：ρ0为介质密度；T为绝对温标；S&(n)为熵产的不可逆部分，可表示为

ρ&(n)=σ&(n)T0TSijεij−ρ0ψ&−ρ&Δ0TS−T

qr

≥0 (2) 式中：

ψ为Helmholtz自由能，S为单位质量熵，ΔT为温度梯度，qr为热通量，ε&(n)

ij为变形速度。

通常假定热传导过程不取决于局部的热力学过程，这样，可以把式(2)分成两个独立的不等式：

局部熵产不等式：

ρ0TS&(n)=σijε&(n)ij−ρ0ψ&−ρ0T&S≥0 (3) 热传导不等式：

ΔTT

qr

≥0 (4) 把应变速度分为可逆部分与不可逆部分之和：

ε&=ε&(e)+ε&(n) (5) 上式等号右边第一项为可逆部分，第二项为不可逆部分。

假设Helmholtz自由能为状态函数，依赖于可逆变形、内变量αi及温度T，即ψ=ψ(ε(e)，αi，T)，

则由式(1)可得

⎛

⎜σ∂ψ⎞⎜ij−ρ0(e)⎟⎟ε&(e)+σρ∂ψijε&(n)−0⎝

∂εij⎠

∂αα&i− i⎛⎜∂ψ⎞⎝ρ0

∂T+S⎟⎠

T&−ΔTr

Tq≥0 (6) 第22卷 第9期 戚承志等. 弹粘塑性孔隙介质在冲击荷载作用下的一种本构关系 • 1407 •

对可逆过程，αi=0，ε(n)=0，考虑到上式括

号中的表达式为状态函数，不取决于ε&(e)及S，因此，可得如下本构关系：

σ∂ψ⎫

ij=ρ0

∂ε(e)ij⎪⎪

⎬ S=−ρ∂ψ (7) ⎪

0

∂T⎪⎭

引入与αi共轭的广义力：

β∂ψi=−ρ0

∂α (8) i

则式(6)变为

σΔTijε&(n)ij−βiα&i−T

qr

≥0 (9) 要组成完备的描述物质不可逆热力学行为的方程组，除式(7)外，还需补充内变量αi的演化方程。内变量可以是损伤、不可逆塑性变形等。

3 弹粘塑性孔隙介质的自由能函数及本构关系

从上述第2节可以看到，要确定介质的本构关系，关键的一步是确定单位质量的自由能。

很多作者根据已有的弹塑性本构关系，凭经验及直觉，以应变、损伤变量、塑性硬化参数等为变量，构造Helmholtz 自由能。某些实验表明，损伤仅与弹性变形有关，而与塑性变形无关，这样，可以把Helmholtz 自由能表示成与可逆变形ε(e)及损伤变量D有关的项，即与ψe(ε(e)，D，T)及与塑性变形有关的量ψp(εp，T)：

ψ=ψe(ε(e)

，D，T)+ψp(εp，T) (10)

首先要找出介质基体材料的Helmholtz 自由能的表达式。

在文[11]中，借鉴了流体的爱因斯坦形式的自由能。但这种形式的自由能采取了爱因斯坦假设，即晶体的谐振子系统所有振动模式都相同。德拜模型考虑了谐振子系统振动模式的分布，因而能更好地反映现实[14]。这样式(10)中的ψe可表示成如下形式：

ρm0ψe=ρm0ψmv(εv，T)+2G(εv，T)J′2

(11) ψmv=CmT[3ln(1−e−x)−D(x)]+em0(εv)−ηm0T (12)

式中：εv=εij，为介质的体积应变，以压为正；

ρm0为密实介质的密度；

G为剪切弹性模量；

Cm为密实介质比热的最大值； em0为介质的零点能，依赖于εv；

J′2

为第二应变偏量不变量； ηm0为零温零弹性畸变形时的熵常量；

D(x)=3∫ x

y3

x

3

0ey−1dy，为德拜函数； x(ε)=θD(εv)

vT

，θD(εv)为依赖于εv的温度

参数；

G(εv，T)为材料的剪切模量，依赖于εv及

T。

介质的应力可表示为静水压力P与偏应力Sij

之和，即

σij=−Pδij+Sij (13)

而压力由下式决定：

P=(1−φ)Pm (14) Pm=Pmv+Pm

′ (15) PΨmv

mv=−ρ∂m0

∂ε= v

−ρ⎡dθDde⎤m0⎢⎣3CmTθD(x)+m0dε⎥=

Ddεvv⎦−ρ⎡γm(εv)

de0⎤m0⎢⎣3CmTJD(x)+m⎥ (16)

mdεv⎦P′=− 2∂G

m

∂εJ′2

(17) v

式中：Jm=1−εv；φ0为初始孔隙率；γm为基体材料的Gruneisen系数，

dθDγ(εθdε=mv)

(18) DvJm

偏应力张量为

Sij=2Gεij(1−φ) (19)

熵由下式决定：

η=ηm=ηmv+η′ (20)

ηmv=−

∂Ψmv

∂T

=− Cm[3ln(1−e−x)−D(x)]+3D(x)+η0 (21)

η′∂G

m

=−2∂T

J′ (22) • 1408 • 岩石力学与工程学报 2003年

内能由下式决定：

e=Ψ+Tη=emv+e′m

(23) emv=Ψmv+Tηmv=em0+3CmTD(x) (24)

e′1⎛m

=ρ⎜2G−T∂G⎞

⎟⎠

J′ (25) m0⎝∂TG(εv，T)按文[15]取下式： G(ε2

v，T)=G0[1+A1/1J

mPmv

−A2(T−T0)] (26)

式中：G0为材料参考构形的剪切模量；A1，A2为材料常数。

Pmv和emv中含e0(εv)，因此，需要确定。可利

用密实介质的Hugoniot冲击绝热线数据PmH，emH来确定e0(εv)。PmH，emH依赖于εv。

以冲击绝热线为参考线的Mie-Gruneisen 物态方程为

Pmv=PmH+ρmγm(emv−emH) (27)

把式(15)，(24)代入上式，有

dem0γm(εv)γ(ε)dε−eP

m0=SH−mvemH (28) v1−εvρm01−εv

这样，联立式(17)，(28)即可确定θD及e0。

为简化计算，方程(17)中的γm与参考构形中的

γm0之间可取如下关系

[16]

：

γm=γm0(1−εv) (29)

代入式(18)，有

dθD

θ=γm0 Ddεv

θD=θD0exp(γm0εv) (30)

θD0=θD(0)，可由参考构形中的零应变时的比

热确定。

F(x)可由表查出[17]。这样，由Cm0可确定θD0。 C=∂emv∂T=3C⎡⎢θD(0)⎤

m0

mF(x0)=3CmF⎣T⎥⎦

(31) 式中：F(x)=3

x

y4ey

x

3

∫ 0(ey−1)2dy。

这样，只需解方程(28)，就可确定e0(εv)，进而确定其他值。

把式(29)代入方程(28)得

dem0dε−γP

m0em0=mH−γm0emHv

ρ (32)

m0其初始条件为

em0(0)=emH(0)−3CmTD(x0) (33)

如果忽略强度效应，并取无应力状态为参考态，则内能由下式确定：

emH=

12P11

mH(v0−v)=2PmHv0εv=2

PHεv/ρ0 (34) 这样，式(32)变为

de0dε−γeP⎛1⎞

m00=mH⎜1−γm0εv⎟ vρ⎝2⎠

(35)

m0可用数值法求解方程(35)。

计算表明，在干孔隙材料零应力条件下由式(14)所计算出来的体积模量数值比实验所测到的体积模量数值要大，也即由式(14)所预测的材料性质比实际材料的刚度大，材料有额外的压缩性。在文[11]中，假定在弹性响应中孔隙率会变化，给出了孔隙率随体积变形的变化函数，以模拟这种材料的额外弹性可压缩性对材料响应的影响，过程比较复杂，需要的参数较多。模拟结果大应力区域与实验结果相符很好，而开始阶段差别较大。

本文采用有效应变法模拟这种现象，即对于一定的名义应变ε，由于孔隙材料受压而使孔隙材料的骨架产生变形(包括弹性变形及非弹性变形)而引起局部应力集中，当所受的应力足够大时，会引起材料骨架塌陷，从而使材料的有效应变比名义应变要小。为此，采用如下的有效应变公式：

εε⎧a⎡b(1−φ20)εv⎤

⎫e=v⎨1−⎩

bexp⎢⎣−φ0(1−ε2⎥⎬+εe0 (36) v)⎦⎭式中：a，b为常数，

a

b

≤1.0；φ0为参考构形的材 料孔隙率；εe0为附加项。

在模拟材料初始阶段无骨架塌陷时，有效应变相对较大时，采用下列函数：

⎧=⎪c⎨φn

ε0εv，

0≤ ε v≤ ε 0e0

⎪φn

(37) ⎩c0ε0，

εv≥ ε 0 式中：c为常数，ε0为初始段有效应变修正范围，n为实数。

函数(36)，(37)的特点是当孔隙率为零时，名义应变与实际应变一样；随着体积应变的增大，有效

第22卷 第9期 戚承志等. 弹粘塑性孔隙介质在冲击荷载作用下的一种本构关系 • 1409 •

应变趋近于体积应变。对孔隙材料，可以通过调整a，b，c，ε0来拟合实验数据。c，ε0对压缩曲线的

初始阶段影响大，而a，b对应力较大区域影响大。这样，对于一定的名义应力，先计算出有效应变，然后把有效应变代入到式(14)～(35)计算应力。

通过有效应变法，可以非常直观地用名义应变与有效应变之差：

δε=εv−εe

来近似估算孔隙率的变化。但是，在这一差值之中，包含由于骨架弹性变形所引起的孔隙率变化，因此，要从此差值中去除由于骨架弹性变形所引起的孔隙率变化。现用下式来近似模拟变形的不可逆部分：

δε⎡

n=δε⎢1−exp⎜⎛⎜⎝−

β(1−φ0)εv⎞⎣

φ0(1−εv)⎟⎤

⎟⎠⎥ (38) ⎦

式中：β为常数，可由卸载实验数据来确定。这样，可以用下式来计算剩余孔隙率：

ϕ=εv−δεn (39)

卸载曲线可在P-

ε

v曲线上用直线连接卸载起

始点与εv轴上残余变形点近似获得。

P=P∗(ε−δεn)(ε∗−δεn) (40)

式中：P∗，ε∗分别为卸载起始时的压力及体积变形的大小。这种处理方法显然能简化计算。

假设孔隙率的变化主要是由非弹性变形引起的，那么，介质的体积变形率就可与孔隙率的变化率相联系，即

ε&v=φ&(1−φ) (41) 这样，在冲击静水压力作用下，公式(9)变为

−Pφ&(1−φ)−ΔTT

qr≥0 (42)

上式左侧两项在冲击静水压力作用下都大于零，故公式成立，因此，本模型满足热力学第二定律。

4 计算实例

确定弹塑性孔隙介质的本构关系需要求解方程(32)，以确定e0(εv)。要确定e0(εv)必须知道方程(32)

中的参数。取文[11，12]中Mt. Helen山的凝灰岩的参数。

材料参考构形孔隙率φ0=0.38，密实状态密度

ρm0=2.32×103 kg/m3，参考构形的Gruneisen系数

γm0=0.69，Cm=1.247×103 J/kg/K，θD0=1 000 K，参考构形的温度T0=300 K，参考构形的德拜函数值D(x0)=D(3.33)=0.249。

采用文[3，7]中的模拟Hugoniot冲击绝热线，及温度变化规律，计算中用多项式进行拟合。不考虑剪切对压力的影响，即取式(26)中的系数A1和A2为零。方程(35)用改进的欧拉法求解。方程(36)，(37)中的参数取如下值：a=0.397，b=0.4，c=0.058 85，ε0=0.05，n=2，此时，εe0=0.005。

从计算结果(图1)可以看出，在加载的初始阶段，拟合的程度比文[11]的结果要好。从图2看出，应力较高区精度与文[11]相当。孔隙率演化如图3所示，亦与文[11]非常接近，尽管所采用的参数比文[11]的少许多。说明本模型在精度上是满足要求的，而复杂性方面，比文[11]的模型更易于实际应用。孔隙率随压缩变形的演化如图3所示，与文[12]结果相比较可以看出，二者非常接近，这说明本文采用的计算方法是可取的。另外，本方法物理意义明确，计算简单，故便于实际应用。

图1 材料压缩曲线的初始段

Fig.1 The initiative stage of compression curves of material

图2 静水压缩曲线

Fig.2 Hydrostatic compression curves

• 1410 • 岩石力学与工程学报 2003年

图3 孔隙率随体积压缩变形的演化

Fig.3 Evolution of porosity with volumetrical compression

deformation

5 结 论

在冲击荷载作用下固体介质的变形与破坏过程

为一个不可逆的耗散过程，因此，固体介质的本构关系应反映这些不可逆特征。本文利用弹塑性介质的不可逆热力学的普遍关系，对弹粘塑性孔隙介质在冲击静水压力作用下的力学行为进行了研究。德拜形式的自由能对于固体更合理，因此，在本构关系中使用了它。利用等效应变法考虑了孔隙率对孔隙介质在冲击静水压力作用下的力学行为的影响，并得到了孔隙率演化方程。最终得到的本构关系符合不可逆热力学定律，形式相对简单，物理意义明确。数值计算结果与已有结果吻合很好。

参 考 文 献

1

李守巨，刘迎曦. 爆炸荷载作用下裂纹扩展过程分析[J]. 岩石力学与工程学报，1998，17(增)：888～891 2

杨 军，金乾坤，高文学等. 岩石爆破损伤模型研究的几个问题[J]. 岩石力学与工程学报，1999，18(3)：255～258 3

黄明利，唐春安，朱万成. 岩石破裂过程的数值模拟研究[J]. 岩石力学与工程学报，2000，19(4)：468～471 4 Качановл. M. Oсновы механики раз –рушения [M]. Mосква ：

Hаука，1974

5

Perzyna P. Constitutive modeling of dissipative solid for localization and fracture[A]. In：Localization and Fracture Phenomena in Inelastic Solids[M]. New York：Springer，1998，99～241 6

Caroll M M，Holt A C. Static and dynamic pore collapse relations for ductile porous materials[J]. J. Applied Physics，1972，43(2)：1 626～1 635 7

Goodman M A，Cowin S C. A continuum theory for granular materials[J]. Arch. Rational Mech. Anal.，1972，44：249～266

8 Киселев С П. Структура ударных волн в пористых упругопла

стических материалах[J]. ПМТФ，1998，39(6)：27～32 9 Герасимов А В，Кректулев Р А，Модельдефор -мирования и. Pазрушения много компонентной пористой упругопластической среды с непрер -ывным измене нием физикомеханических харак -теристик [J]. Проблемып Рочности ，1999，(2)：139～150

10 Drumheller D S. A theory for dynamic compaction of wet porous

solid[J]. Int. J. Solid and Structures，1987，23(2)：211～237 11

Rubin M B，Elatta D，Attia A V. Modeling additional compressibility of porosity and the thermomechanical response of wet porous rock with application to Mt. Helen tuff [J]. Int. J. Solid and Structures，1996，33：761～793 12

Rubin M B，Vorobiev O Yu，Glenn L A. Mechanical and numerical modeling of porous elastoviscoplastic material with tensile failure[J]. Int. J. Solid and Structures，2000，37：1 841～1 870 13

Кукуджанов В Н，Сантаойя К. . Термодинами-ка вязкопластических сред свнутренными параметрами [J]. Мех -аника Твердого Тела，1997，(2)：115～125

14

Ландау ЛД，Лифшиц Е. М. Статистическая физика，ч.1[M]. Москва：Наука，1995

15 Steinberg D J，Cochran S G，Guinan M W. A constitutive model for metals applicable at high-strain rats[J]. J. Applied Physics，1980，51(3)：1 498～1 504

16 高占鹏. 爆炸与冲击中的实际物态方程[J]. 力学进展，1991，21(2)：176～188

17 经福谦. 实验物态方程导论[M]. 北京：科学出版社，1999