Numerical methods have become an indispensable tool in modern geotechnical engineering, providing powerful capabilities to analyze complex soil-structure interaction problems that are often intractable using traditional closed-form analytical solutions. By discretizing a continuum into a finite set of elements, techniques like the Finite Element Method (FEM), Material Point Method (MPM), Hybrid Element Particle Method (HEPM), Finite Difference Method (FDM), Peridynamics method (PD), and Finite Discrete Element Method (FDEM) enable engineers to model the inherently nonlinear, inelastic, and heterogeneous behavior of soil and rock under various loading and environmental conditions. These methods allow for the simulation of intricate processes such as excavation, consolidation, seepage, dynamic liquefaction, and progressive failure, facilitating a more realistic prediction of deformations, stresses, and stability. Consequently, the application of numerical modeling leads to more informed, efficient, and safe geotechnical design, optimizing solutions while providing a deeper understanding of potential failure mechanisms.
In geotechnical engineering, the deformation behavior of soil can range from the subtle strains experienced during design to the extensive displacements observed under extreme conditions. To address these challenges, we have developed an enhanced Smoothed Particle Finite Element Method (Smoothed-PFEM) framework that includes three variants: ES-PFEM, NS-PFEM, and, most notably, SNS-PFEM. SNS-PFEM eliminates the need for variable mapping, significantly improving both computational efficiency and simulation accuracy. This advanced framework further extends its capabilities to handle contact problems, soil-water coupling, and thermo-hydro-mechanical interactions, thereby offering a robust and comprehensive toolset for tackling the complex multi-physics problems encountered in modern geotechnical engineering.
1.1 Development of Edge-based smoothed PFEM
We proposed a novel edge-base strain-smoothed particle finite element method (ES-PFEM) that incorporates an effective edge-based strain smoothing technique within the original PFEM framework, effectively mitigating volumetric locking even with lower-order elements. This method exhibits super-convergence and higher accuracy compared to traditional node-based approaches, and its explicit time integration with adaptive time stepping ensures reliable computational stability.
1.2 Development of implicit Edge-based smoothed PFEM for hydromechanical problem
We developed a fully hydro-mechanically coupled edge-based strain smoothed particle finite element method (ES-PFEM) for modeling the delayed failure of cut slopes in soft sensitive clays. This approach integrates a low-order edge-based strain smoothing element for large deformation analysis, an efficient stabilization method for pore water pressure, and an elastoviscoplastic model that accounts for time dependency, anisotropy, and destructuration. Validation through consolidation tests, strip footing simulations, and slope failure case studies demonstrates its capability to accurately predict failure timing under varying slope angles and clay sensitivities.
1.3 Development of stable node-based smoothed PFEM (SNS-PFEM)
We developed a stable nodal integration particle finite element method (SNS-PFEM) that overcomes the “overly soft” problem of conventional NS-PFEM for large deformation geotechnical problems, utilizing a novel stabilization term and strain smoothing over 3-node triangular elements. Additionally, we extended this method to tackle large deformation hydromechanical problems by incorporating soil–water coupling based on Biot’s mixture theory and an innovative water pressure stabilization scheme. Extensive benchmark tests, including static and dynamic consolidation and slope failure scenarios, confirm the robustness and accuracy of our approach in geotechnical engineering applications.
1.4 Development of SNS-PFEM for incompressible fluid simulation
We developed a stabilized NS-PFEM-FIC formulation that mitigates spatial and temporal instabilities in simulating incompressible free-surface flows by employing a gradient strain field, finite increment calculus stabilization, and a novel slip boundary predictor–corrector algorithm. Validated against classical benchmark tests, this approach effectively handles pressure and velocity updates, demonstrating high accuracy compared to experimental and numerical results.
1.5 Development of Two-Layer Two-Point SNS-PFEM
We developed an improved two-phase, two-layer SNS-PFEM that overcomes previous limitations by using two overlapping mesh layers—one representing the water phase (solved as an incompressible Newtonian fluid via PFEM-T3) and the other representing the soil phase (modeled using an elastoplastic constitutive law with a modified SNS-PFEM). A fractional step algorithm is employed to efficiently and accurately couple these phases, and extensive tests on both quasi-static and large deformation scenarios confirm the method’s robustness and versatility in simulating complex water-soil interactions in geotechnical engineering.
1.6 Development of semi-implicit hybrid Eulerian-Lagrangian PFEM (HELPFEM)
We developed a novel semi-implicit hybrid Eulerian-Lagrangian PFEM (HELPFEM) that partitions the computational domain into subdomains using a refined Lagrangian approach for free-surface regions and a coarser Eulerian method elsewhere to overcome the high computational costs and remeshing issues of traditional PFEM. By integrating the Characteristic-Based Split method, Stabilized Pressure Gradient Projection technique, and an efficient volume control algorithm within a nodal-integration framework, our method delivers enhanced stability and computational efficiency for simulating incompressible free-surface and multiphase flow problems.
1.7 Implicit bES-FEM Formulation for THM Analysis in Saturated Porous Media
In the proposed framework, a cubic bubble function is first incorporated into ES-FEM to efficiently mitigate pressure oscillations that breach the Inf-Sup condition, and then the Streamline Upwind Petrov-Galerkin (SUPG) scheme is adopted in bES-FEM to effectively reduce the spurious oscillations in convection-dominated heat transfer scenarios. Moreover, the simulations of open-loop ground source energy systems demonstrate the proposed method’s exceptional capability in tackling complex THM challenges in real-world applications.
1.8 Enhanced THM coupling for anisotropic geomaterials and smoothed-FEM simulation
An SNS-FEM THM coupling simulator is developed based on an improved thermo-hydro-mechanical coupling theory. The numerical results are useful in shedding light onto the physical meaning and rigorous derivation of the governing conservation equations, in particular, how to incorporate anisotropy and how to derive an analogous term. By ignoring this analogous term and dissipation, the temperature evolution equation downgrades to its standard form used in many textbooks. Numerical examples are provided for elucidating the role of mechanical anisotropy and cross-coupling term through the THM biaxial test, as well as the role of temperature on pipeline penetration resistance.
1.9 Two-phase two-layer SNS-PFEM for THM-Coupling
A multi-physics two-layer SNS-PFEM involving THM analysis in granular mass-water interaction is developed. The governing equations for the water and granular mass were independently solved on their respective Lagrangian meshes. In the overlapping regions representing the mixture, the interaction between water and soil was modeled by drag forces, while heat exchange was handled using the Robin boundary condition. A series of benchmark tests were conducted to verify the accuracy of the multi-physics two-layer SNS-PFEM, covering cases ranging from simple two-physics coupling to complex three-physics coupling problems.
1.10 HM coupled contact method within the SNS-PFEM framework
A novel numerical algorithm is developed to solve frictional contact problems in two-phase porous media, accounting for the large deformation and interface fluid flow, within the framework of the Stable Node-based Smoothed Particle Finite Element Method (SNS-PFEM). The proposed approach builds on our previous work that was only designed for single-phase large deformation contact problems. Within this work, the hydro-mechanical coupling system is discretized using equal-order interpolation, with a polynomial pressure projection method to eliminate the pore pressure instability. Different from existing interface fluid elements, the proposed element is constructed based on the contact algorithm, thereby enabling it to handle non-conforming interfaces and tackle challenges associated with large sliding problems.
Geotechnical and marine geotechnical environments often present a wide range of deformation behaviors—from the small strains encountered during design to the large displacements observed in failure or extreme events. The Material Point Method (MPM) offers an excellent solution for these challenges due to its hybrid Eulerian–Lagrangian framework. This approach allows MPM to accurately capture small deformations while maintaining robustness against mesh distortions during large-scale deformations. Moreover, its particle-based formulation is particularly well-suited for simulating multi-phase and multi-physics coupling problems that are common in marine geotechnical scenarios, such as soil–structure interactions, fluid-structure interaction, porous media flow, and sediment transport.
2.1 Improvements of Explicit MPM and coupling scheme with FEM
We proposed an improved explicit penalty formulation based on particle positions to accurately impose Dirichlet boundary conditions and tie-contact constraints in both MPM and FEM-MPM coupling. By introducing boundary and tied reference points and using a dimensionless penalty factor scheme, our method ensures energy conservation, numerical stability, and compatibility with complex geometries.
2.2 The development of Mortar contact in MPM
Contact handling in MPM has long been challenging. To address this, we propose a novel mortar segment-to-segment frictional contact approach. We introduce boundary vertices and a new kinematic update scheme for accurate boundary representation. A weak form of contact constraints is formulated using the mortar method to enable stable segment-to-segment detection. To enforce non-penetration, we adopt the energetic barrier method within the MPM framework.
2.3 The development of two-Layer two-point MPM
We developed a high-performance semi-implicit two-phase, two-layer MPM framework based on mixture theory to model granular mass-water dynamic interactions. Our approach combines an explicit solution for the solid phase on a co-located grid with a fractional-step method for the fluid phase on a MAC staggered grid, leveraging efficient GPU programming with Taichi, advanced solvers, and optimized contact algorithms to ensure computational efficiency and robustness.
2.4 The application of MPM in Fluid-Structure-Soil Interaction
We applied the Material Point Method (MPM) to Fluid-Structure-Soil Interaction (FSSI) studies, enabling the detailed analysis of complex phenomena. This approach allows us to investigate non-linear wave generation, simulate porous media flows, and capture the intricate non-linear interactions among fluid forces, structural responses, and mooring dynamics. By leveraging MPM in the FSSI framework, we can accurately model and predict the large deformations and transient behaviors inherent in these coupled systems.
Large deformation analysis has long been a challenging topic in computational mechanics. Traditional mesh-based approaches frequently encounter issues with element distortion, while particle methods often suffer from low computational accuracy and efficiency. To address these challenges, we developed a novel numerical framework, termed the Hybrid Element Particle Method (HEPM), which combines the strengths of both mesh-based and particle approaches.
3.1 HEPM for Large Deformation Analysis in Solid Mechanics
HEPM discretizes the computational domain into two complementary components: a set of elements and a set of particles. The elements serve as auxiliary structures that can be regenerated freely to prevent mesh distortion, while the particles are used for storing history variables and establishing equilibrium equations. Physical fields are first discretized using FEM technique, then transferred to the particle level via kernel reconstruction. All subsequent computations are performed on the particles, with elements serving only as transient supports.
3.2 HEPM for Computational Fluid Dynamics
HEPM has been further extended to computational fluid dynamics by discretizing the Navier–Stokes equations through its unique particle interpolation (which is reconstructed from element shape functions). Using the updated Lagrangian formulation, HEPM naturally avoids complex treatments of convective terms and free-surface tracking required in Eulerian approaches. Moreover, HEPM can ensure strict volume/mass conservation, allowing particles to be added or removed dynamically during the computation (even in the Newton–Raphson iteration), a feature uncommon in particle methods.
Relative to its solid mechanics formulation, the present fluid dynamics formulation incorporates the following modifications:
1. Incompressibility condition added to governing equations.
2. Finite calculus (FIC) method used for pressure stabilization.
3. LU decomposition–based fractional step method employed to solve non-positive definite systems.
Efficient and reliable prediction of large-strain consolidation in soft soils is essential for ground improvement, long-term settlement control, and environmental remediation. To address these engineering challenges, we have advanced finite difference method (FDM)-based modelling that captures the key features of consolidation under complex hydro-mechanical-chemical conditions. First, a one-dimensional formulation considering self-weight was developed for large-strain consolidation analysis, which provides a fundamental basis for interpreting settlement and pore pressure dissipation. Building on this, a quasi-2D FDM scheme was established for long-term consolidation with Prefabricated Vertical Drains (PVDs) and Prefabricated Horizontal Drains (PHDs), efficiently accounting for radial-vertical and horizontal-vertical seepage interaction in practical ground improvement projects. Finally, the modelling framework has been extended to include chemical effects, such as chemico-osmotic loading and contaminant migration, enabling predictive analysis for soil improvement under coupled mechanical and environmental influences.
Large deformation analysis has long been a chaDeformation and failure in hyperelastic materials exhibit distinct characteristics not observed in purely elastic materials. The presence of pre-existing cracks and pores further complicates their mechanical behavior. Moreover, in viscoelastic materials, time-dependent deformation and damage introduce additional complexity to the failure process. To address these challenges, we developed an implicit hyperelastic non-ordinary state-based peridynamics (NOSB-PD) framework for modelling finite deformation and damage. Based on this method, we investigated the effects of crack patterns, loading conditions, and material properties on damage evolution.
Dynamic damage in solids exhibits unique characteristics that are absent in static analysis. The presence of pre-existing cracks further increases the complications. To study this phenomenon, we developed an impact damage model based on an improved state-based peridynamic model, focusing on dynamic simulation in cracked solids during impact. Based on the proposed model, we explored the mechanisms by which material properties and crack patterns affect damage. This study offers insights into damage evolution during impact and provides guidelines for assessing dynamic damage in structures with pre-existing defects.
A comprehensive bond-slip model is also introduced into the FDEM, enabling accurate reproduction of the interaction behaviour between the concrete matrix and steel fibre. Owing to its merits in quantitative analysis, the FDEM can serve as a feasible potential tool for mesoscale modelling analysis and steel fibre-reinforced concrete (SFRC) design.
Uncertainty presents a persistent challenge in geotechnical engineering, especially in cold regions where soil heterogeneity, freeze–thaw cycles, and construction deviations interact in complex ways.
This study develops a multi-scale modeling framework that explicitly accounts for uncertainty across data, parameters, and processes. The approach integrates statistical characterization of in-situ soil properties, random field theory for spatial variability, and parameter representations that consider both natural and construction-induced uncertainties, such as drilling deviation.
By simulating the thermo–hydro-mechanical (THM) behavior of frozen soils under varied conditions, the framework reveals how deviations in design or site conditions can significantly alter the evolution of temperature fields and performance of frozen barriers. The results further inform risk evaluation through seepage analysis and system reliability assessment, enabling more resilient and informed design decisions. This framework is broadly applicable to geotechnical problems and offers a generalizable strategy for uncertainty quantification in complex subsurface systems.
8.1 Vortices Generation Mechanism / Islands and headlands
The complex coastlines and the presence of numerous islands with different sizes were responsible for the generation of hundreds of coherent macro-vortices that dominated the coastal hydrodynamics. Tidal currents played a dominant role in most of the Pearl River Estuary and the Hong Kong waters.
8.2 Vortices Generation Mechanism / Compound channels
We report the results of a numerical investigation of the flow structure and mechanism of macro-vortex generation in estuarine compound channels. The main mechanism of the formation of the macro-vortices was explained in terms of vorticity production related to strong gradients of the flow depth.
8.3 Dredging Risk Assessment Strategy / K-means clustering
We formulated a new risk assessment framework for dredging operations in coastal areas based on different hazard indexes that incorporates the variability of the dredging conditions (methodology, working load and duty cycles) and the climatology the local hydrodynamics.
8.4 Dredging Risk Assessment Strategy / Marine animal protection
This study aims to apply a new risk assessment approach based on computational fluid dynamics (CFD) and unsupervised machine learning technique for quantitatively estimating the risk of high turbidity to marine mammals. As an extension, the risk to the natural habitat in PRE was statistically compared among five possible dredging locations according to HK reclamation projects.
