2.4 岩石流变模型的数值分析方法
2.4.1 引言
一般地说，黏性物体的时变效应不仅与当时的应力水平有关，而且还取决于整个应力历史过程。因此，为了分析黏性体的时变，必须采用增量法，并以合适的时间步长逐步地进行计算。由此可见，每一时步的黏性应变增量的计算十分重要。
把时间划分为一系列时段，设在tn时刻和tn+1时刻的黏性应变率分别为
和
，则Δtn=tn+1-tn间隔内的黏性应变增量
为
其中β为不大于1的常数。
（1）β=0时，则
。即黏性应变增量决定于Δtn间隔开始时的黏性应变速率，这是欧拉时进格式，或称为前差分式或完全显式。
（2）β=1时，则
。即黏性应变增量决定于Δtn间隔终了时的黏性应变速率，这是后差分式或完全内含式。
（3）β=1/2时，则
。即黏性应变增量决定于Δtn间隔开始和终了时的黏性应变速率的平均值，这是平均式或 “内含梯度”式。
2.4.2 显式流变分析
2.4.2.1 弹-黏塑性的力学模型
由2.3.1节的分析，将一维状态下的应力σ和应变ε改为复杂应力状态下的应力σ和应变ε，则总应变速率为
式中：
为总应变速率；
为弹性应变速率；
为黏塑性应变速率。
应力速率由式（2.4.3）计算
式中：D为弹性矩阵。
在单向受力时，开始产生黏塑性变形的条件为σp-σs=0。在复杂应力状态下，开始产生黏塑性变形的屈服条件可表示为
式中：σp为圣维南元件（或称摩擦滑块）承受的应力；σs为屈服应力；F为屈服函数；F0为常数。
当式（2.4.4）成立时，便会产生黏塑性应变。至于黏塑性应变速率的数值，与当时的应力状态有关。对于硬化材料，还可能与黏塑性应变值有关。为了简化计算，通常假定黏塑性速率只与当时的应力状态有关，用函数关系表示为
据波兹纳 （P.Perzyna）假设，黏塑性应变率
与瞬时应力之间有如下关系式
由式（2.4.2）得到波兹纳本构方程为
式中：Q为塑性势，Q=Q（σ，εvp，k）；γ为流动系数；Φ（F）为屈服函数的函数。
符号〈 〉的意义如下：
曾经有人提出过不同的函数Φ，目前最常用的两种为
和
式中：α、N分别为两个常数；F、F0意义同前。
有时简单地取为
如果采用关联流动法则，即假定Q≡F，则有
其中
对于屈瑞斯卡屈服函数
对于米赛斯屈服函数
对于莫尔-库仑屈服函数
对于德鲁克-普拉格（D-P准则）屈服函数：C1=α；C2=1.0；C3=0
其中，θ为洛德角，其为
这样，式（2.4.6a）的黏塑性应变速率完全可以求出。例如，对于符合德鲁克-普拉格屈服准则的材料，式（2.4.6a）可以写成
其中
2.4.2.2 弹-黏塑性的有限元公式和计算步骤
由2.4.1节的分析可知，用有限单元法计算黏性应变增量有前差分式（或称显式）和后差分式（或称隐式）等解法。显式解法，又称初应变法。该解法程序编制简单，所需计算机内存较小，它的缺点是时间步长受到稳定条件的限制。隐式解法虽然可以采用较大的时间步长，但程序编制复杂。本节研究显式解法。
弹-黏塑性的有限元公式为
根据材料性质，选定某一屈服准则，例如德鲁克-普拉格屈服准则，就可以从式 （2.4.13）中计算得到
，然后在时间上进行离散，使
，应用时间步长法就可计算出每一时步的εvp，从而连续解方程式 （2.4.14）。具体计算步骤归纳如下：
（1）设已知第n次时步的σn，δn，
，Rn，由式 （2.4.6a）和式（2.4.11）计算出
。
（2）近似确定εvp的增量改变值
如果外荷载Rn在该时刻也有变化，取为Rn+1；又
其中
（3）解方程（2.4.14），计算出δn+1
（4）由几何方程确定εn+1，由物理方程
求得σn+1 ，则tn+1时步的σn+1、δn+1、
都已知，再可重复计算下一时步。依次类推，可求得应力和应变随时间变化的发展曲线。
计算的精度主要取决于式（2.4.16）的近似积分。为了改善精度，可以增加一次迭代循环。计算出
后，重新令
再重复步骤（2）～（4），计算出新的
与原来的值比较，若达不到规定的精度，可以继续迭代，直至满足要求为止。不过在实际编写程序中，迭代会增加许多附加的存储量，因而增加计算时间。在实用上，可限定时间步长Δt，使计算结果得以满足收敛和稳定条件，并保证一定的精度。
时段Δt的取值是直接影响精度的重要因素，倘若Δt的取值过大，可能导致解的漂移或不稳定。为此必须给出Δt的限制条件，基于差分的稳定性理论，对于显式算法，Cormeau （1975）曾给出时段Δt上限的理论值。对于关联黏塑性流动问题Q≡F和线性函数Φ（F）=F/F0，时段长度的限制如下：
屈瑞斯卡材料
米赛斯材料
莫尔-库仑材料
德鲁克-普拉格材料，取下面两式中的小值
式中：γ为流动系数；μ为泊松比；φ为内摩擦角。
弹-黏塑性问题有限元分析方法也适用于黏弹性问题，只需要将有关公式中的εvp视为εv即可。此外，弹-黏塑性问题的有限元法也可用于弹塑性有限元分析。进行计算时，施加一个常荷载并对时间进行积分，直至所有的应变率均为零，得到的便是弹塑性解（Owen等，1980；卓家寿，1996）。也可以利用这一规律性对弹-黏塑性程序的合理性进行验证。
2.4.2.3 基于等效流变应变率的显式算法
假设固体的总应变可以分为弹性应变和与时间相关的蠕变应变，即
由胡克定律可以得到微分形式的应力与应变的关系为
在蠕变过程中的平衡方程为
式中：BT为应变位移矩阵；σ为应力张量；f为节点荷载矢量。
当采用小应变假设时，将式（2.4.24）对时间微分，并忽略惯性的影响，有
将式（2.4.23）代入式（2.4.25），得
如果在tm时刻位移值u（tm）、应力σ（tm）、外力f（tm）和蠕变应变率
）已知，平衡方程可以在Δtm时间增量内积分得到位移增量和应力增量等。
类比塑性流动理论的概念，假设蠕变应变率为
式中：Q为能量耗散函数，通常与屈服函数F相同，即采用关联流动法则。μ由蠕变应变与时间曲线的斜率确定，其中
确定方向。
这样式（2.4.28）可改写为
式中
为等效蠕变应变率
式中：A为常系数；
、T、t分别为等效应力、等效蠕变应变、温度和时间，其一般由试验得到。
此时，对于式（2.4.27）中的Δεc可以由式（2.4.29）通过下式得到
式中：
分别为tm时刻的等效应力、等效蠕变应变。
注意在式（2.4.27）的左端采用的劲度矩阵与线弹性分析时的一样。蠕变项完全处理成伪荷载，因此是一种显式方法，又称为初应变法。显然，此法的优点是不需要重新组装和求解劲度矩阵，而此法的缺点是步长不能太大，否则不能得到稳定的解。实践证明，为了保证计算得到稳定的解，蠕变增量一般不能超过总弹性应变的一半。
2.4.3 隐式流变分析
在隐式蠕变分析中，与显式蠕变分析类似，蠕变应变率与等效蠕变应变率的关系为式（2.4.29）。对等效蠕变应变率进行泰勒（Taylor）展开，此时蠕变应变率为
式中：
为等效应力。
假设在时间间隔Δt内应力σ和时间存在线性关系，蠕变应变率对时间进行积分得到
胡克定律的增量表达式为
利用Sherman-Morrison关系，应力增量可以进一步表示为
其中
此时由平衡方程［即式（2.4.25）］，容易得到隐式蠕变有限元求解的控制方程。按时间步逐步求解，即可得到应力应变随时间变化的曲线。