# Turbulent-Fluid-Based Simulation of Dynamic Liquefaction Using Large Deformation Analysis of Solid Phase

Rohollah Taslimian

College of Engineering and Science, Louisiana Tech University, Ruston, United States

## Article history

Received: 14-12-2023

Revised: 17-01-2024

Accepted: 23-01-2024

Emails: rta017@latech.edu; r.taslimian@gmail.com

**Abstract:** To accurately simulate the cyclic liquefaction phenomenon, it is necessary to establish a well-defined coupling between the solid particles and the high-velocity percolating pore fluid. For this purpose, the existing literature extensively examines large deformation liquefaction analyses for the solid phase, regardless of the fluid flow nonlinearity. On the other hand, the transient fluid pressure gradient governed by the turbulence flow has been previously addressed using the infinitesimal strain tensor for the solid phase. However, an integration of these two aspects has not been considered yet. Therefore, the nonlinear terms of not only the strain tensor but also the momentum balance of the fluid should be taken into account simultaneously. This study introduces these critical inherent nonlinearities related to the material softening caused by pore pressure buildup during seismic events as an open area of research.

**Keywords:** Large Deformation, Non-Darcy Flow Model, Soil Liquefaction, Soil-Pore Fluid Interaction, Wave Propagation, Numerical Modelling

## Introduction

Liquefaction is a general phenomenon where a mixture of fluid and granular material subjected to vibration ultimately leads to circulation patterns similar to both fluid convection and granular convection. Indeed, liquefaction is fluid-granular convection with circulation patterns which are known as sand boils or sand volcanoes in the study of soil liquefaction. Such nonlinear behavior requires an appropriate model to describe.

When saturated porous media experiences ground motion, it induces movement of pore fluid inside the solid skeleton because of soil settlement. In fact, during an excitation, the excess pore pressure gradient derives a vertical fluid flow inside saturated soil towards the free surface. The temporal displacement of pore fluid can notably result in the redistribution of pore fluid pressure, influenced by factors such as soil permeability, loading rate, boundary conditions and pressure gradient. According to this note, liquefaction simulation approaches should explicitly consider the significance of nonlinear fluid flow and large deformation of solid skeleton at the same time. The current state-of-the-art in the simulation of liquefaction aimed at high seepage velocity of the fluid as well as large deformation of the solid skeleton has been reviewed in this study.

## Governing Formulations

### Fluid Flow Laws in Saturated Porous Media

The equations for Darcy (linear) and non-Darcy (nonlinear) flow laws are provided below.

### Darcy Flow Law

Fluid flow in porous media can be characterized by the linear Darcy regime when the Reynolds number is adequately low (Chai et al., 2010), a condition that can be expressed as follows:

$$-\nabla p = \frac{\rho_f g}{k} \dot{w} \qquad (1)$$

$\nabla p$ is the fluid flow pressure gradient, $\rho_f$ is the density of the pore fluid, $\dot{w}$ is the average relative velocity of the fluid, $g$ is the gravity acceleration and $k$ is the permeability coefficient. In this linear regime, the viscous forces dominate over the inertia forces, but only the pore structure affects the fluid flow behavior.

### Non-Darcy Flow Law

By increasing the Reynolds number, the fluid flow behavior is determined by nonlinear laws such as the quadratic Forchheimer regime (Scheidegger, 1957).
$$-\nabla p = \frac{\rho_f g}{k} \dot{w} + S \rho_f \dot{w}^2 \quad (2)$$

S is the non-darcy constant and it is related to the shape, pore size and porosity. In the non-darcy regime, the impact of inertia is more prominent than viscous effects, resulting in the governing influence of inertia forces on fluid flow behavior.

## Fully Coupled Formulations for Saturated Porous Media

The fundamental laws and formulations to model coupled nonlinear fluid flow and large solid skeleton deformation within continuous and fully saturated granular particles have been presented in this section. The subsequent four expressions delineate the equilibrium of momentum for the solid-fluid mixture, the momentum equilibrium of the fluid, effective stress in the soil mass and the conservation of fluid mass (Taslimian et al., 2012; 2015; 2023).

### Momentum Balance of the Mixture

$$\sigma_{ij,j} + \rho b_i = \rho \ddot{u}_i + \rho_f [\ddot{w}_i + \dot{w}_j \dot{w}_{i,j}] \quad (3)$$

where, $\sigma \rho_{ij,j}$ is the divergence of the total stress tensor, and $\ddot{u}_i$ and $\dot{w}_i$ respectively are the solid skeleton and average relative fluid accelerations. $b_i$ is the body acceleration. $\rho$ and $\rho_f$ are the density of the total mixture and the pore fluid. $\dot{w}_j \dot{w}_{i,j}$ the convective terms of relative fluid acceleration.

### Momentum Balance of the Fluid

$$-p_i + \rho_f b_i = \rho_f \ddot{u}_i + \rho_f [\ddot{w}_i + \dot{w}_j \dot{w}_{i,j}] / n + \frac{\rho_f g}{k_{ji}} \dot{w}_j + S \rho_f \dot{w}_i^2 \quad (4)$$

where, $n$ is porosity, $k_{ji}$ is the permeability tensor equal to $k \delta_{ji}$ in the isotropic case, $\delta_{ji}$ is the Kronecker delta, $w_i$ is the average relative fluid velocity and $p_i$ is the pore fluid pressure gradient, $\frac{\rho_f g}{k_{ji}} \dot{w}_j$ the viscous drag forces and $S \rho_f \dot{w}_i^2$ the inertia forces.

### Effective Stress

$$\sigma'_{ij} = \sigma_{ij} + \alpha \delta_{ij} p \quad (5)$$

where, $\alpha$ is the Biot alpha and $\sigma'_{ij}$ is the modified effective stress. If the soil and fluid particles are incompressible then $\alpha = 1$.

### Conservation of Fluid Mass

$$\dot{p}/Q + \alpha \dot{E}_{ii} + \dot{w}_{i,i} + \frac{n \rho_f}{\rho_f} + \dot{s}_0 = 0 \quad (6)$$

$$\frac{1}{Q} = \frac{n}{\kappa_f} + \frac{\alpha - n}{\kappa_s} \quad (7)$$

where, $k_f$ and $k_s$ represent the fluid bulk modulus and bulk modulus of the soil grains, respectively. The strain tensor is formulated by either Lagrangian or Eulerian tensors as follows:

$$E_{IJ} = 0.5 (u_{i,j} + u_{j,i} + u_{K,i} u_{K,j}) \quad (8) \\ \text{(Lagrangian strain tensor)}$$

$$e_{ij} = 0.5 (u_{i,j} + u_{j,i} + u_{k,i} u_{k,j}) \quad (\text{Eulerian strain tensor}) \quad (9)$$

### Constitutive Model

The general elastoplastic behavior can be characterized by the correlation between increments in stress and strain, as outlined by Zienkiewicz et al. (1999):

$$d\sigma'_{ij} = D_{ijkl} (d\varepsilon_{kl} - d\varepsilon_{kl}^{\theta}) + \sigma'_{ik} dw_{kj} + \sigma'_{jk} dw_{ki} \quad (10)$$

Here, $D_{ijkl}$ denotes the tangent matrix of the solid skeleton, contingent on the current state and history of stresses and strains, as well as the direction of strain increments. $d\varepsilon_{kl}^{\theta}$ is the increment of the thermal strain or similar autogenous strains. The last two terms account for the simple rotation of existing stress components. Liquefiable sands can be modeled utilizing a generalized plasticity composed of a yield surface along with a non-associated flow rule called the Pastor-Zienkiewicz Mark III model (Pastor et al., 1990) under dynamic loads.

## Numerical Solutions

To numerically solve the governing equations, the finite element method and the generalized Newmark method are primarily utilized for spatial and time discretization respectively (Zienkiewicz and Taylor, 1991; Katona and Zienkiewicz, 1985). To solve large deformation problems in continuum mechanics, Lagrangian finite element methods derived from Cauchy stresses introduced by improving the small strain tensor to the Total Lagrangian (TL) and Updated Lagrangian (UL) analyses, assuming that the original geometry is unchanged (Bathe, 1996). The equations of TL and UL include the second Piola-Kirchhoff stress tensor and the Lagrangian strain tensor (Eq. 8). In these analyses, the nonlinear term of the Lagrangian strain tensor is added to the infinitesimal strain tensor leading to changes in global stiffness as geometry changes (geometrical nonlinearity). To predict the behavior of saturated soils under large strain conditions, a mixed finite element and finite difference methods based on the UL method are presented (Di and Sato, 2004). Although the finite element method has been employed widely, its inherent defects cannot be neglected easily because the severe mesh distortions caused by large displacement, large rotation, or both, may result in interruption of the calculations suddenly or loss of fidelity. To tackle this problem, some methods including decoupled Arbitrary Lagrangian-Eulerian (ALE)
and coupled Eulerian-Lagrangian (multi-material ALE) methods have recently gained popularity in computational geomechanics (Liu et al., 2022; Bakroon et al., 2020) to model the possible mesh distortion in the large deformation analysis of the liquefaction. The Excess Pore Water Pressure (EPWP) under the ALE method which considers the mesh distortion and the conventional UL method are shown in Fig. 1. As displayed in the picture, when the deformation is small (before 10 s), the EPWP in the two methods have no obvious dissimilarity. However, the model liquefies after 10 seconds and large deformation appears, causing more difference between the two models and more liquefaction potential when using the ALE method. In fact, the average error in the solutions increases using the UL method after 10 s. Therefore, by preserving a finite mesh under healthy conditions when modeling large deformations, liquefaction potential increases.

## Evaluation of Liquefaction Potential

Large post-liquefaction shear distortion of granular material and its progressive development are well-known (Chiaro et al., 2013), but a generally accepted simulation for this phenomenon still needs more investigation. In a study, typical stress and strain results have been presented (Fig. 2b) for pre-liquefaction and post-liquefaction. The results of the study illustrate the cumulative decrease of effective stress toward post-liquefaction and the “butterfly orbit” of the stress path with a large deformation (Fig. 2a). In addition, the stress-strain curves in Fig. 1b show the occurrence of large but bounded shear strains whenever the stress state triggers liquefaction (Wang et al., 2016). In another effort, two large and small strain approaches have been employed to predict the liquefaction potential. The results of the predictions reveal that liquefaction takes place earlier in seabed deposits and at greater depths in the soil surrounding an embankment during strong earthquake motion using the large deformation analysis in comparison to the small strain theory (Di and Sato, 2003).

Additionally, the fluid flow law in porous media is nonlinear for an adequately high seepage velocity, particularly when the permeability coefficient is high enough and when an external water pressure gradient is applied to the porous media. To examine the effect of the nonlinear fluid flow on the liquefaction potential of sand, the stress path and the cyclic shear stress-strain relationship during the cyclic loading are depicted in Figs. 3a-b to compare linear (Darcy) flow and nonlinear (non-darcy) flow. The stress paths as well as the cyclic shear stress-strain plots show a more liquefaction potential for nonlinear flow in comparison to linear flow (Taslimian et al., 2012; 2015; 2023).

<figure><img src="image_1.png" /><figcaption>Fig. 1: Excess pore water pressure of an Element (Liu et al., 2022), Under the terms of the Creative Commons CC BY-NC-ND license, copyright 2022, Soils and Foundations published by Elsevier</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Fig. 2: A comparison between large deformation and small deformation liquefaction; (a) The stress path response; (b) The cyclic shear stress- shear strain curve (Wang et al., 2016), Under the terms of the Creative Commons CC BY license, copyright 2016, Acta Geotechnica published by Springerlink.com</figcaption></figure>

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3:** A comparison between nonlinear fluid flow and linear fluid flow in liquefaction simulation; (a) The stress path response; (b) The cyclic shear stress-strain curve (Taslimian, 2012)</figcaption></figure>

## Conclusion

The results of the studies demonstrate that the liquefaction behavior of saturated soils under dynamic loading, to a great extent, depends on the nonlinear terms of the strain tensor and the inertia forces in the fluid momentum balance formulation. The nonlinear fluid flow analyses indicate a greater liquefaction potential in comparison to linear fluid flow models. Similarly, the granular material tends to liquefy more using the large deformation analyses compared to the small strain analyses. Thus, in order to accurately evaluate liquefaction potential, these two important factors should not be considered separately because both of them cause an increase in the liquefaction potential. Achieving such a high level of accuracy in modeling porous media, without expensive computational costs, may require a novel numerical method.

## References

Bakroon, M., Daryaei, R., Aubram, D., & Rackwitz, F. (2020). Investigation of mesh improvement in multimaterial ALE formulations using geotechnical benchmark problems. *International Journal of Geomechanics*, 20(8), 04020114.
https://doi.org/10.1061/(ASCE)GM.1943-5622.0001723

Bathe, K. J. (1996). Finite element procedures prentice-hall. *New Jersey*, 1037(1).
https://lccn.loc.gov/95001231

Chai, Z., Shi, B., Lu, J., & Guo, Z. (2010). Non-Darcy flow in disordered porous media: A lattice Boltzmann study. *Computers and Fluids*, 39(10), 2069-2077.
https://doi.org/10.1016/j.compfluid.2010.07.012

Chiaro, G., Kiyota, T., & Koseki, J. (2013). Strain localization characteristics of loose saturated Toyoura sand in undrained cyclic torsional shear tests with initial static shear. *Soils and Foundations*, 53(1), 23-34.
https://doi.org/10.1016/j.sandf.2012.07.016

Di, Y., & Sato, T. (2003). Liquefaction analysis of saturated soils taking into account variation in porosity and permeability with large deformation. *Computers and Geotechnics*, 30(7), 623-635.
https://doi.org/10.1016/S0266-352X(03)00060-0

Di, Y., & Sato, T. (2004). A practical numerical method for large strain liquefaction analysis of saturated soils. *Soil Dynamics and Earthquake Engineering*, 24(3), 251-260.
https://doi.org/10.1016/j.soildyn.2003.11.004

Katona, M. C., & Zienkiewicz, O. C. (1985). A unified set of single step algorithms part 3: The beta-m method, a generalization of the Newmark scheme. *International Journal for Numerical Methods in Engineering*, 21(7), 1345-1359.
https://doi.org/10.1002/nme.1620210713

Liu, S., Tang, X., & Li, J. (2022). A decoupled Arbitrary Lagrangian-Eulerian method for large deformation analysis of saturated sand. *Soils and Foundations*, 62(2), 101110.
https://doi.org/10.1016/j.sandf.2022.101110

Pastor, M., Zienkiewicz, O. C., & Chan, A. (1990). Generalized plasticity and the modelling of soil behaviour. *International Journal for Numerical and Analytical Methods in Geomechanics*, 14(3), 151-190.
https://doi.org/10.1002/nag.1610140302

Scheidegger, A. E. (1957). *The physics of flow through porous media*. University of Toronto press.
https://doi.org/10.3138/9781487583750

Taslimian, R. (2012). Numerical simulation of liquefaction in saturated porous media with nonlinear flow considering irregular layers. M. Sc. Thesis, *University of Tehran*, Iran.
https://doi.org/10.13140/RG.2.2.25887.33442

Taslimian, R., Noorzad, A., & Maleki Javan, M. R. (2015). Numerical simulation of liquefaction in porous media using nonlinear fluid flow law. *International Journal for Numerical and Analytical Methods in Geomechanics*, 39(3), 229-250.
https://doi.org/10.1002/nag.2297

Taslimian, R., Noorzad, A., & Maleki Javan, M. R. (2023). Numerical analysis of liquefaction phenomenon considering irregular topographic interfaces between porous layers. *Journal of Earthquake Engineering*, 27(5), 1095-1109.
https://doi.org/10.1080/13632469.2022.2038727
Taslimian, R., Noorzad, A., & Noorzad, A. (2012). Modeling saturated porous media with elasto-plastic behavior and non-Darcy flow law considering different permeability coefficients. *In Proceedings of the 15th World Conference on Earthquake Engineering*, 24-28.
https://hal.science/hal-04362537

Wang, R., Fu, P., Zhang, J. M., & Dafalias, Y. F. (2016). DEM study of fabric features governing undrained post-liquefaction shear deformation of sand. *Acta Geotechnica*, 11, 1321-1337.
https://doi.org/10.1007/s11440-016-0499-8

Zienkiewicz, O. C., Chan, A. H. C., Pastor, M., Schrefler, B. A., & Shiomi, T. (1999). *Computational Geomechanics*, 613. Chichester: Wiley.
https://doi.org/10.1002/1096-9845(200008)29:83.0.CO;2-3

Zienkiewicz, O. C., & Taylor, R. L. (1991). The finite element method-Vol 2: Solid and Fluid Mechanics. *Dynamics and Non-Linearity*, Fourth ed., London.