Page 61 - 《爆炸与冲击》2025年第6期
P. 61
第 45 卷 胡学龙,等: 考虑动态拉压比影响的岩石损伤本构模型 第 6 期
因此,模型的本构关系为:
[ ( )]
dσ ij = (1−ω)D i jmn (1− Q 1 )dε mn − Q 2 + D ijmn ε mn −ε p dω (23)
mn
(1−ω)(∂F /∂σ ij )D ijmn (∂G/∂σ mn )
′
Q 1 = ( ) (24)
(1−ω) ∂F /∂σ ij D ijmn (∂G/∂σ mn )− Q∂F /∂γ p
′
′
[ ( )]
(1−ω) ∂F /∂ω− D ijmn ε mn −ε p D ijmn (∂G/∂σ mn )
′
Q 2 = ( ) mn (25)
(1−ω) ∂F /∂σ ij D ijmn (∂G/∂σ mn )− Q∂F /∂γ p
′
′
2 模型数值实现
为了求解和验证本构模型以及增强模型的
F′ , tri
应用性,需要对本构模型进行数值实现。由于本
F′ , k+1
文中建立的本构模型是一个塑性与损伤耦合的
F′ , k Δσ ij tri, k+1
模型,对模型进行求解比较复杂,因此本文采用
σ ij k Δσ ij crector
图 4 所示的完全隐式后向欧拉算法的简化形式
k+1
σ ij
半隐式应力返回算法 [32] 实现数值计算。模型数
值实现的流程如下。
tri
k
( 1 ) 弹 性 预 测 : 计 算 试 探 应 力 σ = σ +
ij
ij
D i jmn ∆ε k+1 ,其中,k 为迭代步数。
mn 图 4 应力返回算法示意图
F (σ ,γ ,ω )≤0 是
tri
p,k
k
′
ij
(2) 判断强度屈服函数 Fig. 4 Schematic diagram of stress return algorithm
σ k+1 = σ tri ;若不存在,继续
否存在:若存在,则有 ij ij
进行第 (3) 步。
F (σ ,γ ,ω )>0 时,试探应力位于屈服面以外,需要把试探应力拉回到屈服面上。
p,k
tri
k
′
ij
(3) 应力返回:当
根据塑性力学可知,在屈服面上有:
F ′,k+1 ( σ k+1 ,γ p,k+1 ,ω k+1 ) = 0 (26)
ij
对式 (26) 在试探应力处进行一次泰勒展开有:
tri tri
∂F ′ ∂F ′
′,k+1
′,tri
F = F + ∆σ crector + ∆ω = 0 (27)
ij ∂ω
∂σ ij
[ ( ) ] tri
∆σ crector = σ k+1 −σ = − 1−ω D ijmn ∆ε −σ ij /(1−ω )
p
tri
k
k
mn
i j
ij
ij
联立式 (17)、(21)、(26) 和 (27) 可得
F ′,tri
dλ = tri tri (28)
p tri
k
k
′
′
′
′
′
(∂F /∂σ ij )(1−ω )D ijmn ∂G/∂σ mn | +{Q [(∂F /∂σ ij )σ ij /(1−ω )−∂F /∂ω]}| − Q∂F /∂γ |
ß [ ]
(Q+δ i j ∂G/∂σ ij )/ D 1 (p +a) D 2 p≤0
∗
Q = −η 2 γ p
′
η 1 η 2 e Q p>0
σ k+1 p,k+1 γ p,k+1 ω k+1 :
(4) 更新应力 ij 、塑性应变 ε i j 、有效塑性应变 和损伤变量
σ k+1 = σ +∆σ tri,k+1 +∆σ crector (29)
k
ij ij ij ij
p,k+1 p,k k+1 crector
ε ij = ε ij −C ijmn ∆σ mn (30)
p,k
γ p,k+1 = γ +∆γ p,k+1 (31)
k
ω k+1 = ω +∆ω k+1 (32)
式中:C ijm n 为弹性柔度张量。
061412-6