2.3　大体积岩土混凝土宏观力学的计算机随机模拟方法研究

2.3.1　大体积岩土混凝土宏观力学性质

随机模拟又名蒙特卡罗方法，有时亦称随机抽样技术或统计试验方法。其求解问题的基本思想是，先构造一个与问题相联系的概率模型，然后通过对该模型的观察或随机抽样试验来确定问题的解。

对于大体积岩土、混凝土，它的小尺度力学性质随机场为

式中　——小尺度岩土、混凝土弹性模量、泊松比、摩擦系数、粘结力。

由的统计均匀各向同性以及相互独立的性质，的均值与协方差矩阵分别为：

式中　——两点间的距离；

——的协方差函数；

——与的互协方差函数，余类推。

显然，欲在已知的情况下，在计算机上实现岩体的模拟，首先需将宏观岩土、混凝土特征单元аV离散成一个个代表小尺度岩体特征单元的小立方体，并且所模拟的小尺度特征单元数有限（设总数为N）。离散之后，如果将各小尺度特征单元（小立方体）的力学性质统一构成一个4 N维的随机向量z

其中（i＝1，2，…，N）为аV中第i个小尺度单元的位置矢径，则宏观岩土、混凝土特征单元аV可近似用多维随机向量W来表示。

由式（2.3-2）、式（2.3-3），W的均值νW与协方差矩阵DW分别为：

式中

［ka］表示取ka的整数为值（k＝i，j；a＝E，μ，f，c；i，j＝1，2，…，4 N）。

根据工程实践，即使我们已经测得同一材料分区内一处或若干处材料一些小尺度力学性质参数，对于该处或其中一处材料的其他小尺度力学性质参数，或者材料分区内其他任一处的材料的小尺度力学性质参数，也还是无法完全确定的。因此，对于W1-ν1，W2-ν2，…，W4N-ν4N，不存在不完全为零的常数k1，k2，…，k4N，使得：

否则，譬如设km≠0（1≤m≤4N），则由式（2.3-6）可得：

由于在一个材料分区内ν1，ν2，…，ν4N是确定的［见式（2.3-6）］，因而上式表明，当我们测得小尺度力学性质W1，W2，…，Wm-1，Wm＋1，…，W4N之后，即可由它们完全决定小尺度力学性质Wm。这显然与前述工程实践事实相矛盾。故对于任意一个非零的普通向量，有

上式表明W的协方差矩阵DW是正定的。

于是，根据线性代数理论，正定对称的协方差矩阵DW可分解为

式中H为下三角矩阵，即

其中

这样，若假设小尺度材料力学性质服从正态分布，从而随机向量W也服从正态分布，则由随机模拟方法（Monte Carlo）中关于协方差矩阵为正定的多维正态随机向量的抽样方法即可得

式中S＝［S1　S2　…　S4N］T，S1，S2，…，S4N为相互独立的N（0，1）分布的随机变量。事实上，由于

故式（2.3-11）是正确的。

式（2.3-11）表明，对多维随机向量W的模拟，亦即对宏观材料特征单元аV的模拟，最后取决于随机向量S的抽样，而S的抽样方法则由于S1，S2，…，S4N相互独立和均服从N（0，1）分布，是相当简单的：可分别对其各个向量Si独立按下式进行抽样

或

式中ui，νi是相互独立的［0，1］区间上均匀分布的随机变量。由于在一般计算机上均设有［0，1］区间上均匀分布的随机数，故ui，νi可按标准函数直接调用，无需自编算法。

因此，根据式（2.3-12）、式（2.3-13）产生的一组相互独立的N（0，1）分布的随机变量S1，S2，…，S4N，即可由式（2.3-11）产生的实现，亦即宏观材料特征单元аV的小尺度力学性质的实现。

为了求得大尺度岩土混凝土的宏观力学性质，将代表小尺度材料特征单元的小立方体作为空间有限元，并施加于аV以均布应力荷载（图2.3-1），就可以对相应于随机向量W的任一实现的宏观材料аV进行弹塑性有限元计算。由于弹塑性有限元的计算可按增量加载的方式进行，通过计算分析，不仅可以确定大尺度材料的宏观强度，同时还可以知道大尺度材料的受力—变形关系（值得注意的是，后者对工程同样具有重要意义）。

以上对大尺度材料的整个模拟过程见图2.3-1。

2.3.2　深厚覆盖层高压旋喷凝结体的宏观力学参数研究

如何利用现有的地基勘测试验资料来确定设计所用的高喷体变形模量、抗剪断摩擦系数与粘结力等力学性质参数，是工程建设的关键技术问题之一。本小节结合弹塑性有限元和随机模拟两方法，来建立高喷体宏观力学性质的计算机模拟方法。

2.3.2.1　高压旋喷凝结体的统计力学模型

由于高喷体钻孔抽芯是芯样获得率较低，室内试验成果不能宏观全面反映施工质量，仅代表胶结良好的凝结体的各种力学指标，而不代表胶结不良的凝结体或整个凝结体的各种力学指标。另一方面，室内各芯样试验成果也不均一，分散较大。针对这一试验情况，为研究高喷体的宏观力学参数，对河床深厚覆盖层地基各分层的高喷体提出如下统计力学模型：

（1）将高压旋喷凝结体分成两组材料组成：一组为成芯凝结体，另一组为不成芯凝结体。用芯样采取率p来反映成芯凝结体在所分层高喷体中所占比例。据此，不成芯凝结体所占比例q为1-p。

（2）成芯凝结体与不成芯凝结体在所分层高喷体中的空间位置相互独立地呈统计均匀分布。即这两组凝结体分别以比例p和1-p随机均匀分布与整个高喷体中。

（3）在同一分区中，成芯凝结体与不成芯凝结体的小尺度力学参数的概率特性（如均值、方差等）与空间位置和坐标轴旋转无关，即具有统计均匀性和各向同性。

（4）成芯凝结体与不成芯凝结体的小尺度力学参数呈正态或对数正态分布。

（5）从高喷体中取一包含足够多芯样尺度大小的立方体区域为高喷体的宏观特征单元（图2.3-2），该宏观特征单元的力学性质代表高喷体的宏观力学性质。

2.3.2.2　高喷体宏观力学性质的弹塑性有限元随机模拟方法

1.高喷体小尺度力学参数的随机抽样方法

用坐标面将高喷体宏观特征单元离散成一个个代表小尺度凝结体特征单元的微元体。对于每一微元体，其小尺度力学参数有E、μ、f′、c′等。当小尺度力学参数服从正态分布时，其抽样公式为

或

式中x1为小尺度力学参数抽样值；，其中E′、分别为小尺度力学参数E、μ、f′、c′的均值；σX＝σE、σμ、σf′、σc′，其中σE、σμ、σf′、σc′分别为小尺度力学参数E、μ、f′、c′的均方差；u1，u2是［0，1］区间上均匀分布的相互独立的随机变量。

当小尺度力学参数服从对数正态分布时，其抽样公式为

或

当考虑微元体内小尺度抗剪断强度参数f′与c′的相关性时，若设f′与c′的相关系数为ρ，则f′与c′的协方差矩阵C为

协方差矩阵C为正定对称矩阵，根据线性代数理论，C可分解为上下三角矩阵之积，即

式中H为下三角矩阵，即

其中hij（i，j＝1，2）分别为

于是，由蒙特卡罗方法中关于协方差矩阵为正定的多维正态随机变量的抽样方法（徐钟济，1985），可得小尺度抗剪断强度参数f′与c′相关时的抽样公式为

式中s1、s2为相互独立的标准正态随机变量［服从N（0，1）分布］，其抽样公式为

因此，根据式（2.3-14）～式（2.3-24），即可产生高喷体宏观特征单元内的各小尺度力学参数的实现或抽样。

2.弹塑性有限元方法

在弹塑性有限元计算中，材料屈服准则采用Drucker-Prager准则，塑性流动法则采用正交流动法则。所用Drucker-Prager准则的屈服圆锥面通过Mohr-Coulomb屈服六边形的外顶点，相应屈服函数F为

式中　I1——应力张量的第一不变量；

J2——应力偏张量的第二不变量；

α、K——材料参数。

I1、J2、α、K的计算公式分别为

其中σx、σy、σz、τxy、τyz、τzx为应力分量；φ＝arctanf′。

3.高喷体弹塑性有限元随机模拟方法的计算步骤

（1）用坐标面将高喷体宏观特征单元离散成一个个代表小尺度凝结体特征单元的微元体，其中每一个微元体在有限元计算中用一个空间六面体有限单元来反映（图2.3-3）。

（2）按空间位置相互独立地服从统计均匀分布和成芯凝结体在所分层高喷体中所占比例p，确定成芯凝结体与不成芯凝结体在所分层高喷体中的空间位置。具体方法为：利用［0，1］区间上均匀分布的随机变量u，对每一个微元体相互独立地进行［0，1］区间上均匀分布的随机变量u的抽样。若设第i个微元体的均匀随机变量的抽样值为ui，当ui≤p时，该微元体属于成芯凝结体；当ui≥p，该微元体属于不成芯凝结体。

（3）对成芯材料和不成芯材料的微元体，分别按式（2.3-14）～式（2.3-24）进行高喷体宏观特征单元内的各小尺度力学参数的抽样。

（4）在无水平围压（即σ2＝σ3＝0）而铅直均布面力为σ1的情况下，对小尺度力学参数抽样后的高喷体宏观特征单元进行三维有限元应力变形计算。根据宏观特征单元的铅直、水平变形计算值，可得高喷体的宏观弹性模量E宏和宏观泊松比μ宏分别为，。式中分别为宏观特征单元沿铅直向和水平向的平均压缩（拉伸）变形；ΔL为宏观特征单元的边长。

（5）对小尺度力学参数抽样后的高喷体宏观特征单元，以增量加载的方式进行渐进破坏过程的弹塑性有限元数值模拟。增量加载的方式为先施加与宏观特征单元以等量的均布荷载σ1、σ2、σ3，然后再逐步增加σ1，直至宏观特征单元完全破坏为止。其中均布荷载σ1、σ2、σ3实际为宏观特征单元的宏观主应力。

（6）根据第五步所得若干不同围压σ2、σ3下的σ1，可画出一系列强度莫尔圆。由相应于工程应力范围内的强度包络线，即可得高喷体的宏观摩擦系数f′宏和宏观粘结力c′宏。