Biot多孔弹性模型的公式推导过程
流体连续性方程达西定律
流体质量守恒方程一般形式的biot方程Compressible
流体的压缩性
β
\beta
β定义为：为流体密度
ρ
f
\rho_f
ρf 对压力 𝑝 的变化率
β
=
∂
ρ
f
∂
p
\beta=\frac{\partial \rho_f}{\partial p}
β=∂p∂ρf
我们关心的是随时间的变化率，因此我们将变量
ρ
f
\rho_f
ρf 和压力 𝑝 都对时间
t
t
t 求导：
∂
ρ
f
∂
t
=
∂
ρ
f
∂
p
∂
p
∂
t
\frac{\partial \rho_f}{\partial t}=\frac{\partial \rho_f}{\partial p}\frac{\partial p}{\partial t}
∂t∂ρf=∂p∂ρf∂t∂p
则带入
β
\beta
β得到：
∂
ρ
f
∂
t
=
β
∂
p
∂
t
\frac{\partial \rho_f}{\partial t}=\beta \frac{\partial p}{\partial t}
∂t∂ρf=β∂t∂p
考虑孔隙流体的质量守恒，流体连续性方程描述了流体在多孔介质中的流动情况。以下是详细推导过程：
流体连续性方程
在多孔介质中，流体的质量守恒方程可以表示为：
∂
(
n
ρ
f
)
∂
t
+
∇
⋅
(
ρ
f
v
f
)
=
0
\frac{\partial (n \rho_f)}{\partial t} + \nabla \cdot (\rho_f \mathbf{v_f}) = 0
∂t∂(nρf)+∇⋅(ρfvf)=0 其中：
n
n
n 是孔隙度。
ρ
f
\rho_f
ρf是流体密度。
v
f
\mathbf{v_f}
vf 是流体的平均速度。
为了简化计算，假设流体密度
ρ
f
\rho_f
ρf是常数，即流体是不可压缩的。则方程简化为：
ϕ
∂
ρ
f
∂
t
+
ρ
f
∂
n
∂
t
+
ρ
f
∇
⋅
v
f
=
0
\phi \frac{\partial \rho_f}{\partial t} + \rho_f \frac{\partial n}{\partial t} + \rho_f \nabla \cdot \mathbf{v_f} = 0
ϕ∂t∂ρf+ρf∂t∂n+ρf∇⋅vf=0
ρ
f
∂
ϕ
∂
t
+
ρ
f
∇
⋅
v
f
=
0
\rho_f \frac{\partial \phi}{\partial t} + \rho_f \nabla \cdot \mathbf{v_f} = 0
ρf∂t∂ϕ+ρf∇⋅vf=0
将流体密度 (\rho_f) 提取出来，得：
∂
n
∂
t
+
∇
⋅
v
f
=
0
\frac{\partial n}{\partial t} + \nabla \cdot \mathbf{v_f} = 0
∂t∂n+∇⋅vf=0
达西定律
达西定律描述了流体在多孔介质中的流动：
v
f
=
−
k
μ
∇
p
\mathbf{v_f} = -\frac{k}{\mu} \nabla p
vf=−μk∇p 其中：
k
k
k是渗透率。
μ
\mu
μ 是流体的黏度。
p
p
p 是孔隙压力。
将达西定律代入流体连续性方程中：
∂
n
∂
t
−
∇
⋅
(
k
μ
∇
p
)
=
0
\frac{\partial n}{\partial t} - \nabla \cdot \left( \frac{k}{\mu} \nabla p \right) = 0
∂t∂n−∇⋅(μk∇p)=0
流体质量守恒方程
单位体积内流体质量的变化率等于流体密度变化率和孔隙度变化率的和
这一点与多元函数求偏导加和是一致的
n
n
n为孔隙度，
ρ
f
\rho_f
ρf为流体密度，单位体积多孔介质中的流体质量为
n
ρ
f
n\rho_f
nρf
则：
∂
(
n
ρ
f
)
∂
t
=
n
∂
ρ
f
∂
t
+
ρ
f
∂
n
∂
t
\frac{\partial (n\rho_f)}{\partial t}=n\frac{\partial \rho_f}{\partial t}+\rho_f\frac{\partial n}{\partial t}
∂t∂(nρf)=n∂t∂ρf+ρf∂t∂n
将
∂
ρ
f
∂
t
=
β
∂
p
∂
t
\frac{\partial \rho_f}{\partial t}=\beta \frac{\partial p}{\partial t}
∂t∂ρf=β∂t∂p带入
得：
∂
(
n
ρ
f
)
∂
t
=
n
β
∂
p
∂
t
+
ρ
f
∂
n
∂
t
\frac{\partial (n\rho_f)}{\partial t}=n\beta\frac{\partial p}{\partial t}+\rho_f\frac{\partial n}{\partial t}
∂t∂(nρf)=nβ∂t∂p+ρf∂t∂n
将达西定律中得
∂
n
∂
t
\frac{\partial n}{\partial t}
∂t∂n带入得：
∂
(
n
ρ
f
)
∂
t
=
n
β
∂
p
∂
t
+
ρ
f
∇
⋅
(
k
μ
∇
p
)
\frac{\partial (n\rho_f)}{\partial t}=n\beta\frac{\partial p}{\partial t}+\rho_f\nabla \cdot \left( \frac{k}{\mu} \nabla p \right)
∂t∂(nρf)=nβ∂t∂p+ρf∇⋅(μk∇p)
一般形式的biot方程
（Zienkiewicz, et al., 1980; M. B. C. Ulker and M. S. Rahman, 2009） σ_ij=σ_ij^‘-δ_ij p σ_ij^‘是有效应力，σ_ij是总应力，孔压p，假设ui是土骨架的位移向量，小应变能被定义为： ε_ij=1/2(u_(i,j)+u_(j,i)) 有效应力的变化定义为： dσ_ij^’=D_ijkl (dε_kl-dε_kl^0) D_ijkl是弹性张量，ε_kl^0是蠕变和不均匀应变 用拉梅常数λ和剪切模量G表示的有效应力表达式为： σ_ij^’=λδ_ij ε_kk+2Gε_ij λ是拉梅第一常数 G是拉梅第二常数 其与泊松比和弹性模量（杨氏模量）的关系为：
G=E/(2(1+μ)), λ=Eμ/((1+μ)(1-2μ)), λ=2Gμ/(1-2μ)
则： σ_ij^'=2Gμ/(1-2μ) δ_ij ε_kk+2Gε_ij
根据Terzaghi原理，有效应力与孔隙压力之间的关系来表征土骨架二维平面应变模型的平衡方程，由以下公式给出 (∂σ_x’)/∂x+(∂τ_xz’)/∂z=∂p/∂x （2） (∂∂τ_xz’)/∂x+(∂σ_z’)/∂z=∂p/∂z （3）
对于流体位移w_i，其时间偏导是平均速度，w ̇_i=(∂w_i)/∂t 体积流量除以总面积 单位总体积的平衡方程被写为 σ_(ij,j)+ρg_i=ρu ̈_i+ρ_f w ̈_i ρ是混合密度, ρ=(1-n) ρ_s+nρ_f，n是孔隙度， 对于流体相的平衡方程可以写为：： -p_(,i)+ρ_f g_i=ρ_f u ̈_i+ρ_f/n w ̈_i+(ρ_f g)/k w ̇_i k是渗透系数， 如果流体具有体积压缩性K_f 则： ε ̈_ii+w ̇_(i,i)=-(p_(,i) n)/K_f
密度变化可以表征为压力变化的函数 (dρ_f)/dt=ρ_f/K_f dp/dt 则 〖dρ〗_f=ρ_f/K_f dp 〖dρ〗_f/dp=1/K_f ρ_f
流体相对固体骨架的速度 w=n(v_f-u ̇) ∇w=-n(dn/dt+n/ρ_f (dρ_f)/dt) dn/dt≈(dε_v)/dt=ε ̇_ii=∇u
∇w=-n(ε ̇_ii+n/K_f p ̇) w ̇_(i,i)=-n(ε ̇_ii+n/K_f p ̇)
总应力σ_ij与有效应力σ_ij^‘和孔压p的关系为： σ_ij=σ_ij^’+δ_ij p 对于各向同性线弹性材料：
总应力σ_ij与应变ε_ij的 σ_ij=2με_ij+λε_kk δ_ij 则有效应力σ_ij^‘与应变ε_ij的关系为： σ_ij^’=2με_ij+λε_kk δ_ij-δ_ij p
ε_kk=ε_xx+ε_yy+ε_zz
写成非张量形式为： （1）正应力分量 σ_xx=2με_xx+λ(ε_xx+ε_yy+ε_zz)-p σ_yy=2με_yy+λ(ε_xx+ε_yy+ε_zz )-p σ_zz=2με_zz+λ(ε_xx+ε_yy+ε_zz)-p （2）切应力分量 对于切应力，i≠j, δ_ij=0，则： σ_xy=2με_xy σ_yz=2με_yz σ_zx=2με_zx
Compressible
如果孔隙水中完全不含空气，则