基于真实微结构驱动的砂岩有效弹性模量数值模拟

Numerical Simulation of Effective Elastic Moduli for Sandstone Based on Real Microstructure-Driven Mechanisms

摘要: 砂岩主要由碎屑颗粒和胶结物构成，属于典型的非均匀介质。本文基于X射线断层扫描技术获取两种砂岩样品的真实微结构，在此基础上，结合孔隙率–弹性模量半经验模型和有限元方法对砂岩样品的有效弹性模量开展数值模拟，并与超声波法测试进行比对。结果表明，本文提出的基于真实微结构驱动的有效弹性模量数值模拟方法切实可行，且具有较高的精度。

Abstract: Sandstone is primarily composed of detrital grains and cement, representing a typical heterogeneous medium. Based on X-ray computed tomography technology, this study obtained the real microstructures of two sandstone samples. Subsequently, numerical simulations of the effective elastic moduli of the sandstone samples were conducted by integrating a porosity-elastic modulus semi-empirical model with the finite element method. The results were compared with ultrasonic testing measurements. The findings demonstrate that the proposed numerical simulation method for effective elastic moduli, driven by real microstructures, is both feasible and highly accurate.

1. 引言

砂岩是自然界地表岩层中分布最广的沉积岩之一，约占沉积岩总量的三分之一，其分布范围仅次于泥岩。它主要由碎屑颗粒和胶结物组成，属于典型的非均匀介质，碎屑颗粒直径通常在0.05~2 mm之间。石英是其主要矿物成分，并含有少量长石、白云母等其他矿物，因而具有较高的结构稳定性。作为工程中常见的沉积岩，砂岩也是人类使用最广泛的天然石材之一，广泛应用于建筑装饰、雕刻艺术以及基础设施领域，如铺路和水利工程。此外，砂岩还可作为玻璃原料、磨削工具，并在能源开发中作为油气储集层发挥重要作用。

在岩石研究领域，X射线断层扫描(X-ray Computed Tomography, XCT)技术以其无损、高分辨率的特点，成为研究岩石真实微结构的有力工具[1]。通过XCT动态扫描，研究人员可以获取岩石内部孔隙分布、裂纹演化等的三维图像，再结合理论推导与数值计算，从而深入了解岩石的破裂机理以及物理性能。通常将岩石实施XCT扫描成像数据进行三维重构得到的数字化岩石模型，称为“数字岩心”。数字岩心被认为可以准确地反映岩石内部的空间结构和物质组成，为后续的仿真模拟和理论分析提供基础。例如，针对油田普遍面临的低采收率问题，很大程度上是由于注水过程中剩余油被困在多孔介质的孔隙中所导致的，Hu等[2]利用XCT技术重构数字岩心，通过数值模拟方法，系统研究了注水速度、油水黏度比、润湿性条件等因素对孔隙尺度油水流动特性及采收率的影响。Elkhoury等[3]基于XCT技术，对来源于印第安纳州的石灰岩和砂岩进行系统分析，研究中使用不同的分割算法提取孔隙结构，并分析影响结论的各个因素。Lei等[4]针对火山岩复杂的裂隙分布情况，提出交互式图像分割技术处理XCT图像，并将传统阈值分割方法获得的结果与使用深度学习技术获得的结果进行了比较。Qajar等[5]利用XCT成像技术获取数字岩心，进而依据Laplace方程估算岩石的渗透率，同时对岩石的非均质性和各向异性进行评估。值得指出的是，应用普通台式XCT扫描仪，其空间分辨率基本能够满足常规砂岩的建模需求，然而对于孔隙结构极其致密的非常规岩石，如页岩、碳酸盐岩等，台式XCT扫描仪的空间分辨率常常难以满足孔隙结构准确识别的要求，这时需要应用同步辐射XCT扫描仪才能构建孔隙尺寸更小、结构更复杂储层的三维数字岩心，其空间分辨率达到纳米级别[6]-[8]。

弹性模量属于材料力学基本性能参数，用以衡量其在承受荷载时抵抗弹性变形的能力。作为典型的多孔非均质介质，砂岩的结构–性能关系不仅取决于固体骨架，还取决于孔隙网络特征，包括孔隙尺寸与形状、空间分布与连通性等因素[9]-[19]。关于岩石多孔结构–弹性模量关系的现有研究方法主要分为三类，即经验法、解析法和数值法。经验法是通过大量的试验建立起岩石多孔结构特征参数如孔隙率与弹性模量的经验关系。解析法以有效介质理论为代表，通过对岩石孔隙或固体颗粒的尺寸与形状、空间分布与连通性等特征做出合理的设定，使其满足特定的统计规律，然后基于等价效应对土工介质整体的应力应变线弹性行为进行数学求解。如Liu等[20]基于矿物微力学性质，采用广义有效介质理论研究了岩石单轴抗压强度与弹性模量，实现了矿物–岩体力学性质的跨尺度关联。数值法主要依赖现代计算机的强大运算能力，以有限元法为代表，通过对土工介质的内部结构进行离散化处理，进而施加虚力或虚位移等，计算系统的能量，使用迭代计算直至系统的能量达到最小值，最后依据广义胡克定律确定土工介质整体的有效弹性模量。如Tang等[21]使用矿物分析仪、纳米压痕测试和原子力显微镜来测量每个样品的微观结构以及造岩矿物的杨氏模量，借助有限元方法建立了考虑矿物空间分布特征的颗粒尺度数值模型，该模型能够准确预测花岗岩的宏观弹性性能。

2. 材料与实验

2.1. 样品来源

本文砂岩样品取自四川省自贡市，采集深度处于距离地表120~150 m的范围，包括红砂岩以及灰砂岩两种，两者均属于细砂岩，最大颗粒粒径不超过0.1 mm。两种砂岩采集试样的标准尺寸为直径50 mm，高度100 mm。样品的矿物组分如表1所示。

Table 1. Mineral components of sandstones

表1. 砂岩矿物组分

砂岩类型	石英	钠长石	方解石	蒙脱石	微斜长石	伊利石	绿泥石	赤铁矿	钾长石
红砂岩	61%	14%	9%	6%	4%	4%	1%	1%	⁄
灰砂岩	50%	18%	12%	9%	⁄	6%	4%	⁄	1%

2.2. XCT实验

为便于进行XCT测试，进一步从标准试样的中心部位钻取直径8 mm、高度10 mm的样品，再手工打磨至直径4 mm、高度5 mm。岩石样品测试使用的是SKYSCAN 1272型X射线显微镜。对尺寸为直径4 mm、高度5 mm的岩石样品以100 kV的X射线峰值能量和100 μA的电流进行。样品距射线源和探测器分别为60 mm、212 mm。每次扫描获得1200帧二维XCT图像。对岩石样品，每个投影的曝光时间为2.3 s，XCT图像像素的空间分辨率为4 μm。对不同砂岩图像顺序裁取128帧二维XCT图像(512 × 512 μm²)，将其重构成三维XCT图像。图1所示为岩石样品的三维XCT重构图像(512 × 512 × 512 μm³)，其像素灰度值的变化范围为0~255。

(a) (b)

Figure 1. 3D XCT image of sandstones (a) red sandstone; (b) gray sandstone

图1. 砂岩三维XCT图像(a)红砂岩；(b)灰砂岩

2.3. 超声波实验

