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.
Metamodeling for computational mechanics
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. Many machine learning models used in mechanical computer simulations employ supervised and unsupervised learning techniques to make predictions with little or no interpretability. This disconnection with human knowledge makes it difficult to assess and manage risk when using such blackbox models for engineering applications. In this talk, we present a new framework in which the AI is not simply tasked to make predictions on a mechanical process, but to mimic the actions of a human engineer to create, validate and test handcrafted models that make predictions. As such, we introduce the concept of metamodeling in which two AI agents, an experimentalist that generates data and a modeler that interprets data, are tasked with building a directed graph that represents the knowledge of the physical process with the least amount of data required. In this metamodeling framework, all data accessible by the agents is modeled as labelled vertices (e.g. porosity, permeability) and all possible ways to connect those vertices are modeled as directed edges (e.g. mathematical expression, neural nets, support vector machine). Therefore, all the possible actions of the agents are mathematically represented by a labeled directed multigraph. The learning process is therefore recast as the procedure to find a directed graph that links the input and output data with the optimal set of edges that maximize the reward (the objective function that measures the prediction qualities). With welldefined agents, reward, action space and environment, the AI modeler will then improve its ability to write the models through selfpracticing in a deep reinforcement learning framework. Like humans, it will attempt to take different actions to write a model, then estimate the reward of each action via a deep neural network in a Monte Carlo tree search. By leveraging the ability of the AI modeler to repetitively practice and improve its modeling skills, we demonstrate that the proposed algorithm is able to both discover new hidden hierarchical structures of mechanics knowledge and rediscover wellknown mechanics knowledge without any human intervention. In the case where availability of data is limited, the metamodeling algorithm also explores the weakness of links in the constitutive laws and explores the optimal set of experiments to yield the best forward predictions under a limited budget. 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. Originally written/last updated 1/31/2019. Team Members:

Score distribution of the selfderived elastoplasticity models from five families 1: Generalized plasticity with critical state, 2: Generalized plasticity w/o critical state, 3. Classical nonassociative plasticity with critical state, 4. Classical nonassociative plasticity w/o critical state, 5. Others.

Modified material point method for complex geometries undergoing finite deformation
One key drawback of classical discrete element methods for particulate mechanics is that they cannot generate any kinetic and kinematics information at the subgrain scales. Nevertheless, many dissipative mechanisms of granular materials are actually originated from pathdependent behaviors at the subgrain scales where wear, damage, fracture and fragmentation occurs. Nevertheless, recovery of these subscale information while capturing the geometry and form of each individual grains is inherently a difficult task.
The goal of this research is to establish an efficient imagetosimulation workflow to overcome this technical barrier. We introduce a mathematical framework designed to enable a simple imagetosimulation workflow for solids of complex geometries undergoing large deformation in the geometrically nonlinear regime. In particular, we adopt the integration scheme of the material point method to resolve the convergent issues for Lagrangian meshes due to mesh distortion, while using a shifted domain technique originated from Scovazzi and coworker to represent the boundary condition implicitly via a level set or signed distance function. Consequently, this method completely bypasses the need to generate highquality conformal mesh to represent complex geometries and therefore allows modelers to select the space of the interpolation function without the constraints due to the geometrical need. This important simplification enables us to simulate deformation of complex geometries inferred from voxel images. Verification examples on deformable body subjected to finite rotation have shown that the new shifted domain material point method is able to generate frameindifference results. Meanwhile, simulations using microCT images of a Hostun sand have demonstrated that this method is able to reproduce the quasibrittle damage mechanisms of single grain without the excessively concentrated nodes commonly displayed in conformal meshes that represent 3D objects with local fine details.
Originally written/last updated 1/31/2019.
Team Members:
Published Work:
The goal of this research is to establish an efficient imagetosimulation workflow to overcome this technical barrier. We introduce a mathematical framework designed to enable a simple imagetosimulation workflow for solids of complex geometries undergoing large deformation in the geometrically nonlinear regime. In particular, we adopt the integration scheme of the material point method to resolve the convergent issues for Lagrangian meshes due to mesh distortion, while using a shifted domain technique originated from Scovazzi and coworker to represent the boundary condition implicitly via a level set or signed distance function. Consequently, this method completely bypasses the need to generate highquality conformal mesh to represent complex geometries and therefore allows modelers to select the space of the interpolation function without the constraints due to the geometrical need. This important simplification enables us to simulate deformation of complex geometries inferred from voxel images. Verification examples on deformable body subjected to finite rotation have shown that the new shifted domain material point method is able to generate frameindifference results. Meanwhile, simulations using microCT images of a Hostun sand have demonstrated that this method is able to reproduce the quasibrittle damage mechanisms of single grain without the excessively concentrated nodes commonly displayed in conformal meshes that represent 3D objects with local fine details.
Originally written/last updated 1/31/2019.
Team Members:
 Chuanqi Liu, 09/2018current
Published Work:
 C. Liu", W.C. Sun, Shift domain material point method for solids in the finite deformation range, under review.
Particle damage occurs during indentation simulated by a nonlocal damage model. Note that classical DEM cannot recover any stress field or energy flux information near the singular tip and hence cannot accurately predict the propagation of crack and crack branching consistent with fracture mechanics theory.
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. Originally written/last updated 1/31/2018. 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.
Originally written/last updated 12/31/2018.
Originally written/last updated 12/31/2018.
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  12/2018
 Ran Ma, 1/2019current
 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). Originally written/last updated 12/31/2018. 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.
Originally written/last updated 12/31/2018.
Team Members:
Originally written/last updated 12/31/2018.
Team Members:
 SeonHong Na, 09/2015  12/2018 (now assistant professor at McMaster University)
 Alessandro Milleri, 09/2018  current
 Hyoung Suk Suh, 90/2018  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/201712/2018 (now assistant professor at McMaster University)
 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  12/2018 (now assistant professor at McMaster University)
 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: