Research Theme
The objective of our group's research is to develop theoretical and numerical models that aid more accurate predictions for engineering design and applications. We act as a bridge between mathematical science and engineering industry. As a result, we emphasize on both the mathematical foundations and the practicality when deriving models for specific engineering problems. The figure below summarizes some of the current and previous work of this research group. Examples of previous and ongoing research conducted by the research group are listed below.
Datadriven discovery of closure for computational geomechanics models with incomplete knowledge
In this research, our goal is to three major machine learning techniques (supervised, unsupervised and reinforcement learning) to discovery new physics of pathdependent materials that cannot be otherwise determined from handcrafted models. For a given set of data, we seek to find way to gain knowledge not necessarily from excess curvefitting from training materials, but from finding the right relationships among physical quantities.
In a nutshell, the key departure of our work from the other ML applications for constitutive laws and computational mechanics is that we are not attempting to write a surrogate model or constitutive laws from data. Instead, we are attempting to modeling the actions and behaviors of a mechanician who is tasked with writing a constitutive laws for a given set of data. As such, we introduce a game whose objective is to write a predictive constitutive law and then recast a modeler as a player in this game with mathematically defined reward, action, state and environment. The constitutive models are conceptualized as information flow in directed graphs. The process of writing constitutive models is simplified as a sequence of forming graph edges with the goal of maximizing the model score (a function of accuracy, robustness and forward prediction quality). Thus metamodeling can be formulated as a Markov decision process with welldefined states, actions, rules, objective functions and rewards. By using neural networks to estimate policies and state values, the computer agent is able to efficiently selfimprove the constitutive model it generated through selfplaying, in the same way AlphaGo Zero (the algorithm that outplayed the world champion in the game of Go) improves its gameplay. Our numerical examples show that this automated metamodeling framework not only produces models which outperform existing cohesive models on benchmark tractionseparation data, but is also capable of detecting hidden mechanisms among microstructural features and incorporating them in constitutive models to improve the forward prediction accuracy, which are difficult tasksto do manually. We then employ this metalmodeling machine learning techniques to perform offline homogenization for multiphase multipermeability porous media. This model allows one to break down the computational barrier commonly exhibited in DEMFEM and FEM2 models by taking the homogenization procedure offline. As such, we may combine subscale simulation results and experimental data to form a hybridized database that employs one to train and validate deep neural network. A more important implication is that it helps discovering mechanics knowledge that is difficult to see manually without the assistance from the augmented intelligence. For instance, in our published work, we have used the deep learning to help us identifying the importance of multiple microstructural attributes on the macroscopic responses and identify the relative importance of (1) fabric tensor, (2) coordination number of particles and (3) other measurement. Based on the results of the deep learning studies, we then propose models that, for the first time, incorporates the fabric tensors as input for the noncoaxial tractionseparation law and yields a model that is suitable to capture the complex damageplasticity mechanisms of interfaces undergoing a combination of tensile, compressive and shear kinematics (see Figures in the RHS). Our discovery does not only show that the fabric tensor is critical for anisotropic responses of granular materials, but we also find a way to automatically incorporate those influences without handcrafting a new model every times we want to test a new hypothesis. Another related development is the usage of unsupervised learning to generate reducedbasis models (see Figures below). The idea behind this is to use previous simulations of similar natures (e.g. periodic responses, repeated procedures in manufacture, mixing) to determine a Galerkin projection that enables simulations carried with very limited degree of freedoms associated with the orthogonal basis rather than the nodal responses of a standard finite element or discrete element models. We also introduce a number of new ideas in the training procedure. For instance, we introduce the usage of longshortterm memory neuron to predict historydependent responses and employ spectral decomposition to represent evolutions of tensorial input and output information such that the framedependent issues exhibited in RNN constitutive laws are corrected. We use a dualpermeability problem as a test bed to demonstrate how this new modeling approach works. More information can be found in the published articles. Team Members:


Micromechanics of wetted and dry granular materials
The reason we can build sand sculpture is due to the presence of liquid bridges among particles that are otherwise cohesionless. From a micromechanical perspective, liquid bridge and the solid grains are two interacting networks that transmits force and moment. As a result, when subjected to external loadings, granular assemblies with a a very small volume fraction of water can behave profoundly different than the dry counterpart. The objective of this research is to predict how this solidfluid interaction in grain scale affect the macroscopic outcomes. This research is important for a number of engineering applications, including landslides, vehicle mobility and seismic signature of earth materials.
Currently, the research team is working on extending the framework to analyze granular materials beyond the pendular regime using twophase flow model to simulate the airwatersolid interaction in the void space. To make the multiscale model more accessible, the research team has also studied and proposed new methodology to identify material parameters from microCT images. Through collaboration with Prof. Salager and his colleagues, we construct inverse problem that takes account both the macroscopic constitutive responses and microstructural attributes captured by microCT and digital image correlation in the objective function. We are currently expanding this work to consider higherorder (e.g. micropolar) kinematics obtained from microCT. Our preliminary results show that this method may provide a new way to introduce physical length scale directly from the grain particles, without the need to explicitly defining a nonlocal radius or other equivalent nonlocal length scale parameter. Team Members:


Single and polycrystal plasticity of rock salt for nuclear waste disposal
The thermohydromechanical behavior of rock salt is strongly depending on the microstructural attributes and micromechanics of the halite crystal that forms the rock salt. The desired engineering properties, such as the high thermal conductivity and the low permeability, and selfhealing of rock salt are both attributed to the crystalline nature. In a geological system composed of salt or other crystalline rocks, individual grains or crystals are the fundamental building blocks that forms polycrystal assembles. The deformation of these assemblies are then related to both the deformation of each individual crystal and the grain boundary motion. As a result, thermomechanical responses at the microscopic scale, such as intergranular and intragranular fractures, as well as plastic slips, dislocation and the creeping due to precipitation may have significant implications for large scale engineering applications. While there are a multitude of macroscopic phenomenological model attempted to capture the complicated responses of rock salt, a multiscale approach that rightfully treated rock salt as polycrystalline materials may allow more robust predictions on the macroscopic anisotropic responses of rock salt in the finite deformation range. The goal of this work is to derive, formulate, predict, verify and validate polycrystalline model for rock salt and propose multiscale method to bridge the spatial and temporal scales.
Currently, we have already completed the theoretical and algorithmic framework for smallstrain crystal plasticity specifically designed for halite crystal in high pressure high temperature environment. In particular, we adopt the approach in which an elastic region in the stress space is defined and the sequence of the slip activation of the single crystal halite is considered. The analysis includes the thermomechanical creeping and microcracking for singly crystals under different loading conditions (orientation, loading rate, temperature, etc.). As a first step, we are interested at the anisotropic damageplasticity of the slip system in Halite. Using the effective stress hypothesis, we employ a multiphasefield model to capture the evolution of anisotropic damage of the single crystal, while using the crystal plasticity model to capture the plastic slip developed under different temperature and loading rate. The figures describe the axial stressstrain behavior under different orientations. The orientations of the single crystal of halite are described by Euler angles defining crystal axes (xc, yc, zc) relative to the reference coordinate system (x, y, z). Our next goal is to extend it to polycrystal problems and incorporate more mechanisms such as healing and migrations along grain boundaries.
The equivalent plastic strain developed during biaxial compression test for single crystal of different orientations are shown above.
The anisotropic phase field that represents damage of the single crystal are shown above.
Team Members:
Team Members:
 SeonHong Na, 1/2017  current
 S. Na, W.C. Sun, Computational thermomechanics of crystalline rock salt for nuclear waste disposal. Part I: a combined multiphasefield/crystal plasticity approach for single crystal simulations, Computer Methods in Applied Mechanics and Engineering, 2018.
