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 among mathematical principles, data 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.
“What can be shown, cannot be said.” ― Ludwig Wittgenstein
Geometrical symbolic regression for interpretable constitutive modeling
Conventionally, symbolic regression is often limited to function with one or two variables (see the benchmark). For mechanics models, especially the one that requires a large number variables to describe complex phenomena, symbolic regression is often out of reach as a tool.
Meanwhile, conventional neural network is highly expressive but lacks the interpretability of a concise analytical expression. In the case where the neural network used to express the function is large, the calculation via neural network also takes up more memory and computational time and hence not suitable for largescale solid mechanics simulations.
In this line of research, we introduce a twostep machinelearning approach that returns mathematical models interpretable by human experts.
In particular, we hypothesize that a mathematical model (e.g., yield surfaces or hyperelastic energy functionals) can be expressed in terms of a set of singlevariable feature mappings obtained from supervised learning. This idea of feature space can be linear (also see the Neural Additive Model) but may also be extended to polynomials (see Figure below).
Meanwhile, conventional neural network is highly expressive but lacks the interpretability of a concise analytical expression. In the case where the neural network used to express the function is large, the calculation via neural network also takes up more memory and computational time and hence not suitable for largescale solid mechanics simulations.
In this line of research, we introduce a twostep machinelearning approach that returns mathematical models interpretable by human experts.
In particular, we hypothesize that a mathematical model (e.g., yield surfaces or hyperelastic energy functionals) can be expressed in terms of a set of singlevariable feature mappings obtained from supervised learning. This idea of feature space can be linear (also see the Neural Additive Model) but may also be extended to polynomials (see Figure below).
A postprocessing step is then used to reinterpret the set of singlevariable neural network mapping functions into mathematical form through symbolic regression (see below). This divideandconquer approach provides several important advantages. First, it enables us to overcome the scaling issue of symbolic regression algorithms. From a practical perspective, it enhances the portability of learned models for partial differential equation solvers written in different programming languages.
Related publications:
 B. Bahmani, W.C. Sun, Discovering interpretable elastoplasticity models via the neural polynomial method enabled symbolic regressions, under review. (Arxiv first submitted on 24th July 2023.) [Arxiv]
B. Bahmani, W.C. Sun, Physicsconstrained symbolic model discovery for polyconvex incompressible hyperelastic materials, under review. (Arxiv first submitted on 6th Oct, 2023.) [Arxiv]
 Bahador Bahmani, 01/2023  current
 Huijian Cai, 09/2023  current
 Nhon Phan, 10/2023  current
Geometric Learning for computational solid mechanics (with manifold data)
Distanceminimization modelfree data approach introduced by Kirchdoerfer & Ortiz provides a completely new paradigm to run physical simulations of deformable without constitutive laws.
Yet, key technical barriers include the high demand on the data in the phase space, the combinatorial search of the discrete data set and the arbitrary notion of norms equipped for the phase space. In this work, we introduce a deep neural network to globally map data from the constitutive manifold onto a lowerdimensional Euclidean vector space. As such, we establish the relation between the norm of the mapped Euclidean vector space and the metric of the manifold and lead to a more physically consistent notion of distance for the material data. This treatment in return allows us to bypass the expensive combinatorial optimization, which may significantly speed up the modelfree simulations when data are abundant and of high dimensions. Meanwhile, the learning of embedding also improves the robustness of the algorithm when the data is sparse or distributed unevenly in the parametric space.
Another idea we have attempted earlier, which is now published in CMAME, is to create a hybridized modelbased/modelfree approach to handle multiphysics problems where data of different constitutive responses are often of different fidelities or if one of the data set does not provide sufficient data for modelfree approach (i.e. Darcy's flow vs. stressstrain relation). Furthermore, to deal with the costly combinatorial search, we have adopted a KD tree search to speed up the local search.
Related publications:
Yet, key technical barriers include the high demand on the data in the phase space, the combinatorial search of the discrete data set and the arbitrary notion of norms equipped for the phase space. In this work, we introduce a deep neural network to globally map data from the constitutive manifold onto a lowerdimensional Euclidean vector space. As such, we establish the relation between the norm of the mapped Euclidean vector space and the metric of the manifold and lead to a more physically consistent notion of distance for the material data. This treatment in return allows us to bypass the expensive combinatorial optimization, which may significantly speed up the modelfree simulations when data are abundant and of high dimensions. Meanwhile, the learning of embedding also improves the robustness of the algorithm when the data is sparse or distributed unevenly in the parametric space.
Another idea we have attempted earlier, which is now published in CMAME, is to create a hybridized modelbased/modelfree approach to handle multiphysics problems where data of different constitutive responses are often of different fidelities or if one of the data set does not provide sufficient data for modelfree approach (i.e. Darcy's flow vs. stressstrain relation). Furthermore, to deal with the costly combinatorial search, we have adopted a KD tree search to speed up the local search.
Related publications:
 B. Bahmani, W.C. Sun, Manifold embedding datadriven mechanics, Journal of Mechanics and Physics in Solids, accepted, 2022. [Video][Preprint]
 B. Bahmani, W.C. Sun, A kdtree accelerated hybrid datadriven/modelbased approach for poroelasticity problems with multifidelity multiphysics data, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2021.113868, 2021. [Video] [PDF]
 B. Bahmani, W.C. Sun, Distancepreserving manifold denoising for datadriven mechanics, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2022.115857, 2022. [URL]
 Bahador Bahmani, 09/2019  current
Geometric Learning for computational solid mechanics (with graph data)
In this work, our goal is to incorporate nonEuclidean microstructural data and create interpretable descriptors to create a graphbased elastoplasticity model. While traditional hyperelasticity models often incorporate homogenized measures of microstructural attributes, such as porosity, averaged orientation of constitutes, these measures cannot reflect the topological structures of the attributes. We fill this knowledge gap by introducing the concept of a weighted graph as a new means to store topological information, such as the connectivity of anisotropic grains in assembles. Then, by leveraging a graph convolutional deep neural network architecture in the spectral domain, we introduce a mechanism to incorporate these nonEuclidean weighted graph data directly as input for training and for predicting the elastic responses of materials with complex microstructures. To ensure smoothness and prevent nonconvexity of the trained stored energy functional, we introduce a Sobolev training technique for neural networks such that stress measure is obtained implicitly from taking directional derivatives of the trained energy functional. By optimizing the neural network to approximate both the energy functional output and the stress measure, we introduce a training procedure the improves efficiency and generalize the learned energy functional for different microstructures.
Related publications:
Related publications:
 N. Vlassis, W.C. Sun, Sobolev training of thermodynamicinformed neural network for interpretable elastoplasticity models with level set hardening, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2021.113695, 2021. [Video][preprint]
 Y. Heider", H.S. Suh, W.C. Sun, An offline multiscale unsaturated poromechanics model enabled by selfdesigned/selfimproved neural network, International Journal for Numerical and Analytical Methods in Geomechanics, doi:10.1002/nag.3196, 2021.
 N. Vlassis, R. Ma", W.C. Sun, Geometric deep learning for computational mechanics Part I: Anisotropic Hyperelasticity, Computer Methods in Applied Mechanics and Engineering, 2020. [arxiv]
 C. Cai, N. Vlassis, L. Magee, R. Ma", Z. Xiong, B. Bahmani, TF Wong, Y. Wang, W.C. Sun, Equivariant geometric learning for digital rock physics. Part I: Estimating formation factor and effective permeability tensors, International Journal for Multiscale Computation and Engineering, accepted, 2021.[arxiv]
 Nikolaos Vlassis, 09/2017  current
 Ran Ma, 01/2019  current
 Bahador Bahmani, 09/2019  current
Sobolev training for interpretable machine learning plasticity model
To circumvent the lack of interpretability of the classical blackbox neural network, we introduce a higherorder supervised machine learning technique that generates components of elastoplastic models such as elasticity functional, yield function, hardening mechanisms, and plastic flow. The geometrical interpretation in the principal stress space allows us to use convexity and smoothness to ensure thermodynamic consistency. The speed function from the HamiltonJacobi equation is deduced from the DNS data to formulate hardening and nonassociative plastic flow rules governed by the evolution of the lowdimensional descriptors. By incorporating a noncooperative game that determines the necessary data to calibrate material models, the machine learning generated model is continuously tested, calibrated, and improved as new data guided by the adversarial agents are generated. A graph convolutional neural network is used to deduce lowdimensional descriptors that encodes the evolutional of particle topology under pathdependent deformation and are used to replace internal variables. The resultant constitutive laws can be used in a finite element solver or incorporated as a loss function for the physicalinformed neural network run physical simulations.
Related publications:
 N. Vlassis, W.C. Sun, Sobolev training of thermodynamicinformed neural network for interpretable elastoplasticity models with level set hardening, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2021.113695, 2021. [Video][preprint]
 N. Vlassis, W.C. Sun, Componentbased machine learning paradigm for discovering ratedependent and pressuresensitive levelset plasticity models, Journal of Applied Mechanics, doi:10.1115/1.4052684, 2021.
 Nikolaos Vlassis, 09/2017  current
 Mian Xiao, 09/2020  current
 Bahador Bahmani 09/2019  current
Highstrainrate responses of energetic materials
The thermomechanical response of the energetic material under impacting is closely related to its safety and reliability during processing, storage, and transportation. The goal of this project is to establish a multiscale and multiphysics model to predict the shock loading resistance of the energetic material. In our previous work, we developed a pressuredependent anisotropic hyperelasticity and crystal plasticity model to incorporate pressure sensitivity while predicting the plastic dissipation. It is found that the external pressure has a strong influence on the plastic dissipation under shock loading, which cannot be captured by traditional crystal plasticity model. The material point method (MPM) is a hybrid EulerianLangrangian method, which is immune to mesh distortion and hence suitable for large deformation simulation. Also, MPM provides automatic nopenetration noslip contact algorithm which enables the porecollapse simulation. In our preliminary results with JohnsonCook model, it is found that the heatinduced strain softening and shear band can be reproduced with MPM. We are currently working on introducing the anisotropic phase field fracture model into the MPM to simulate the fracture near the collapsing pores and related heat generation.
Related publications:
Related publications:
 R. Ma, W.C. Sun, C. R. Picu, Atomisticmodel informed pressuresensitive crystal plasticity for crystalline HMX, International Journal of Solid and Structures, 232:111170, doi:10.1016/j.ijsolstr.2021.111170, 2021.
 R. Ma, W.C. Sun, A finite microrotation material point method for micropolar solid and fluid dynamics with threedimensional evolving contacts and free surfaces, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2021.114540, 2022.
 Ran Ma, 09/2019  current
 Nick Vlassis, 09/2020current
 Mian Xiao, 09/2021current
Metamodeling for computational mechanics: modeling the process of writing a material model
Our goal in this research is to not write material models for simulation but simulate how modelers write a model by interpreting data. In many machine learning applications for physics problems, the majority of time is spent on solving the optimization problems (tuning hyperparameters, balancing constraints for boundary and governing laws). Here our objective is to take a more human approach where we also attempt to build an interpretable explanation abstracted by knowledge graphs.
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. We use a dualpermeability problem as a test bed to demonstrate how this new modeling approach works. 

Team Members:
 Qing Yin 03/2021  current
 Kun Wang, 09/2014  09/2018
 Nikolaos Vlassis, 09/2017  current
 Xinran Zhong, 09/2016  08/2018
 Yousef Heider 10/201812/2020
 Jarett Poliner, 2020  current
 K. Wang, W.C. Sun, A multiscale multipermeability poroplasticity model linked by recursive homogenizations and deep learning, Computer Methods in Applied Mechanics and Engineering, accepted, 2018. [PDF]
 K. Wang, W.C. Sun, Metamodeling game for deriving theoreticalconsistent, microstructuralbased tractionseparation laws via deep reinforcement learning, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2018.11.026, 2019. [URL]
 K. Wang, W.C. Sun, An updated Lagrangian LBMDEMFEM coupling model for dualpermeability porous media with embedded discontinuities, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2018.09.034, 2010.
 K. Wang, W.C. Sun, Q. Du, A cooperative game for automated learning of elastoplasticity knowledge graphs and models with AIguided experimentation, Computational Mechanics, special issue for DataDriven Modeling and Simulations: Theory, Methods and Applications, 64(2):67–499, doi:,10.1007/s00466019017231, 2019. [PDF]
 Y. Heider", K. Wang, W.C. Sun, SO(3)invariance of graphbased deep neural network for anisotropic elastoplastic materials, Computer Methods in Applied Mechanics and Engineering, 2020.
 A. Fuchs, Y. Heider", K. Wang, W.C. Sun, M. Kaliske, DNN2: A hyperparameter reinforcement learning game for selfdesign neural network elastoplastic constitutive laws, Computer and Structures, 249:106505, doi:10.1016/j.compstruc.2021.106505, 2021.
 X. Sun, B. Bahmani, N. Vlassis, W.C. Sun, Y. Xu, Datadriven discovery of interpretable causal relations for deep learning material laws with uncertainty propagation, Granular Matter, doi:10.1007/s1003502101137y, 2021.
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:
 Mian Xiao, 09/2020current
 Chuanqi Liu, 09/201809/2020
 Ran Ma, 01/2021current
Published Work:
 C. Liu, W.C. Sun, Shift domain material point method for solids in the finite deformation range, special thematic issue for Meshfree and Particle Methods for Modeling Extreme Loadings, Computational Particle Mechanics, doi: 10.1007/s4057101900239y, 2019. [PDF]
 C. Liu, W.C. Sun, ILSMPM: an unbiased implicit levelsetbased material point method for frictional particulate contact mechanics of deformable particles, Computer Methods in Applied Mechanics and Engineering, accepted, doi:10.1016/j.cma.2020.113168, 2020. [PDF]
 R. Ma, W.C. Sun, A finite microrotation material point method for micropolar solid and fluid dynamics with threedimensional evolving contacts and free surfaces, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2021.114540, 2022.
 M. Xiao, W.C. Sun, C. Liu, DPMPM: Domain partitioning material point method for evolving multibody thermalmechanical contacts and fragmentation, Computer Methods in Applied Mechanics and Engineering, 385:114063, doi:10.1016/j.cma.2021.114063, 2021. [Video]
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 porous 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:
 Kun Wang, 09/2014  09/2018
 Hyoung Suk Suh, 09/2018  current
 Zeyu Xiong, 09/2020  current
 Qing Yin, 03/2021  current
 K. Wang, W.C. Sun, A micromechanical study on anisotropy of capillary stress in wetted granular matters, accepted, Journal of Engineering Mechanics, 2015.
 M.R. Kuhn, W.C. Sun, Q. Wang, Stressinduced anisotropy in granular materials: fabric, stiffness, and permeability, Acta Geotechnica, 2015.
 K. Wang, W.C. Sun, A semiimplicit discretecontinuum coupling method for fluidsaturated porous media at finite strain, Computer Method in Applied Mechanics and Engineering, 2016. [PDF]
 K. Wang, W.C. Sun, S. Salager, S. Na, G. Khaddour, Identifying material parameters for a micropolar plasticity model via Xray microCT images: lessons learned from the curvefitting exercises, accepted, International Journal of Multiscale Computational Engineering, 2016. [DRAFT]
 W.C. Sun, TF. Wong, Prediction of hydraulic and electrical transport properties of sandstone with multiscale lattice Boltzmann/finite element simulation on microtomographic images, International Journal of Rock Mechanics and Mining Science, 2018. [PDF]
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. In a more recent work, our postdoc research scientist Dr. Ran Ma has introduced two FFT solvers to simulate brittle fracture with strongly anisotropic fracture energy. By leveraging the global continuity of the Fourier series, we introduce a generalized frequency vector to compute the higherorder derivative required for the higherorder phase field. This method is compared with the multiphasefield formulation where multiple phase field evolution equations and the mechanics problems are solved in a staggered update scheme. Results obtained from both phase field models are compared. The FFT provides a simple way to solve the problems without specific mesh generation and treatments to resolve the higherorder terms.
 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.
 R. Ma, W.C. Sun, FFTbased solver for higherorder and multiphasefield fracture models applied to strongly anisotropic brittle materials and polycrystals, Computer Methods in Applied Mechanics and Engineering, accepted, 2019.
 R. Ma, W.C. Sun, Computational thermomechanics for crystalline rock. Part II: chemodamageplasticity and healing in strongly anisotropic polycrystals, under review.
In a more recent work, our postdoc research scientist Dr. Ran Man extends the FFT solver to incorporate the effect of healing, the precipitation creeping due to the transport, the permeability evolution and the crystal plasticity of polycrystalline salt (see right Figure). To capture the healing process that often begins at the crack tips, an indicator function is used to locate the crack tip and the interfacial diffusion is captured explicitly along the grain boundaries to faithfully capture the solutionprecipitation creeping of polycrystals. The resultant ratedependence and the effect the evolution of damage and plastic deformation are shown below. A manuscript will be submitted in the near future (est. Oct 2019) to document this work.
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, 09/2018  current
 Qing Yin 03/2021  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]
 H. Suh, W.C. Sun, Asynchronous phase field fracture model for porous media with thermally nonequilibrated constituents, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2021.114182, 2021. [Preprint]
 H. Suh, W.C. Sun, Multiphasefield microporomechanics model for simulating ice lens growth and thaw in frozen soil, under review.
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:
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/20159/2019
 Eric Bryant, 09/201708/2020
 S. Na, E.C. Bryant, W.C. Sun, A configurational force for adaptive remeshing of gradientenhanced poromechanics problems with historydependent variables, Computer Methods in Applied Mechanics and Engineering, doi:10.1016/j.cma.2019.112572, 2019.
 A. Mota, W.C. Sun, J.T.Ostein, J.W. Foulk III, K.N. Long, LieGroup interpolation and variational recovery for internal variables, Computational Mechanics, 52:12811299, doi:10.1007/s0046601308761, 2013. [PDF] [Bibtex]
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:
Our research at Columbia is proudly sponsored by:
Disclaimer: Electronic copies of the preprints and articles are provided for convenience only. Before downloading, please check to make sure you have the right permission to download them from your institution. Unless specified otherwise, all materials are subjected to the Attribution 4.0 International (CC BY 4.0) license. If you use/incorporate/reproduce/borrow/develop any new derivation, data, figure, idea or information from this archive into your own work, either consciously or subconsciously, my students and I would greatly appreciate your decency by giving the appropriate credit where credit's due. We appreciate any comment, suggestion or critique, please contact wsun@columbia.edu. Thank you for your interest in our work.