DEM analysis of macro- and micro-mechanical behaviors of cemented sand subjected to undrained cyclic loading
-
摘要: 天然或人工胶结的存在能够提高砂土的抗液化能力,从宏微观尺度对其动力学性质进行研究具有重大意义。将已有的三维完整胶结接触模型引入到三维离散元程序中,对胶结砂土不排水循环三轴剪切试验进行三维离散元模拟,研究颗粒间胶结、循环应力比对离散元试样宏微观力学性质的影响。研究结果表明,胶结的存在能够抑制轴应变和孔压的发展,提高砂土的抗液化强度,循环应力比与液化振次之间具有指数函数关系,证实了本文离散元模拟能够反映胶结砂土的宏观动力学性质。在微观尺度上,当循环应力比较小时,胶结试样内部仅有极少量胶结发生破坏,力学配位数基本不变,外界输入功主要用于增加颗粒和胶结弹性能。对于特定胶结程度的试样,在初始液化发生之前,随循环应力比增加,试样内部胶结破坏更为剧烈,力学配位数下降速率更快,颗粒和胶结弹性能更快地趋向于0,颗粒摩擦耗能、弯转耗能、扭转耗能更快地达到最大值,而破坏胶结接触点、胶结接触点和无胶结接触点法方向的空间分布更快地趋向于各向同性性质。Abstract: Natural or artificial cementation formed between sand particles can strengthen the liquefaction resistance of sand. Hence, it is significant to investigate the dynamic behavior of cemented sand at the macro- and micro-scale. By introducing an existing three-dimensional (3-D) complete bond contact model into a 3-D distinct element method (DEM) code, a series of undrained cyclic triaxial shear tests on the cemented sand are simulated, where the effects of cementation and cyclic stress ratio are studied. The results show that the inter-particle cementation can restrain the development of axial strain and pore pressure, and increase the liquefaction resistance. In addition, there is an exponential relationship between the cyclic stress ratio and the number of cycles to trigger the initial liquefaction. These findings confirm the reliability of the DEM simulation in this study. When the value of cyclic stress ratio is relatively small, within the cemented specimen with given degree of cementation, a very small amount of bonds break, the mechanical coordination number remains almost unchanged, and the input work is mainly used to increase the elastic energy at the particle contacts and bond contacts. For the cemented sand before the initial liquefaction, as the cyclic stress ratio increases, the inter-particle bonds break more intensely, the mechanical coordination number declines faster. Likewise, as the cyclic stress ratio increases, the elastic energy at the particle and bond contacts tends to zero faster, and the dissipated energy reaches the maximum value in a shorter period of time. In addition, the contact normal orientation tends to be isotropic more rapidly.
-
0. 引言
地下水渗流问题涉及到石油工程、采矿工程、岩土工程、水利工程和环境工程等多个工程应用领域，其中，在给定渗流参数的情况下精确求解达西流速场（正向问题）和基于渗流实测数据合理识别渗流参数（反向问题）对于渗流模拟分析具有重要意义[1-3]。
对于正向问题，鉴于物理试验的局限性，达西渗流流速的求解主要采用数值模拟，具体方式分为两种[2]：①首先对地下水水头进行分析，再基于达西定律求解渗流流速[1-2, 5]；②建立达西渗流定律和渗流模型的耦合方程，通过联合求解同时得到渗流流速和水头[6-7]。两种方式都已被证明能够进行准确的流速计算，但却难以依据渗流数据对渗流参数进行反演。关于渗流参数的反演，目前主要基于水头或流速数据，采用最速下降法、高斯-牛顿法、遗传算法等优化算法[8]进行求解，此类方法虽对均质渗流参数具有较好计算表现，但当地下水流为非均质流情况时（渗流参数呈现空间变异性），却难以获得渗流参数的空间分布特点。因此，发展一种既能计算非均质渗流模型的流速场，又能反演其渗流参数的新算法具有重要意义。
近年来，随着人工智能的快速发展，基于深度学习万能逼近定理[10]建立的物理信息神经网络（Physics-Informed Neural Networks, PINNs）算法受到了学者的广泛关注，其核心原理：利用自动微分技术(Automatic Differentiation, AD)[11]将偏微分方程嵌入到神经网络的损失函数当中，再通过优化算法对损失函数进行寻优，从而实现神经网络对方程解的近似，具有不受数据量限制、求解格式稳定等特点。此外，相对于有限元法（FEM）、有限差分法（FDM）等，PINNs算法可通过在损失函数中加入实测数据的损失项，对方程的未知物理参数（如渗流方程中的均质或非均质渗流参数）进行反演，甚至可通过较少数据直接挖掘物理规律并构建出相应的物理方程，故在岩土工程研究中具有一定的应用前景。目前，PINNs算法已在一些领域得到了较广泛的应用，如Haghighat等[12]将PINNs算法引入到固体力学中，基于Python语言的Keras开发了SciANN框架，对带缺陷的平板拉伸变形行为进行分析，得到了平板内各点的位移。基于SciANN框架，唐明健等[13]对矩形薄板的正反问题进行了分析，准确给出了矩阵薄板各点的挠度，并利用薄板挠度反演得到了边界的受力情况。Raissi等[14]通过将PINNs算法应用于流体力学当中，准确得到了Burgers方程、退化Navier-Stokes方程的数值解，并实现了方程中流动参数的识别。Chen等[15]基于实测数据，利用PINNs算法在大量备选的偏微分算子中，寻找合适的算子，实现了Burgers方程的智能建立。PINNs算法在岩土工程领域也有一些探索，Zhang等[16]引入了独热编码向量，与PINNs算法相结合，对多孔介质的大变形固结问题进行了分析。Bekele等[17]采用PINNs研究了多孔介质中的流体流动与固体变形之间的相互作用。兰鹏等[18]采用PINNs算法对连续排水边界的一维固结问题进行了分析，得到了土体的超孔隙水压力消散规律，并实现了连续排水边界界面参数的准确反演。然而，上述的PINNs算法的应用存在如下限制，一方面，PINNs算法在进行反问题分析时，学者采用的是单物理场神经网络模型，此模型只能对未知参数为定值的情况进行反演，当反演的参数呈现空间变异性时，则无法进行识别。因此，对于非均质地下水渗流问题，采用单物理场神经网络模型PINNs算法无法获得非均质渗流参数分布，此时需引入多物理场神经网络方法[19-22]；另一方面，上述应用中，PINNs算法在进行求解时，对于边界的处理主要是软边界形式。
事实上，PINNs算法软硬约束边界的设置会对计算结果产生较为显著的影响，原因在于PINNs算法在进行求解时，需要构造方程、初始和边界条件的损失函数。初始和边界条件的损失函数主要是通过初始和边界上随机选取样本点及其取值获得，所以在进行优化时，只能保证在选取的样本点处拟合结果较好，无法保证在初始和边界上均满足条件，故而属于一种软约束形式。针对软约束这一问题，Lu等[23]和Sun等[24]引入辅助函数对其进行修正，使PINNs计算的结果自动满足边界或者初始条件，此时边界约束形式被称为硬约束形式，相比于软约束，具有更高的计算精度。目前，PINNs算法对于硬约束的应用较少，主要集中在如下方面，Lagari等[25]总结了不同类型的边界条件的PINNs硬约束构造方法，随后，Lu等[26]和陆至彬等[27]分别将硬约束引入到了超材料设计和平面板导热问题分析中，结果表明，硬约束可以获得更高的求解精度。但是在地下水渗流领域，PINNs硬约束应用暂未形成一定的研究规模。
本文将基于PINNs算法，通过施加硬约束边界条件，解决渗流模型的渗流流速计算和渗流参数反演问题。在正向问题中，对比了两种PINNs方式求解渗流流速（渗流方程和达西定律耦合，先求解水头再进行自动微分）的计算精度；在反向问题中，分别构造了单物理场或多物理场神经网络模型对均质或非均质渗流参数进行反演。最后，为了验证PINNs算法在地下水渗流问题应用的可行性，对均质和非均质渗流的若干算例进行分析，表明硬约束PINNs算法对水头进行自动微分求解渗流流速计算精度更高，且能够准确反演不同分布方式的渗流参数。
1. 原理和算法
1.1 渗流控制方程
对于计算域为Ω⊂Rn的稳定渗流问题，其控制模型为[4]
−∇T(K∇H(x))−W=0 (x∈Ω)。 (1) 式中，∇为梯度算子，x=（x1,x2,⋯,xn）为空间坐标，K为渗流参数张量（对于非均质渗流可表示为x的函数），H为地下水水头值，W为源汇项。
相关的边界条件（BC）可以表示为
B(H(x))=0(x∈∂Ω)。 (2) 式中，∂Ω为计算域的边界。基于达西定律，可得到V(x)的表达式：
V(x)=−K∇H(x)。 (3) 1.2 PINNs算法求解达西流速分布和渗流参数
（1）神经网络
运用PINNs算法进行渗流流速求解时，首先需构造神经网络（常见的为前馈神经网络）对水头或流速做近似估计。前馈神经网络由神经元连接而成，分为输入层、隐藏层、输出层3部分。各神经元输入间通过相应的权重和偏置进行线性连接，进而传输到非线性激活函数σ中得到其输出，具体计算方式如下所示，
g(z)=σ(Σpj=1wjzj+b)。 (4) 式中，z=[z1,z2,⋯,zp]为神经元的输入，g(z)为神经元输出，wj为输入zj的所占的权重，b为对该神经元所引入的偏置。一般地，常见的激活函数σ包括Tanh=(ex-e-x)/(ex+e-x)、ReLU=max{0, x}和Leaky ReLU=max{ax, x}(a为任意常数)等。本文所求解的控制方程，需要计算二阶偏导数，而ReLU和Leaky ReLU的二阶导数为0，因此选择能够无限次可微Tanh作为激活函数。
（2）自动微分（AD）
运用PINNs算法对渗流流速进行求解时，主要利用AD技术计算网络输出函数对输入变量的导数。其中，AD是一种介于符号微分和数值微分的方法，其计算原理是将复杂解析函数分解成一系列初等运算的组合，通过符号微分进行求导，再代入相关数值，通过计算机存储中间结果，最后利用链式法则得到期望导数值[11]。因此，与数值微分相比，AD能够替代网格尺度的差分运算，避免了数值微分在求解过程中方程离散带来的计算误差，提高了计算精度，且自动微分的方式主要针对神经网络的输出进行，这使梯度的计算是精确的。
目前常见的深度学习框架TensorFlow[28]和PyTorch[29]等都内置自动微分模块，本文主要采用TensorFlow2.3进行PINNs计算。
（3）PINNs算法求解达西流速
PINNs算法进行渗流流速求解时，有两种不同的求解策略。方法一是通过PINNs算法求解达西渗流定律和渗流模型的耦合方程，其求解原理如图 1所示，首先构造输入为x，输出为ˆH(x;θ)和ˆV(x;θ)的FNN分别对水头H和流速V做近似估计，其中，θ为神经网络中所有权重和偏置参数的集合。其次，通过AD对式（1），（3）和BC构造如下损失函数MSE：
![]()
图 1 PINNs算法求解渗流流速示意图Figure 1. Schematic diagram of PINNs algorithms for deriving Darcy velocityMSE=ωPMSEPDE+ωBMSEBC, (5) 其中，
MSEPDE=1|τp|∑x∈τp[‖ (6) \operatorname{MSE}_{\mathrm{BC}}=\frac{1}{\left|\tau_{\mathrm{B}}\right|} \sum\limits_{x \in \tau_{\mathrm{B}}}\|\mathcal{B}(\hat{H}(\boldsymbol{x} ; \boldsymbol{\theta}))\|_2^2。 (7) 式中： {\left\|\cdot\right\|_2} 为欧几里得范数， {\tau _\mathrm{P}} ， {\tau _\mathrm{B}} 分别为计算区域内部和边界区域中用于训练神经网络的残差点集合，一般可以通过在定义域内随机抽样[30]或进行指定位置设置； {\omega _{\text{p}}} ， {\omega _\mathrm{B}} 分别为修正损失项MSEPDE，MSEBC的权重。最后，通过Adam[31]、L-BFGS-B[32]等优化算法对损失函数进行寻优即可得到最优的神经网络参数 \hat{\boldsymbol{\theta}} ：
(\hat{\boldsymbol{\theta}})=\arg \min _{\boldsymbol{\theta}} \operatorname{MSE}(\boldsymbol{\theta}) 。 (8) 同理，方法二则是通过PINNs算法构造输入为x，输出为 \hat H(x;\boldsymbol{\theta} ) 的神经网络计算地下水头分布 \hat H(x;\hat {\boldsymbol{\theta}}) ，再基于达西定律，利用AD计算求解 \hat V(x;\boldsymbol{\theta} ) 。值得一提的是，由于方法二求解的控制方程不同，此时MSEPDE变为
\mathrm{MSE}_{\mathrm{PDE}}=\frac{1}{\left|\tau_{\mathrm{P}}\right|} \sum\limits_{\boldsymbol{x} \in \tau_{\mathrm{p}}}\left\|-\nabla^{\mathrm{T}}(\boldsymbol{K} \nabla \hat{H}(\boldsymbol{x} ; \boldsymbol{\theta}))-W\right\|_2^2。 (9) （4）PINNs算法识别渗流参数
PINNs算法不仅能够实现达西流速模拟，还能根据水头的实测或模拟数据对渗流参数进行反演，其原理如图 2所示。
![]()
图 2 PINNs算法反演渗流参数示意图Figure 2. Schematic diagram of PINNs algorithms for identifying seepage parameters在进行均质渗流参数反演时，基于正向水头求解的损失函数及实测点 {\boldsymbol{x}_D} 的H值（H( {\boldsymbol{x}_D} )），加入衡量实测数据与预测值间偏差的MSEData损失项及其权重 {\omega _D} ：
\omega_{\mathrm{D}} \mathrm{MSE}_{\mathrm{Data}}=\omega_{\mathrm{D}} \frac{1}{\left|\tau_{\mathrm{D}}\right|} \sum\limits_{\boldsymbol{x} \in \tau_{\mathrm{D}}}\left\|\hat{H}(\boldsymbol{x} ; \boldsymbol{\theta})-H\left(\boldsymbol{x}_{\mathrm{D}}\right)\right\|_2^2。 (10) 式中： {\tau _D} 为实测数据组成的集合。同理，通过优化算法对总损失函数进行优化可得到最终的神经网络参数 \hat {\boldsymbol{\theta}} 和均质渗流参数 \hat {\boldsymbol{K}} 。
PINNs算法进行非均质渗流参数反演时，由于K随x变化，需要构造多物理场神经网络对渗流参数进行反演，即构造两个神经网络 \hat H(x;{\boldsymbol{\theta} _{\text{1}}}) 和 \hat {\boldsymbol{K}}(x;{\boldsymbol{\theta }_{\text{2}}}) 分别对水头 H(\boldsymbol{x}) 和 K(\boldsymbol{x}) 做近似估计，此时 {\boldsymbol{\theta} _{\text{1}}} 和 {\boldsymbol{\theta} _{\text{2}}} 为各自神经网络的待优化参数，总的神经网络参数可表示为 \boldsymbol{\theta}= \{ {\boldsymbol{\theta }_{\text{1}}}, {\boldsymbol{\theta} _{\text{2}}}\} 。相应地，MSEPDE、MSEBC和MSEData分别变为
\mathrm{MSE}_{\mathrm{PDE}}=\frac{1}{\left|\tau_{\mathrm{P}}\right|} \sum \limits_{\boldsymbol{x} \in \tau_{\mathrm{P}}}\left\|-\nabla^{\mathrm{T}}\left(\hat{\boldsymbol{K}}\left(\boldsymbol{x} ; \boldsymbol{\theta}_2\right) \nabla \hat{H}\left(\boldsymbol{x} ; \boldsymbol{\theta}_1\right)\right)-W\right\|_2^2， (11) \operatorname{MSE}_{\mathrm{BC}}=\frac{1}{\left|\tau_{\mathrm{B}}\right|} \sum \limits_{\boldsymbol{x} \in \tau_{\mathrm{B}}}\left\|\mathcal{B}\left(\hat{H}\left(\boldsymbol{x} ; \boldsymbol{\theta}_1\right)\right)\right\|_2^2{，} (12) \operatorname{MSE}_{\text {Data }}=\frac{1}{\left|\tau_{\mathrm{D}}\right|} \sum \limits_{\bf{x} \in \tau_{\mathrm{D}}}\left\|\hat{H}\left(\boldsymbol{x} ; \boldsymbol{\theta}_1\right)-H\left(\boldsymbol{x}_{\mathrm{D}}\right)\right\|_2^2 。 (13) 通过优化算法对MSE优化可得到最终的参数 \hat {\boldsymbol{\theta}} ：
(\hat{\boldsymbol{\theta}})=\arg \min \limits _{\boldsymbol{\theta}_1, \boldsymbol{\theta}_2} \operatorname{MSE}\left(\boldsymbol{\theta}_1, \boldsymbol{\theta}_2\right) 。 (14) （5）硬约束PINNs算法
PINNs算法进行偏微分方程求解或参数反演时，边界的约束主要是软约束的形式，导致计算精度会受到一定的影响。
为了解决上述问题，可以引入合适的函数A(x)和Z(x)，使PINNs计算的结果
u(\boldsymbol{x} ; \boldsymbol{\theta})=A(\boldsymbol{x})+Z(\boldsymbol{x}) a^{\mathrm{L}}(\boldsymbol{x} ; \boldsymbol{\theta})。 (15) 能够自动满足边界条件[25]，式中， {a^L}(\boldsymbol{x};\boldsymbol{\theta} ) 为神经网络的输出值。以第三类边界条件为例，假定偏微分方程满足如下边界条件：
\hat{B}_j^{\mathrm{L}} f(\boldsymbol{x}) \equiv\left[\lambda_j^{\mathrm{L}} f(\boldsymbol{x})+\mu_j^{\mathrm{L}} \frac{\partial f(\boldsymbol{x})}{\partial x_j}\right]_{x_j=a_j} , (16) \hat{B}_j^{\mathrm{R}} f(\boldsymbol{x}) \equiv\left[\lambda_j^{\mathrm{R}} f(\boldsymbol{x})+\mu_j^{\mathrm{R}} \frac{\partial f(\boldsymbol{x})}{\partial x_j}\right]_{x_j=b_j} 。 (17) 式中：f(x)为偏微分方程的解，xj为x第j个空间坐标； \hat B_j^L ， \hat B_j^R 分别为xj=aj和xj=bj的左、右边界条件算子； \lambda 和μ为已知常数。A(x)和Z(x)对于第三类边界条件，可以进行如下表示[25]：
A(\boldsymbol{x}) = \left[ {1 - \prod\limits_{j = 1}^n {\left( {1 - {\mathcal{M}_j}} \right)} } \right]f(\boldsymbol{x}) \text{，} (18) Z(\boldsymbol{x}) = \prod\limits_{j = 1}^n {{{({x_j} - {a_j})}^2}} {({x_j} - {b_j})^2} \text{，} (19) 式中， {\mathcal{M}_j} 为作用在f(x)的算子，其定义为
{\mathcal{M}_j}f(\boldsymbol{x}) \equiv \mathit{\Phi} _j^{\text{R}}{\text{(}}{x_j}{\text{)}}\hat B_j^{\text{L}}f(\boldsymbol{x}) - \mathit{\Phi} _j^{\text{L}}({x_i})\hat B_j^{\text{R}}f(\boldsymbol{x}) \text{，} (20) 其中，
\left. \begin{array}{l}{\mathit{\Phi} }_{j}^\text{R}({x}_{j})=\frac{{\lambda }_{j}^\text{R}({x}_{j}-{b}_{j})-{\mu }_{j}^\text{R}}{({a}_{j}-{b}_{j}){\lambda }_{j}^{\text{L}}{\lambda }_{j}^\text{R}+{\mu }_{j}^{\text{L}}{\lambda }_{j}^\text{R}-{\lambda }_{j}^{\text{L}}{\mu }_{j}^\text{R}}\text{，}\\ {\mathit{\Phi} }_{j}^{\text{L}}({x}_{j})=\frac{{\lambda }_{j}^{\text{L}}({x}_{j}-{a}_{i})-{\mu }_{j}^{\text{L}}}{({a}_{j}-{b}_{j}){\lambda }_{j}^{\text{L}}{\lambda }_{j}^\text{R}+{\mu }_{j}^{\text{L}}{\lambda }_{j}^\text{R}-{\lambda }_{j}^{\text{L}}{\mu }_{j}^\text{R}}。\end{array}\right\} (21) 2. 数值试验
在本节中，将PINNs算法应用于二维（空间坐标为x，y）均质和非均质无量纲渗流模型（参数量纲均为1）求解中，通过与解析解和实际解进行对比，验证硬约束PINNs算法对于渗流流速的求解（以x方向的渗流流速Vx为例）和渗流参数反演的准确性。采用如下缩写方式：PINNs-H-I、PINNs-H-II分别表示方法一、二的硬约束PINNs算法，PINNs-S-I、PINNs-S-II分别表示方法一、二的软约束PINNs算法。所有的PINNs算法都使用Python编写，修正的损失项权重设置为1，优化算法采用学习率为0.001，迭代步数为20000的Adam算法和L-BFGS-B算法，即使用Adam算法迭代至指定步数后，再通过L-BFGS-B进行优化，这样能够提高计算精度和效率[18]，其他超参数(PINNs算法运行之前设置的参数值)选择如表 1所示。
表 1 算例所用的计算超参数Table 1. Hyperparameters of PINNs algorithms used in all examples of Section 2求解方法 {\tau _{\text{P}}} {\tau _{\text{B}}} {\tau _{\text{D}}} NN(H或H，Vx) NN(K) 正向求解 2000 200 — 50×3 — 反向求解 400 200 400 100×6 60×6 注：NN表示神经网络隐藏层结构。 为了对计算结果的精度进行分析，采用L2相对误差（ \varepsilon ）作为评价指标：
\varepsilon= \frac{{\parallel u({x_i}, {y_j}) - \hat u({x_i}, {y_j}){\parallel _2}}}{{\parallel u({x_i}, {y_j}){\parallel _2}}} \text{，} (22) 式中， ({x_i}, {y_j}) 为计算域内随机分布的点， u({x_i}, {y_j}) 为解析解， \hat u({x_i}, {y_j}) 为PINNs数值计算结果。
2.1 渗流参数为均质的二维稳定流
地下水二维稳定流的水头和Vx的控制方程可由式（1）和式（2）得出，令Kx=Ky=1，求解的计算区域 \mathit{\Omega} =[0, 1]×[0, 1]，控制方程的解析解H= xy(1-x)(1-y)，控制方程满足的Dirichlet条件和源汇项W均由解析解反推得到。
为了进行对比，本例分析了在不同y条件下（y=0.3, 0.5, 0.9），PINNs-H-II、PINNs-S-II、PINNs-H-I、PINNs-S-I和FEM的渗流流速求解精度，其中FEM采用计算域为1024个单元的三角形剖分方法求解H，再基于达西定律求解Vx，如图 3所示。可以看出，4种PINNs计算方式都具有较高的计算稳定性，这表明在不同分布的 {\tau _{\text{P}}} 、 {\tau _{\text{B}}} 下（随机抽样设置），利用优化算法进行神经网络参数寻优时，能够迭代至较优的计算结果；FEM的计算误差（主要来自于划分网格时的离散误差，以及达西渗透流速在节点处不连续导致的误差[2]）高于PINNs-H-II产生的误差（由近似误差、泛化误差和优化误差共同决定[23]），低于其他3种PINNs算法所产生的误差，表明PINNs-H-II可替代FEM作为精确求解地下水达西渗流流速的计算方法；同时，还注意到方法一施加硬约束对计算性能提升不明显，且计算精度低于PINNs-H-II，其原因在于：方法一的MSEPDE是通过对水头和流速控制方程同时定义的，硬约束只对水头施加，因此无法影响流速计算结果，而方法二的MSEPDE仅由水头控制方程构造，优化过程硬约束能够对使优化算法更易迭代至全局最优，并且H进行自动微分运算得到Vx的过程是数值准确的，因此优化算法对损失函数进行优化时，方法一相较于方法二会产生更大的优化误差，且施加硬约束对计算性能提升不明显。图 4即为采用PINNs-H-II求解Vx的计算结果，其中图 4（a）为Vx的求解结果，图 4（b）为求解结果与解析解的绝对误差，可以看出，采用自动微分技术能够得到连续光滑且数值准确的数值计算结果，且绝对误差最大为6.8×10-6，具有较高的计算精度。
![]()
图 3 不同位置y下，PINNs-H-II、PINNs-S-II、PINNs-H-I、PINNs-S-I和FEM所计算Vx的L2相对误差图Figure 3. Relative errors of L2 given by PINNs-H-II, PINNs-S-II, PINNs-H-I, PINNs-S-I and FEM methods under different values of y![]()
图 4 PINNs-H-II求解均质渗流流速Vx的计算结果图Figure 4. Values of Darcy velocity Vx of homogeneous seepage and absolute errors given by PINNs-H-II method为了验证PINNs算法在进行渗流参数反演的准确性，假定Kx=Ky=1是需要反演的量，并设定残差点 {\tau _{\text{P}}} = {\tau _{\text{D}}} = 400 （由x，y方向上均匀选取20个点生成），实测值则根据解析解得到。在进行计算时，在使用方法二进行正向求解所构造损失函数的基础上增加一项关于实测值的损失函数，利用优化算法对损失函数进行全局寻优，最终反演结果如图 5所示。可以看出，随着优化迭代步数的增加，PINNs算法施加硬约束后能够更快收敛至准确值，而最终20000步后软约束反演结果为0.97，需要更多的迭代步数才能收敛至准确值，这表明，通过施加硬约束，PINNs算法能够提高对均质二维稳定流的渗流参数反演精度。
![]()
图 5 均质情况下反演渗流参数随迭代步数变化情况Figure 5. Variation of homogeneous seepage parameter with iteration step numbers in backward problem2.2 渗流参数渐变的二维非均质稳定流
本例控制方程可由式（1）和式（2）得出，令Kx=Ky=\frac{{20}}{{2(1 - x)(1 - y) + 1}}，求解的计算区域 \mathit{\Omega} =[0, 1]×[0, 1]，控制方程的解析解H=xy(1-x)(1-y)，控制方程满足的Dirichlet条件和源汇项W均由解析解反推得到。
图 6展示了基于两种硬约束PINNs算法的渗流流速的求解结果，其中图 6（a），（b）分别为PINNs-H-I、PINNs-H-II与解析解之间的绝对误差，与2.1节中得到的结果类似，PINNs-H-II求解结果的绝对误差比PINNs-H-I低了约3个数量级，表明在非均质渗流条件下，PINNs-H-II的计算精度高于PINNs-H-I。
![]()
图 6 PINNs算法求解非均质渗流流速Vx与解析解的绝对误差云图Figure 6. Absolute errors of Vx of non-homogeneous seepage given by PINNs-H-I and PINNs-H-II methods为了进一步说明软硬约束在边界的作用情况，对比了在y=0下，PINNs-H-II和PINNs-S-II的x与H的计算关系，如图 7所示。可以看出，硬约束PINNs算法的水头计算结果在边界处恒等于0，这与实际水头相吻合，而软约束PINNs算法在边界计算结果存在振荡现象，因此可说明硬约束PINNs算法能够对软约束边界条件进行修正，具有更高的计算精度。
![]()
图 7 y=0时，PINNs-H-II和PINNs-S-II所计算的x与H的关系Figure 7. Relationship between x and H via PINNs-H-II and PINNs-S-II methods at y=0对于渗流参数的反演，在方法二求解水头问题的基础上，假定Kx=Ky=\frac{{20}}{{2(1 - x)(1 - y) + 1}}为待反演函数，构建多物理场神经网络模型，求解结果如图 8所示，图 8（a）为反演的渗流参数，图 8（b）为反演的渗流参数与实际解的绝对误差云图。可以看出，最大绝对误差在4.5×10-2附近，因此可证明硬约束PINNs算法反演具有较高的计算精度，并且从图 8（c）可知，在y=0.5处，硬约束PINNs算法能够准确反演出渗流参数随空间位置x的变化关系，而软约束PINNs算法却无法获得准确的求解结果，说明通过施加硬约束，能够提高反演的准确性。
![]()
图 8 非均质情况下渗流参数反演结果Figure 8. Identified seepage parameters and corresponding absolute errors via soft- and hard-constraint PINNs in non-homogeneous seepage3. 结论
传统的数值算法主要是通过划分网格的方式进行地下水模型的渗流流速计算，较难实现非均质模型的渗流参数反演，本文将PINNs算法应用于地下渗流模型的正向和反向求解当中，通过施加硬约束提高了算法的精度，主要得到以下3点结论。
（1）基于硬约束PINNs算法构造了PINNs-H-I和PINNs-H-II两种渗流流速计算方法，分析表明，两种计算方法都具有较高的求解稳定性，并通过对比可知，PINNs-H-I在进行损失函数构造时，MSEPDE需要增加求解V的损失项，导致损失函数更为复杂，因此优化算法更难收敛至最优点，而PINNs-H-II仅由渗流方程确定，进行寻优时具有更高的计算精度。
（2）对PINNs算法施加硬约束可自动满足偏微分方程的边界要求，并能够提高对渗流模型求解和渗流参数反演的精度。
（3）针对渗流参数的分布特点，采用硬约束PINNs算法分布构造了单物理场或多物理场的神经网络模型，通过最小化神经网络的损失函数，对均质或非均质渗流参数进行了反演。在算例中，通过与实际的渗流参数进行对比，验证了反演参数值的准确性。
-
![]()
图 5 不同循环应力比条件下胶结砂土动力学响应
Figure 5. Cyclic responses of cemented sand under different cyclic stress ratios
![]()
图 8 循环应力比0.425条件下胶结砂土微观力学性质
Figure 8. Micro-mechancial behaviors of cemented sand under cyclic stress ratio of 0.425
![]()
图 9 循环应力比0.5条件下胶结砂土微观力学性质
Figure 9. Micro-mechancial behaviors of cemented sand under cyclic stress ratio of 0.5
表 1 离散元模型参数
Table 1 Material parameters used in DEM simulations
微观参数 数值 颗粒部分 弹性模量/(N·m-2) 7.0×108 法向切向刚度比 5 抗转动系数 0.25 局部压碎系数 2.1 摩擦系数 0.5 胶结部分 弹性模量/(N·m-2) 3.5×107 泊松比 0.32 抗拉强度/(N·m-2) 0.5×106 抗压强度/(N·m-2) 0.5×107 临界长细比 0.05 半径乘子 0.332 
下载: 导出CSV
表 2 数值模拟方案
Table 2 Program of DEM simulation
离散元试样 循环应力比CSR 净砂 0.1, 0.125, 0.15, 0.2 胶结砂土 0.2, 0.4, 0.425, 0.45, 0.5
下载: 导出CSV
-
[1] MASSARSCH K R, FELLENIUS B H. Vibratory compaction of coarse-grained soils[J]. Canadian Geotechnical Journal, 2002, 39(3): 695-709. doi: 10.1139/t02-006
[2] WANG Y H, LEUNG S C. Characterization of cemented sand by experimental and numerical investigations[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2008, 134(7): 992-1004. doi: 10.1061/(ASCE)1090-0241(2008)134:7(992)
[3] PORCINO D, MARCIANÒ V, GRANATA R. Static and dynamic properties of a lightly cemented silicate-grouted sand[J]. Canadian Geotechnical Journal, 2012, 49(10): 1117-1133. doi: 10.1139/t2012-069
[4] DEJONG J T, SOGA K, KAVAZANJIAN E, et al. Biogeochemical processes and geotechnical applications: progress, opportunities and challenges[J]. Géotechnique, 2013, 63(4): 287-301. doi: 10.1680/geot.SIP13.P.017
[5] 刘汉龙, 肖鹏, 肖杨, 等. MICP胶结钙质砂动力特性试验研究[J]. 岩土工程学报, 2018, 40(1): 38-45. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201801003.htm LIU Han-long, XIAO Peng, XIAO Yang, et al. Dynamic behaviors of MICP-treated calcareous sand in cyclic tests[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(1): 38-45. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201801003.htm
[6] SANTAMARINA J C, KLEIN A, FAM M A. Soils and Waves[M]. New York: John Wiley and Sons, 2001: 25-54.
[7] MITCHELL J K, SOGA K. Fundamentals of Soil Behavior[M]. 3rd ed. New York: John Wiley and Sons, 2005: 5-33.
[8] FRYDMAN S, HENDRON D, HORN H, et al. Liquefaction study of cemented sand[J]. Journal of Geotechnical and Geoenvironmental Engineering, 1980, 106(GT3): 275-297.
[9] CLOUGH G W, IWABUCHI J, RAD N S, et al. Influence of cementation on liquefaction of sands[J]. Journal of Geotechnical Engineering, 1989, 115(8): 1102-1117. doi: 10.1061/(ASCE)0733-9410(1989)115:8(1102)
[10] QADIMI A, COOP M R. The undrained cyclic behavior of a carbonate sand[J]. Géotechnique, 2007, 57(9): 739-750. doi: 10.1680/geot.2007.57.9.739
[11] PORCINO D, MARCIANÒ V, GRANATA R. Cyclic liquefaction behaviour of a moderately cemented grouted sand under repeated loading[J]. Soil Dynamics and Earthquake Engineering, 2015, 79(Part A): 36-46.
[12] VRANNA A, TIKA T. Undrained monotonic and cyclic response of weakly cemented sand[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2020, 146(5): 04020018. doi: 10.1061/(ASCE)GT.1943-5606.0002246
[13] RASOULI H, FATAHI B, NIMBALKAR S. Liquefaction and post-liquefaction assessment of lightly cemented sands[J]. Canadian Geotechnical Journal, 2020, 57: 173-188. doi: 10.1139/cgj-2018-0833
[14] CUNDALL P A, STRACK O D L. Discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65. doi: 10.1680/geot.1979.29.1.47
[15] 蒋明镜. 现代土力学研究的新视野——宏微观土力学[J]. 岩土工程学报, 2019, 41(2): 195-254. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201902002.htm JIANG Ming-jing. New paradigm for modern soil mechanics: geomechanics from micro to macro[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(2): 195-254. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201902002.htm
[16] DE BONO J P, MCDOWELL G R. Discrete element modelling of one-dimensional compression of cemented sand[J]. Granular Matter, 2014, 16(1): 79-90. doi: 10.1007/s10035-013-0466-0
[17] ZHANG F G, JIANG M J. Do the normal compression lines of cemented and uncemented geomaterials run parallel or converge to each other after yielding?[J]. European Journal of Environmental and Civil Engineering, 2021, 25(2): 368-386. doi: 10.1080/19648189.2018.1531788
[18] WANG Y H, LEUNG S C. A particulate-scale investigation of cemented sand behavior[J]. Canadian Geotechnical Journal, 2008, 45(1): 29-44. doi: 10.1139/T07-070
[19] JIANG M J, ZHANG F G, THORNTON C. A simple three-dimensional distinct element modeling of the mechanical behavior of bonded sands[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39(16): 1791-1820. doi: 10.1002/nag.2371
[20] 张伏光, 蒋明镜. 胶结砂土强度特性的真三轴试验离散元分析[J]. 岩石力学与工程学报, 2018, 37(6): 1530-1539. doi: 10.13722/j.cnki.jrme.2017.1464 ZHANG Fu-guang, JIANG Ming-jing. Analysis of the strength behaviour of cemented sands in true triaxial test with distinct element method[J]. Chinese Journal of Rock Mechanics and Engineering, 2018, 37(6): 1530-1539. (in Chinese) doi: 10.13722/j.cnki.jrme.2017.1464
[21] 蒋明镜, 孙若晗, 李涛, 等. 微生物处理砂土不排水循环三轴剪切 CFD-DEM模拟[J]. 岩土工程学报, 2020, 42(1): 20-28. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202001005.htm JIANG Ming-jing, SUN Ruo-han, LI Tao, et al. CFD-DEM simulation of microbially treated sands under undrained consolidated cyclic triaxial tests[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(1): 20-28. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202001005.htm
[22] SHEN Z F, JIANG M J, THORNTON C. DEM simulation of bonded granular material: Part I contact model and application to cemented sand[J]. Computers and Geotechnics, 2016, 75: 192-209. doi: 10.1016/j.compgeo.2016.02.007
[23] JIANG M J, ZHANG A, LI T. Distinct element analysis of the microstructure evolution in granular soils under cyclic loading[J]. Granular Matter, 2019, 21: 39. doi: 10.1007/s10035-019-0892-8
[24] JIANG M J, KONRAD J, LEROUEIL S. An efficient technique for generating homogeneous specimens for DEM studies[J]. Computers and Geotechnics, 2003, 30(7): 579-597. doi: 10.1016/S0266-352X(03)00064-8
[25] JIANG M J, SUN Y G, LI LQ, et al. Contact behavior of idealized granules bonded in two different interparticle distances: an experimental investigation[J]. Mechanics of Materials, 2012, 55: 1-15. doi: 10.1016/j.mechmat.2012.07.002
[26] 申志福. 深海能源土力学特性三维多尺度数值模拟[D]. 上海: 同济大学, 2016. SHEN Zhi-fu. Three-Dimensional Multi-Scale Numerical Simulations of the Mechanical Behavior of Methane Hydrate Bearing Sediments[D]. Shanghai: Tongji University, 2016. (in Chinese)
[27] KUHN M R, RENKEN H E, MIXSELL A D, et al. Investigation of cyclic liquefaction with discrete element simulations[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2014, 140(12): 04014075. doi: 10.1061/(ASCE)GT.1943-5606.0001181
[28] THORNTON C. Numerical simulations of deviatoric shear deformation of granular media[J]. Géotechnique, 2000, 50(1): 43-53. doi: 10.1680/geot.2000.50.1.43
[29] BOLTON M D, NAKATA Y, CHENG Y P. Micro- and macro-mechanical behaviour of DEM crushable materials[J]. Géotechnique, 2008, 58(6): 471-480. doi: 10.1680/geot.2008.58.6.471
[30] SAZZAD M M, SUZUKI K. Density dependent macro-micro behavior of granular materials in general triaxial loading for varying intermediate principal stress using DEM[J]. Granular Matter, 2013, 15(5): 583-593. doi: 10.1007/s10035-013-0422-z
[31] MARTIN E L, THORNTON C, UTILI S. Micromechanical investigation of liquefaction of granular media by cyclic 3D DEM tests[J]. Géotechnique, 2020, 70(10): 906-915. doi: 10.1680/jgeot.18.P.267
-
期刊类型引用(11)
1. 李淑娥，陈志明，徐永福，徐宇冉，康峰沂，杜仲宝. 基于颗粒分布分形模型毛细水上升高度计算分析. 岩土工程学报. 2024(10): 2221-2228 . 
本站查看
2. 曲诗章，刘晓明，黎莉，陈仁朋. 基于双分形级配模型参数的粗粒土渗透系数计算公式. 岩土工程学报. 2023(01): 144-152 . 本站查看
3. 韩志洋，曹志翔，黄开放. 基于离散元模拟的土石混合体剪切与变形特性研究. 中国农村水利水电. 2023(05): 238-244 . 百度学术
4. 刘晓义，胡敏，刘大顺. 基于离散元法的砂砾石颗粒破碎特征研究. 低温建筑技术. 2023(12): 24-28 . 百度学术
5. 孟敏强，肖杨，孙增春，张志超，蒋翔，刘汉龙，何想，吴焕然，史金权. 粗粒料及粒间微生物胶结的破碎-强度-能量耗散研究进展. 中国科学:技术科学. 2022(07): 999-1021 . 百度学术
6. 王瑞，郭聚坤，尹斌，雷胜友，魏道凯. 钙质砂颗粒形状及破碎特性试验研究. 海洋工程. 2022(05): 158-166 . 百度学术
7. 陈晓斌，郭云鹏，蔡德钩，尧俊凯，肖源杰. 铁路工程粗颗粒土路基填料研究现状与发展综述. 路基工程. 2021(03): 1-11 . 百度学术
8. 叶阳升，朱宏伟，尧俊凯，蔡德钩，安再展. 高速铁路路基振动压实理论与智能压实技术综述. 中国铁道科学. 2021(05): 1-11 . 百度学术
9. 于玉贞，张向韬，王远，吕禾，孙逊. 堆石料真三轴条件下力学特性试验研究进展. 工程力学. 2020(04): 1-21+29 . 百度学术
10. 王晓帅，王子寒，景晓昆，肖成志. 粗粒土大型直剪试验宏细观研究与离散元模拟. 深圳大学学报(理工版). 2020(03): 279-286 . 百度学术
11. 孟敏强，王磊，蒋翔，汪成贵，刘汉龙，肖杨. 基于尺寸效应的粗粒土单颗粒破碎试验及数值模拟. 岩土力学. 2020(09): 2953-2962 . 百度学术
其他类型引用(20)
计量
- 文章访问数: 430
- HTML全文浏览量: 31
- PDF下载量: 362
- 被引次数: 31