Environmental driven fractures and compaction banding of porous brittle solids
The goal of this study is to analyze the solidfluid interaction of cohesivefrictional and brittle materials across length scale. Our goal is to analyze the relationships between hydraulic properties (such as permeability, hydraulic and mechanical apertures) and microstructural attributes, (such as connectivity of force chain and bonding strengths) before and after the onset of fractures. At the grainscale level, We apply graph theory to unveil how the topology and geometry of force chain and flow network leads to the complex macroscopic outcome, such as fracturing and strain localization.
A recent extension of work includes the derivation of a new variational formulation that incorporate the eigenfracture model for fluiddriven fracture problem and the introduction of eigendeformation mode to simulate the propagation of compaction band and shearenhanced compaction band as a Mode I anticrack. We also propose a new way to compute the hydraulic aperture from the nonlocal phase field fracture and eigenfracture model and validate the numerical simulations with NooruMohamed tests (see Figure on the RHS). Team Members:

Computational modeling of frozen porous media at finite strain
Freezethaw cycle of porous media can lead to significant damage in concrete or pavement. Interestingly, due to the surface energy difference, water thinfilm often develops in between solid grains and ice crystal. At a sufficiently large scale in which RVE exists, the frozen porous media can be regarded as a threephase mixture in which ice, water and solid grains all occupied a fraction of volume. However, depending the temperature and pressure, the ice and water constituents may have phase change and that leads to mass exchange. In particular, the thawing of the ice crystal in frozen soil may lead to a very weak and wet soil layers that are likely to undergo large deformation. Nevertheless, the finite strain behavior of the frozen soil with unfrozen water flowing inside the pores is not well understood. The goal of this study is to formulate a new mathematical framework to capture the mechanics of the frozen porous media that undergoing phase transition at the geometrical nonlinear regime. A climatecontrolled triaxial apparatus will be utilized to study the thermohydromechanical coupling effect of the frozen soil. A preliminary study has been conducted by PhD student SeonHong Na to study the strain localization occurred in frozen soil at finite strain. This work was presented at Engineering Mechanics Institute at Nashville in which SeonHong has won the 2nd place in the best paper competition.
Team Members:
Team Members:
 SeonHong Na, 09/2015  current
 S. Na, W.C. Sun, Computational thermohydromechanics for multiphase freezing and thawing porous media in the finite deformation range, Computer Methods in Applied Mechanics and Engineering, 2017. [URL]
Nonlocal and higherorder continuumdiscrete Coupling for porous media
We propose a nonlocal multiscale framework that couples grainscale discrete element simulations with a macroscopic continuumbased finite element model. The upshot of this nonlocal coupling model is that it retains the simplicity and efficiency of the continuumbased finite element model, while possessing the original length scale of the granular system. In Liu et al. 2015, the collective mechanical responses of particles at material points are homogenized via a staggered nonlocal operator applied on local regions such that the multiscale simulations exhibit no pathological mesh dependence. However, the local regions must be sufficiently large in order to collect enough integration points to compute the nonlocal strain in the finite element model.
As a result, the research group has made other attempts to incorporate length scales, such as introducing secondorder homogenization techniques. For instance, one may consider the incorporation of the higherorder kinematics, such as microrotation, microstretch and microtorsion, and formulate a HillMandel Lemma with more one one energyconjugated pairs. Team Members:


Configurational Poromechanics with Liealgebra internal variable recovery
The motivation of this research is to use the configuration force on poromechanics problems to help solving some of the common problems encountered in modeling material failures. For instance, configurational force can be used as a criterion to detect fracture in brittle materials or used as a criterion to detect needs for remeshing. While this ability has been exploited since the 1990s, the application to poromechanics problems is relatively rare.
To demonstrate this idea, we are currently working on a critical state plasticity model that exhibits sizedependent anisotropy. Because of this sizedependent effect, remeshing is a helpful tool to make simulations results more efficient. Typically, remeshing of problems dealing with pathdependent models is difficult due to (1) the need to mapping history variables and (2) the establishment of new equilibrium after the remeshing. To resolve (1), we use Lie algebra. As demonstrated previously in Wang and Sun, CMAME, 2018 andd earlier in Mota et al. 2014, the quality of projecting tensorial quantity that represent the history of the material is the key for any adaptions. The fields in the formulation are the deformation mapping, the target or mapped internal variables and a Lagrange multiplier that enforces the equality between the source and target internal variables. This formulation leads to an L2 projection that minimizes the distance between the source and target internal variables as measured in the L2 norm of the internal variable space.
To ensure that the target internal variables remain in their original space, their interpolation is performed by recourse to Lie groups, which allows for direct polynomial interpolation of the corresponding Lie algebras by means of the logarithmic map. Once the Lie algebras are interpolated, the mapped variables are recovered by the exponential map, thus guaranteeing that they remain in the appropriate space.
Team Members:
Related Publication:
To demonstrate this idea, we are currently working on a critical state plasticity model that exhibits sizedependent anisotropy. Because of this sizedependent effect, remeshing is a helpful tool to make simulations results more efficient. Typically, remeshing of problems dealing with pathdependent models is difficult due to (1) the need to mapping history variables and (2) the establishment of new equilibrium after the remeshing. To resolve (1), we use Lie algebra. As demonstrated previously in Wang and Sun, CMAME, 2018 andd earlier in Mota et al. 2014, the quality of projecting tensorial quantity that represent the history of the material is the key for any adaptions. The fields in the formulation are the deformation mapping, the target or mapped internal variables and a Lagrange multiplier that enforces the equality between the source and target internal variables. This formulation leads to an L2 projection that minimizes the distance between the source and target internal variables as measured in the L2 norm of the internal variable space.
To ensure that the target internal variables remain in their original space, their interpolation is performed by recourse to Lie groups, which allows for direct polynomial interpolation of the corresponding Lie algebras by means of the logarithmic map. Once the Lie algebras are interpolated, the mapped variables are recovered by the exponential map, thus guaranteeing that they remain in the appropriate space.
Team Members:
 SeonHong Na, 09/2017current
 Kun Wang, 09/2015current
 Eric Bryant, 09/2017current
Related Publication:
 S. Na, E.C. Bryant, W.C. Sun, Capturing size effects of anisotropic fluidinfiltrating materials via a higherorder Camclaytype model with propertypreserving adaptive meshing, in preparation.
 Mota, Sun, Ostien, Foulk, Long, LieGroup interpolation and variational recovery for internal variables, Computational Mechanics, 2013.
Machinesoil Interaction with meshless method
Compaction of road materials, such as soils, aggregate bases are often completed by using vibratory rollers. In recent years, modern technology in sensing, GPS meeting and feedback control has allowed realtime compaction monitoring and adaptive adjustments to the compaction process. Truly predictive computer simulations of this adaptive process may help engineers design more accurate procedure to handle highly heterogeneous soil layers and avoid premature pavement failure. Nevertheless, this process is not suitable for conventional finite element simulations due to the large deformation on the upper layer of the soil. As a result, a stabilized meshless method will be adopted in this study to analyze the rollersoil interaction during the compaction. The simulations shown below is done by former visiting scholar Zhilin Liu during his stay at Columbia.
Team Members:
Team Members:
 Nikolaos Vlassis, Fall 2017  current
 TBA