如图2所示，针对直径50 mm、高度100 mm的圆柱体试件，采用ZBL-U5100非金属超声波探测仪进行砂岩弹性模量的测试。该设备基于超声脉冲技术，主要包括以下四个步骤，即超声波发射：仪器内置的超声波发射器产生高频超声波脉冲，通过单通道或双通道配置将脉冲信号传输至待测试件中；超声波传播：超声波在试件内部传播过程中，遇到介质界面(如缺陷、空洞或密度变化区域)时会产生反射、折射等物理现象；信号接收与处理：接收器捕获反射波信号后，由内置处理单元对超声波的传播时间、幅度变化等特征参数进行分析，据此评估试件的内部结构特征和力学性能；结果显示：经处理后的检测数据通过高亮度彩色液晶屏实时显示，研究者可直观获取检测结果并进行后续分析。

(a) (b)

Figure 2. Ultrasonic detector for testing the elastic modulus of sandstones. (a) Sensor and coupling agent; (b) Sample test diagram

图2. 超声波探测仪测试砂岩弹性模量。(a) 传感器和耦合剂；(b) 样品测试图

超声波法测试固体材料弹性模量E的基本公式如下：

$E=\frac{v_p^2\rho \left(1+\nu \right)\left(1-2\nu \right)}{1-\nu }$ (1)

其中，v_p为超声波测试波速，ρ为试件的密度，$\nu$ 为试件材料的泊松比常数。测试装置在试件一端产生超声波信号，另一端则接收超声波信号，超声波信号的传播距离记作L (即圆柱体试件的高度)，超声波信号的传播时间记作t，则超声波的测试波速为：

$v_p=\frac{L}{t}$ (2)

将式(2)代入式(1)，则有

$\frac{L^2}{t^2}=\frac{E\left(1-\nu \right)}{\rho \left(1+\nu \right)\left(1-2\nu \right)}$ (3)

试件的密度ρ满足

$\rho =\frac{M}{V_s}$ (4)

其中，M、V_s分别为试件的质量、体积。对每一种类型的砂岩进行超声波测试，均采用5个试件为一组，对每一组的测试值取平均，测试误差均不超过5%，结果如表2所列。

Table 2. Ultrasonic velocity of sandstones

表2. 砂岩超声波波速测试值

砂岩类型	超声波波速
红砂岩	17.7 m/s
灰砂岩	10.9 m/s

3. XCT像素灰度值与微孔隙率关联方程

以往基于XCT与有限元法(FEM)的异质材料(如岩石、混凝土等)有效弹性模量研究，多采用矿物相分割的方法，即依据灰度阈值划分不同物相(如石英、长石、方解石、黏土胶结物和孔隙等)，并为各相赋予相应的本构关系。该类方法虽然使模型更接近物理真实情况，然而从数值角度看，其存在一种缺陷，即不同的灰度阈值范围没有统一的标准，阈值选择常依赖人为经验，从而引入主观误差。本文提出并采用了一种直接基于XCT图像灰度值的模拟方法，以规避此种缺陷。

本文中对于砂岩的三维XCT图像，依据像素灰度值可计算微孔隙对应的局部孔隙率。假定每个像素包含一定的微孔隙，其对应的局部孔隙率f_v定义如下：

$f_v=\frac{1-h_v}{1-\langle h_v\rangle }{f}_s$ (5)

式中，$\langle h_v\rangle$ 代表关于图像内所有像素归一化灰度的平均值，f_s为测试样品的整体孔隙率。显然，当h_v = 1时，f_v = 0；当h_v = 0时，$f_v=f_s/\left(1-\langle h_v\rangle \right)$ 以及$\langle f_v\rangle =f_s$ 。对本文所考虑两种砂岩样品的微孔隙率分布进行统计分析，发现其具有典型的峰值特征，并呈近似高斯函数，如图3所示。两种砂岩的微孔隙率分布范围基本一致，但红砂岩的微孔隙率分布区间相对集中，灰砂岩则更为分散，这与二者的矿物组成差异(灰砂岩中黏土矿物含量较高)密切相关。

(a) (b)

Figure 3. Statistical distribution of micro-porosity in sandstones. (a) Red sandstone; (b) gray sandstone

图3. 砂岩微孔隙分布统计。(a) 红砂岩；(b) 灰砂岩

4. 孔隙率–弹性模量半经验模型

从复合介质角度看，多孔介质通常被视为属于最简单的一类复合介质，即其只包含两种组分介质：固体骨架(Solid Skeleton)和孔隙网络(Pore Network)。假定孔隙网络的体积分数即孔隙率为f，则固体骨架的体积分数为1 − f，式(6)即所谓的Voigt极限，表明多孔介质的有效弹性模量E随孔隙率增加而线性降低[22]。

$M=M_0\left(1-f\right)$ (6)

其中，M₀代表固体骨架的弹性模量。基于弹性力学理论，Dewey等给出了低孔隙率极限条件下多孔介质的有效弹性模量公式[23]：

$E=1-\frac{3\left(1-{\nu }_0\right)\left(9+5{\nu }_0\right)}{2\left(7-5{\nu }_0\right)}f$ (7)

其中，${\nu }_0$ 代表固体骨架的泊松比常数。在低孔隙率极限条件下，孔隙之间的相互作用忽略不计，从而导致多孔介质的有效弹性模量与孔隙率呈线性关系。在更普遍范围的中等以及高孔隙率条件下，孔隙之间的相互作用不能忽略不计，从而导致多孔介质的有效弹性模量与孔隙率呈非线性关系。使用较为广泛的非线性关系表达式如下[24] [25]：

$\frac{M}{M_0}=1-\left[M\right]f+\left(\left[M\right]-1\right)f^{ 2}$ (8)

其中，[M]代表固体骨架的本征弹性模量，其定义式如下：

$\left[M\right]=-\underset{f\to 0}{\lim }\frac{M_r-1}{f}$ (9)

其中，M_r = M/M₀代表多孔介质的相对弹性模量。当孔隙形貌为规则球体且固体骨架的泊松比${\nu }_0=0.2$ 时，式(8)具有十分简洁的形式[26]：

$\frac{M}{M_0}={\left(1-f\right)}^2$ (10)

对于实际多孔介质难以满足孔隙形貌为规则球体且固体骨架的泊松比${\nu }_0=0.2$ 时，Gibson-Ashby半经验模型通常被用以描述实际多孔介质的弹性模量–孔隙率非线性关系，即：

$\frac{M}{M_0}=\gamma {\left(1-f\right)}^2$ (11)

其中，γ为经验参数，γ的数学表达式有不同形式可供选择，如：

$\gamma =\frac{3\left(1+{\nu }_0\right)}{4}$ (12)

$\gamma =\frac{5\left(1+{\nu }_0\right)}{6}$ (13)

除式(11)外，Spriggs [27]提出了另一种广泛使用的半经验模型，即：

$\frac{E}{E_0}={\text{e}}^{-\alpha f}$ (14)

其中，α为经验参数。不难发现，当$f\to 0$ 时，式(14)退化成：

$\frac{E}{E_0}=1-\alpha f$ (15)

本文采用Spriggs半经验模型(14)，在前期基于X射线断层扫描和纳米压痕测试的砂岩微孔隙与微弹模分布特征的研究中[28]，已获得的孔隙率–弹性模量半经验模型参数如表3所列。

Table 3. Fitted parameters of Spriggs formula in porosity-modulus relation for sandstones [28]

表3. 基于Spriggs半经验模型的砂岩孔隙–弹模关系拟合参数[28]

砂岩类型	E0	α
红砂岩	162 GPa	8.0
灰砂岩	124 GPa	8.7

5. 有限元方法

本文利用有限元方法模拟砂岩有效弹性模量的理论基础是变分原理。具体地，针对砂岩的XCT三维微观结构施加一定的全场应变，通过迭代计算使系统(三维微观结构)的能量(弹性自由能) E_n (n代表单元数目)最小，由此确定每个单元(三维XCT图像像素点)上的最终位移分布。为使系统能量E_n取极小值，需满足能量对变量u_m (节点位移)的偏导数均为零，如下：

$\frac{\partial E_n}{\partial u_m}=0$ (16)

在数值求解过程中，当系统能量E_n对第m个节点位移的偏导数构成的梯度矢量的平方和小于某一给定允许误差时，可近似认为式(16)成立，由此确定砂岩三维微观结构各单元上的应力和应变分布，并结合广义胡克定律计算有效弹性模量。

利用有限单元法解决实际问题时，通常需要用分割线将计算区域离散成网格，网格越密，即单元设置的越小，离散化后的数据越能准确表示实际研究对象，但是数据量将快速增大，进而影响计算效率，因此选择合适大小的网格对于有限元计算过程是极为重要的。节点是网格线的交汇点，连接相邻单元。针对由XCT技术测定获得的三维微观结构，本身就是由离散的体像素所构成，每个像素即可视为一个有限元单元，因此无需进行网格划分。将每一像素视为一个有限元单元，单元的形状为立方体，通过8个节点(恰好对应着单元立方体的8个角)与相邻单元连接。在下面的讨论中，网格、单元和像素三者表示的意义是相同的，一个像素就代表一个单元。

Figure 4. Schematic diagram of finite element nodes

图4. 有限元单元节点示意图

Table 4. Unit node numbering rules

表4. 单元体节点编号规则

节点编号	Δi	Δj	Δk
1	0	0	0
2	1	0	0
3	1	1	0
4	0	1	0
5	0	0	1
6	1	0	1
7	1	1	1
8	0	1	1

为了推导出每个像素点上的有限元方程，首先需要定义每个像素上8个节点的编号以及顺序。为方便起见，如图4所示，规定节点1的编号同时代表该像素在三维网格体系中的编号(i, j, k)，而另外2至8节点的编号规则如表4所示，其中(Δi, Δj, Δk)均相对于节点1而言。图中的(i, j, k)坐标轴对应着直角坐标系中的(x, y, z)轴。

单一像素有限单元体具有的弹性能EE由下式给出：

$EE=\frac{1}{2}{\displaystyle \int {\epsilon }_{mn}{C}_{mnrs}{\epsilon }_{rs}\text{d}V}$ (17)

式中，m = 1, 2, 3, n = 1, 2, 3, r = 1, 2, 3, s = 1, 2, 3表示三维笛卡尔坐标系中的3个方向分量，ε_mn、ε_rs代表应变分量，C_mnrs代表广义弹性模量。系统的总能量可以通过对所有像素单元的弹性能求和获得。对于各项同性介质，应变张量具有对称性，利用Voigt记号，则有：

$EE=\frac{1}{2}{\displaystyle \int {\epsilon }_a{C}_{ab}{\epsilon }_b\text{d}V}$ (18)

式中，ε_a、ε_b是单元内一点(x, y, z)处的局部应变；a (或b) = 1, ∙∙∙, 6代表分量标签，即(ε_xx, ε_yy, ε_zz, ε_xz, ε_yz, ε_xy)；C_ab是广义弹性模量。它们满足关系式：

${\epsilon }_{mm}=\frac{\partial u_m}{\partial x_m}$ (19)

${\epsilon }_{mn}=\frac{\partial u_m}{\partial x_n}+\frac{\partial u_n}{\partial x_m}$ (20)

为方便起见，x₂ = y，x₃ = z写法通用。u_m (x, y, z)代表该像素内各处的位移变量，其中0 < x, y, z < 1，坐标原点位于标号为1的节点。u_m (x, y, z)可依据各节点处位移通过线性插值方法获得，即：

$u_m\left(x,y,z\right)=F_{mrn}\left(x,y,z\right)u_{rn}$ (21)

式中，u_rn代表像素单元节点位移，F_mrn (x, y, z)代表线性插值系数，r (或s) = 1, ∙∙∙, 8代表有限单元体中8个节点的标签。F_mrn满足：

$F_{1r1}=F_{2r2}=F_{1r1}=F_r$ (22)

$F_{1r2}=F_{1r3}=F_{2r1}=F_{2r3}=F_{3r1}=F_{3r2}=0$ (23)

F_r (x, y, z)的具体数学表达式如下：

$F_1=\left(1-x\right)\left(1-y\right)\left(1-z\right)$ (24)

$F_2=x\left(1-y\right)\left(1-z\right)$ (25)

$F_3=xy\left(1-z\right)$ (26)

$F_4=\left(1-x\right)y\left(1-z\right)$ (27)

$F_5=\left(1-x\right)\left(1-y\right)z$ (28)

$F_6=x\left(1-y\right)z$ (29)

$F_7=xyz$ (30)

$F_8=\left(1-x\right)yz$ (31)

应变ε_a (x, y, z)与位移u_m (x, y, z)满足关系：

${\epsilon }_a\left(x,y,z\right)=S_{arn}\left(x,y,z\right)u_{rn}$ (32)

式中，

$S_{arn}\left(x,y,z\right)=L_{am}{F}_{mrn}\left(x,y,z\right)$ (33)

转换算子L的数学表达式如表5所示。

Table 5. Mathematical expressions of strain-displacement transformation operators

表5. 应变–位移转换算子数学表达式

a分量	m = 1 (x)	m = 2 (y)	m = 3 (z)
1	∂/∂x	0	0
2	0	∂/∂y	0
3	0	0	∂/∂z
4	∂/∂z	0	∂/∂x
5	0	∂/∂z	∂/∂y
6	∂/∂y	∂/∂x	0

重新整理像素单元体的弹性能则有：

$EE=\frac{1}{2}{\displaystyle \int {\left[S_{arm}{u}_{rm}\right]}^{\text{T}}{C}_{ab}\left[S_{bsn}{u}_{sn}\right]\text{d}V}$ (34)

上式可进一步改写成：

$EE=\frac{1}{2}{u}_{rm}^{\text{T}}{D}_{rmsn}{u}_{sn}$ (35)

式中，D_rmsn的表达式为：

$D_{rmsn}={\displaystyle \int {\left[S_{arm}\right]}^{\text{T}}{C}_{ab}\left[S_{bsn}\right]\text{d}V}$ (36)

6. 结果与讨论

根据XCT图像灰度，可以得到任意像素单元体弹性模量的数学表达式，即红砂岩：E = 162e^(−8.0f)以及灰砂岩：E = 124e^(−8.7f)。应用有限元方法，在每个单元上施加均匀的工程应变(ε₁, ε₂, ε₃, ε₄, ε₅, ε₆) = (0.1, 0.1, 0.1, 0.1, 0.2, 0.3)，依据共轭梯度法历经循环迭代求解满足系统总能量最小值的应力及应变分布，如图5所示、图6所示。根据三维微观结构中所有像素单元体具有的应变张量和应力张量，计算算术平均值获得有效应变〈ε〉和应力〈σ〉，进而依据广义胡克定律确定水泥石整体的有效剪切模量、体积模量、以及杨氏模量。

鉴于Spriggs半经验模型的砂岩孔隙–弹模关系拟合参数与实际值可能存在一定的误差，即参数E₀以及α的取值难以完全准确。为了考察参数的取值对岩石弹性模量数值模拟结果的影响，另外针对红砂岩选取α = 7.0进行有限元计算，以及针对灰砂岩选取α = 7.7进行有限元计算。砂岩有效弹性模量数值模拟以及超声波法测试结果如表6所示。

(a) (b)

Figure 5. Red sandstone E₀ = 162 GPa, α = 8.0, dimensionless Strain, Stress unit GPa (a) Stress; (b) Strain

图5. 红砂岩E₀ = 162 GPa，α = 8.0，应变无量纲，应力单位GPa (a)应力分布；(b)应变分布

(a) (b)

Figure 6. Gray sandstone E₀ = 124 GPa, α = 8.7, dimensionless Strain, Stress unit GPa (a) Stress; (b) Strain

图6. 灰砂岩E₀ = 124 GPa，α = 8.7，应变无量纲，应力单位GPa (a) 应力分布；(b) 应变分布

Table 6. Numerical effective elastic modulus of sandstones

表6. 砂岩有效弹性模量数值模拟

砂岩类型	参数设置(E0, α)		弹性模量模拟值	超声波测试值
红砂岩	E0 = 162 GPa	A = 7.0	21.5 GPa	17.7 GPa
		A = 8.0	16.2 GPa
灰砂岩	E0 = 124 GPa	A = 7.7	12.2 GPa	10.9 GPa
		A = 8.7	9.7 GPa

不难发现，红砂岩的有效弹性模量明显高于灰砂岩，这与两种砂岩的矿物组成密切相关。具体地，红砂岩主要以高弹性模量的石英矿物为骨架，而灰砂岩主要以低弹性模量的黏土矿物为骨架。此外，随着参数α的减小，两种砂岩对应的有效弹性模量值均增大，红砂岩增加32.7%，灰砂岩增加25.8%，说明参数α对有效弹性模量E有显著负相关性。

对比两种砂岩有效弹性模量的测试值，发现结合Spriggs经典半经验模型，红砂岩：E = 162e^(−8.0f)，灰砂岩：E = 124e^(−8.7f)和有限元方法的模拟值与测试值吻合较好，其中误差均不超过15%。这说明，Spriggs半经验模型较好地描述了砂岩微孔隙率–微弹性模量之间的数学关系；同时，基于XCT测试技术和有限元方法开展砂岩有效弹性模量的预测是可行的。相较而言，砂岩的模拟值低于超声波测试值，这种现象可能归咎于超声波反映的是材料动态弹性模量，而在模拟方法中获取的弹性模量为静弹性模量。

本文构建基于真实微结构驱动的砂岩有效弹性模量数值模拟框架，即从XCT图像出发，构建灰度值–局部孔隙率关系，利用Spriggs模型，实现有效弹性模量的有限元模拟。当前工作以4 μm的空间分辨率为参考，证明其切实可行，且具有较高的精度。值得注意的是，XCT空间分辨率对该方法具有重要的影响，如文献[29]所述，灰度值–局部孔隙率关系的适用性存在于一定的空间分辨率范围。本研究提出的模拟方法摒弃了传统的物相区分步骤，而是将其视作不同组成及堆积结构的统计平均，从而直接利用XCT图像的灰度值反映微观结构。该方法在牺牲一定细节的同时合理简化模型，显著提升了计算效率，实现了高效性与准确性的良好平衡。

7. 结论

本文基于X射线断层扫描技术获取两种砂岩样品的真实微结构，在此基础上，结合孔隙率–弹性模量半经验模型和有限元方法对砂岩样品的实际有效弹性模量开展数值模拟，并与超声波法测试进行比对。主要结论如下：

1) 红砂岩的有效弹性模量明显高于灰砂岩，这与两种砂岩的矿物组成密切相关。具体地，红砂岩主要以高弹性模量的石英矿物为骨架，而灰砂岩主要以低弹性模量的黏土矿物为骨架。

2) Spriggs半经验模型E = E₀e^(−αf)较好地描述了两种砂岩微孔隙率–微弹性模量之间的数学关系，对于红砂岩满足E = 162e^(−8.0f)，对于灰砂岩满足E = 124e^(−8.7f)。

3) 对比两种砂岩有效弹性模量的测试值，发现结合Spriggs半经验模型和有限元方法的模拟值与测试值吻合较好，其中误差均不超过15%，说明基于XCT测试技术和有限元方法开展砂岩有效弹性模量的预测是可行的。

基金项目

在此衷心感谢国家自然科学基金(项目编号：U2040222、52078123)的资助。