Concurrent and sequential coupling of multiphysics problems at finite strain
We design and implement of iterative sequential and monolithic schemes to simulate a variety of multiphysics problems. Our objective is to derive general framework to simulate a variety of multiphysics processes that involve solid mechanics coupled with hydraulic, heat and chemical transport. Since conventionally equalorder finite element space often lead to numerical instability, we formulated adaptively stabilized schemes that eliminate spurious oscillations while preserving correct boundary layer thickness. In addition, a critical state plasticity model is being implemented to replicate the constitutive responses of sands under various fluid saturation condition.
We have recently written a review on the thermohydromechanical problem for the September issue of the imechanica journal club. The full text can be found at http://imechanica.org/node/17096.
Team Members:
We have recently written a review on the thermohydromechanical problem for the September issue of the imechanica journal club. The full text can be found at http://imechanica.org/node/17096.
Team Members:
 SeonHong Na, 09/2014  current
 Kun Wang, 07/2014  current
 Eric C. Bryant, 09/2016  current
 Federica Ronchi (visiting student from Università degli Studi di Perugia)
 O.I. Ulven*, W.C. Sun, A twoway coupling dualgraph lattice model for fluiddriven fracture, International Journal for Numerical and Analytical Methods in Geomechanics, 2018.
 S. Na, W.C. Sun, Thermohydromechanical coupling effects on dynamic wave propagation in a fully saturated softening porous medium, International Journal for Numerical and Analytical Methods in Geomechanics, doi:10.1002/nag.2505, 2016. [PDF]
 W.C. Sun, A Stabilized finite element formulation for monolithic thermohydromechanical simulations at finite strain, International Journal for Numerical Methods in Engineering, 2015.
 W.C. Sun, Q. Chen, J.T. Ostien, Modeling hydromechanical responses of strip and circular footings on saturated collapsible geomaterials, Acta Geotechnica, 2014.
 W.C. Sun J.T. Ostien, A. Salinger, A stabilized assumed deformation gradient finite element formulation for strongly coupled poromechanical simulations at finite strain, International Journal for Numerical and Analytical Methods in Geomechanics, 2013.
 W.C. Sun, An unified method to predict diffuse and localized instabilities in sands, Geomechanics and Geoengineering, 8(22):6575, 2013.
Multimodel concurrent analysis with domain overlapping method
The macroscopic mechanical failures, such as fracture and strain localization, are to a large extent determined by physical processes which occur at a scale several orders of magnitude smaller than the macroscopic observation. In order to model this inherent multiscale character efficiently, models designed specifically for macro and microscales are often blended together in a handshake domain. Bridging the multimodel, multiscale domain is, however, not a trivial problem. Spurious wave reflection, ghost force, artificial energy dissipation/increase, mesh dependence and incompatible deformation modes are all challenges persisted for this multimodel approach. To allude these drawbacks, we introduced an infsup stable domain overlapping method and extended it for large deformation, pathdependent plasticity and poroplasticity problems. To circumvent mesh dependence due to the loss of ellipticity, we introduced nonlocal plasticity model and concurrently coupled it with classical continuum model for farfield responses. Such a numerical treatment allows one to concentrate computation resource in a highly localized region of interest and thus significantly increases efficiency.
Graduate student(s):
Graduate student(s):
 Kun Wang, 09/2015current
 SeonHong Na, 09/2015current
 Eric Bryant, 09/2016current
 Sun, Mota, A multiscale overlapped coupling formulation for large deformation strain localization, Computational Mechanics, 2015.
 W.C. Sun, Z. Cai, J. Choo, Mixed Arlequin method for multiscale poromechanics problems, International Journal for Numerical Methods in Engineering, 2016.
Multiscale hierarchical analysis with geometrical enhancements from microstructures: material parameter identification across length scales
Analysis from our studies shows that hierarchical multiscale flow simulation that neglects pore connectivity may lead to significant error and reduce efficiency. As a result, we introduced a new method based on graph theory to extract geometrical attributes from Xray CT images. Information of geometrical attributes is used to improve the coupling among the porescale lattice Boltzmann and field scale finite element simulations. By applying this improved scheme, my collaborators and I quantified the relations between grain crushing and hydraulic property changes. Such information is important for analyzing material bifurcation that leads to fracture and deformation band formation. This work helps improving our understanding on mechanical/hydraulic interactions of solid skeleton and porefluid, which is highly relevant to carbon dioxide storage, nuclear waste management and reservoir management.
Team Members:

Estimating coseismic deformation with combined deterministicstochastic Monto Carlos simulations
My previous research project at Stanford University focuses on estimating coseismic deformation during the 1989 Loma Prieta earthquake. In this study, Professor Borja and I investigated the influence of noise on the estimated coseismic deformation using a combined deterministicstochastic approach. We compare the classical equivalent linear method with a deterministic nonlinear finite element scheme with stochastic input. By quantifying the impact of noise and baseline offset presented in the input ground motion, we conclude that the sediment deformation may have affected the accuracy of the previously backfigured fault rupture mechanism.
Related Publications:
Related Publications: