WIXLIB

Computers and Geotechnics 149 (2022) 104855S. Zhao et al.typically employed to simulate the response of a representative volumeelement (RVE) as a surrogate of the conventional constitutive model ata material point to feed a continuum-based method, such as FEM (Guoand Zhao, 2014; Liu et al., 2016; Desrues et al., 2019), SPFEM (Guoet al., 2021) and MPM (Liang and Zhao, 2019), for solution of a boundary value problem. The hierarchical coupling structure offered in theseapproaches enables effective simulations of various engineering problems with complex boundary and initial conditions and yet effortlesslycapture the complicated granular material responses that are highlyloading path and state dependent. The hierarchical multiscale couplingscheme can be a promising alternative for thermo-mechanical modelingof granular media. Indeed, a hierarchical coupling of FEM and DEM forthermo-mechanical modeling was proposed in our previous study (Zhaoet al., 2020), which is, however, limited to thermally induced smalldeformation and stress regime. It is desirable to develop an alternativemultiscale approach that can rigorously tackle TM coupling of granularmedia in the regime of large deformation and flow.This work aims at developing an entirely new hierarchical multiscale framework to simulate coupled thermo-mechanical behavior inlarge deformation and flow of granular media. The framework will bebuilt upon a pure mechanical hierarchical coupling formulation of MPMand DEM (termed as DEMPM) that was developed for modeling largedeformation of granular media (Liu et al., 2017; Liang and Zhao, 2019;Zhao et al., 2021). In DEMPM, a material point serves as an interfacefor information exchange between the macro-scale MPM solver andthe meso-scale DEM solver. Notably, a velocity gradient at a givenmaterial point can be obtained first from the MPM solver. It is thenpassed on to the corresponding RVE assembly as constraints for a DEMsolution which is further used to extract required material responses(e.g., stress) to feedback to the MPM solver. The new development incorporating thermo-mechanical coupling will be termed as TM-DEMPMhereafter. To further consider thermo-mechanical coupling in diversified conditions, the generalized interpolation material point method(GIMP) (Bardenhagen and Kober, 2004) will be adopted in TM-DEMPMas the continuum MPM solver to solve both thermal and mechanicalfields in a staggered manner at the continuum scale, whereas physicsbased thermo-mechanical coupling is fully considered and rigorouslymodeled at the RVE scale in the DEM solver. As will be demonstrated,the proposed new framework presents three outstanding features whengranular media develop deformation en route to large deformation andflow regime: (1) it is able to model thermally induced stress due to grainexpansion/contraction; (2) it can capture mechanically induced heatdue to inter-particle friction dissipation; and (3) it can handle thermallyinduced melting at the contact scale. Consequently, TM-DEMPM offersa computational pathway for physics-based, cross-scale modeling andunderstanding of the thermo-mechanical responses of granular media.The rest of the paper is organized as follows. Section 2 providesthe macroscopic formulation of TM modeling by GIMP, followed by themesoscopic formulation of TM coupling in DEM in Section 3. Section 4introduces the hierarchical multiscale coupling and computing scheme.Numerical examples are showcased in Section 5, where the TM-MPMsolver is benchmarked against analytical and FEM solutions prior to thevalidation of TM-DEMPM, followed by three examples representativeof typical TM problems in granular media further presented to demonstrate the predictive capabilities and features of the proposed approach,including heat generation in biaxial compression, heat generation indischarge of a granular silo, and melting induced column collapse.Concluding remarks are summarized in Section 6. Tensorial indicialnotations and Einstein summation convention are followed in the study,and boldface letters for matrices are used.geostructures during large deformation and failure of granular media.In agricultural or powder industry, continuous discharging of grainsin a silo can accumulate heat to cause remarkable increase of temperature, posing a serious trigger for dangerous silo explosion (Russoet al., 2017). In chemical engineering, impregnation and calcinationperformance depends crucially on heat transfer and flow properties ofgranular media in rotary vessels (Saruwatari and Nakamura, 2022).For frozen granular media in nature, the melting process of interparticle bonds can significantly weaken the strength of granular matrix,triggering large deformation and flow of soils that may cause catastrophic geohazards and environmental disasters (Maslin et al., 2010).For example, the dissociation of gas hydrates in submarine settingsmay trigger potential geohazards during the course of exploitation,ranging from submarine landslides to earthquakes and tsunamis, threatening the operation and safety of offshore infrastructures and coastalcities (Maslin et al., 2010). Thaw-induced landslides in permafrost areasserve both as an awakening warning of global warming and an exacerbating factor by causing carbon release into the atmosphere (Froitzheimet al., 2021). There are evidently pressing needs from both industrial developments and environmental sustainability to advance ourunderstanding of thermo-mechanical response of granular media, inparticular when they reach the regime of large deformation and/orflow.Great efforts have been devoted to the study of thermo-mechanicalcoupling in granular media from diversified communities of granularphysics, thermal processing industry, and chemical and geotechnicalengineering. A prevailing body of these studies have been based oncontinuum mechanics based models and numerical approaches. Theytypically smear the discrete nature of a granular material and consider it as a continuum body with homogenized macroscopic responsesdescribed by assumed phenomenological constitutive models (Na andSun, 2017; Liu et al., 2018). To treat an engineering scale boundary value problem, these models are routinely implemented into acontinuum-based numerical approach such as Finite Element Method(FEM) and its alternative mesh-free ones such as Material Point Method(MPM) (Sulsky et al., 1995) and Smoothed Particle Hydrodynamics(SPH) (Gingold and Monaghan, 1977) which are specifically useful forlarge deformation problems (Bui et al., 2008; Lei et al., 2021). Despitetheir great success, continuum constitutive models may face difficultiesin capturing the physical behavior of granular media under complexloading conditions entangled with complicated coupling. When dealingwith problems involving large deformation and flow, even greaterchallenges are posed to both model developments and their implementations. Micromechanics-based approaches and tools, representedby Discrete Element Method (DEM) (Cundall and Strack, 1979), haverecently triggered growing interests by their capability of reproducingvarious complex mechanical characteristics of granular media from aparticulate perspective, such as anisotropy and liquefaction (Guo andZhao, 2013), strain localization and failure (Chen et al., 2011), noncoaxiality (Li and Yu, 2015) and the rich transitional behavior betweenfluid and solid (Herrmann et al., 2013). Micromechanics approachescan also facilitate the consideration of various coupling processes.Indeed, DEM may accommodate the thermo-mechanical coupling withproper network or pore-scale models (Nguyen et al., 2009; Gan et al.,2012; Choo et al., 2013; Moscardini et al., 2018; Caulk et al., 2020).It can further be coupled with other computational schemes such asComputational Fluid Dynamics (CFD) and Lattice Boltzmann Method(LBM) (Zhang et al., 2016) to simulate coupled hydro-mechanicalbehavior of granular media. However, computational cost may bean outstanding issue for DEM that limits its practical application forlarge-scale problems.Recent spotlights of modeling of granular media have more beenattracted by a class of hierarchical multiscale approaches that leveragethe strengths of both continuum- and discrete-based methods (Andradeet al., 2011; Chen et al., 2011; Guo and Zhao, 2014; Liu et al., 2016).They employ a hierarchical numerical coupling scheme where DEM is2. Thermo-mechanical modeling by GIMP2.1. Governing equations and weak form formulationThe macroscopic thermo-mechanical response of granular mediais modeled in a continuum-based manner, by solving the following2

Computers and Geotechnics 149 (2022) 104855S. Zhao et al.balance equations of energy and linear momentum in an updatedLagrangian framework:𝜌𝑐𝑇,𝑡 𝑞𝑖,𝑖 𝑄𝑡(1)𝜎𝑖𝑗,𝑗 𝜌𝑏𝑖 𝜌𝑎𝑖(2)physical field 𝑓 (𝑥𝑖 ) can be obtained by interpolating the correspondingvalues 𝑓𝑝 carried on the material points in terms of the so-called particlecharacteristic function 𝜒𝑝 (𝑥𝑖 ) (Bardenhagen and Kober, 2004), i.e., 𝑓𝑝 𝜒𝑝 (𝑥𝑖 ),(7)𝑓 (𝑥𝑖 ) 𝑝where 𝜌 is the bulk mass density; 𝑐 is the specific heat capacity; 𝑇 isthe temperature; 𝑡 is the time; 𝑞𝑖 is the heat flux; 𝑄𝑡 is a heat source orsink; 𝜎𝑖𝑗 is the Cauchy stress tensor; 𝑏𝑖 is the body force per unit of masspossibly performed on materials (e.g., gravitational acceleration), and𝑎𝑖 is the acceleration term; 𝑖, 𝑗 are sequential indices in {1, 2, 3} for 3Dor {1, 2} for 2D. Note that a latent heat term can be added in Eq. (1)to further consider the effect of phase transition on heat transfer. Toavoid unnecessary defocusing of the study, the detailed implementationis not presented in this paper. In future considerations, special attentionhas to be paid to the different thermal behaviors of coarse and finesoils during phase transition. Notably, the freezing point is widelyconsidered stationary for coarse soils during phase change, whilst it isregarded as a function of temperature for fine soils (Michalowski andZhu, 2006).According to the Fourier’s law for anisotropic heat conduction, theheat flux 𝑞𝑖 can be written as (Onsager, 1931)such that the volume of a material particle is given by𝑉𝑝 Note that 𝜒𝑝 (𝑥𝑖 ) satisfies the partition of unity property in the undeformed configuration (Bardenhagen and Kober, 2004) and a ‘‘top-hat’’function is employed here. For the background grid, a physical field𝑔(𝑥𝑖 ) can be also interpolated based on the corresponding value 𝑔𝐼 atthe grid node 𝐼, i.e., 𝑔(𝑥𝑖 ) 𝑔𝐼 𝑁𝐼 (𝑥𝑖 )(9)𝐼where 𝑁𝐼 (𝑥𝑖 ) is the grid shape function at node 𝐼. Note that thesubscripts 𝑝 or 𝐼 denote properties or functions associated with particle𝑝 or node 𝐼 hereafter.The weak forms in Eqs. (5) and (6) can be rewritten into thefollowing integration forms over all particles] (𝑇 ) [ (𝑇 ) 𝜔𝐼𝑉𝑝 𝜙𝐼𝑝 𝑄𝑡𝑝 𝑞̄𝑝 ℎ 1 𝜙𝐼𝑝𝜔𝐼𝑚𝑝 𝑐𝑝 𝜙𝐼𝑝 𝑁𝐽 𝑝 𝑇𝐽 ,𝑡 (3)𝑞𝑖 𝑘𝑖𝑗 𝑇,𝑗where 𝑘𝑖𝑗 is the thermal conductivity tensor.The following general boundary and initial conditions are considered for a typical heat conduction and mechanical problem:𝑇 (𝑡) 𝑇̄on 𝛤𝑇(4a)𝑞(𝑡)𝑖 𝑞̄𝑖on(4b)𝛤𝑞𝑢𝑖 𝑢̄ 𝑖𝜎𝑖𝑗 𝑛𝑗 𝑡̄𝑖𝛤𝑢(4d)on 𝛤𝑡on(4e))𝜔(𝑇,𝑖 𝑘𝑖𝑗 𝑇,𝑗 d𝛺 𝛺 𝛺𝜔(𝑢)𝑖,𝑗 𝜎𝑖𝑗 d𝛺 𝜔(𝑇 ) 𝜌𝑐𝑇,𝑡 d𝛺 𝛺𝜔(𝑢)𝑖 𝜌𝑎𝑖 d𝛺 𝛺 𝛺𝜔(𝑇 ) 𝑄𝑡 d𝛺 𝜔(𝑢)𝑖 𝜌𝑏𝑖 d𝛺 𝛤𝑞 𝛤𝑡𝜔(𝑇 ) 𝑞d𝛤̄(5)̄𝜔(𝑢)𝑖 𝑡𝑖 d𝛤(6)𝐼 𝐼𝑝)𝜔(𝑇𝐼 𝑝𝑉𝑝 𝜙𝐼𝑝,𝑖 𝑘𝑝𝑖𝑗 𝑇𝑝,𝑗𝑝(10) 𝜔(𝑢)𝐼𝐼 𝐽𝑚𝑝 𝜙𝐼𝑝 𝜙𝐽 𝑎𝑖𝐽 𝑝 𝐼 𝜔(𝑢)𝐼 𝑝𝜔(𝑢)𝐼𝑚𝑝 𝜙𝐼𝑝 𝑏𝑖𝑝 𝐼 𝜔(𝑢)𝐼 𝑚𝑝 𝜙𝐼𝑝 𝑡̄𝑖𝑝 ℎ 1𝑝𝐼𝑉𝑝 𝜙𝐼𝑝,𝑗 𝜎𝑝𝑖𝑗𝑝(11)where 𝛤𝑇 and 𝛤𝑞 are the prescribed temperature and heat flux boundaries of the problem domain 𝛺, respectively; 𝑇̄ and 𝑞̄𝑖 are the prescribedboundary temperature on 𝛤𝑇 and boundary heat flux on 𝛤𝑞 , respectively; 𝑇𝑟𝑒𝑓 is the ambient or reference temperature; Eqs. (4a) and(4b) are Dirichlet (fixed temperature) and Neumann (fixed heat flux)boundary conditions, respectively, while Eq. (4c) is the initial conditionfor temperature distribution at the initial state. 𝑛𝑗 is the boundaryoutward normal of the domain 𝛺; 𝑢̄ 𝑖 and 𝑡̄𝑖 are the prescribed materialdisplacement on 𝛤𝑢 and boundary traction on 𝛤𝑡 , respectively; Eqs. (4d)and (4e) are Dirichlet and Neumann boundary conditions, respectively.Multiplying Eqs. (1) and (2) by a corresponding test function 𝜔(𝑇 )(with zeros on 𝛤𝑇 ) and 𝜔(𝑢)(with zeros on 𝛤𝑢 ), respectively, and𝑖integrating over the entire domain 𝛺 by application of integration byparts and Green’s formula, the weak form formulations of the governingequations can be obtained as follows 𝛺𝑝𝐽𝐼(4c)𝑇 (0) 𝑇𝑟𝑒𝑓(8)𝜒𝑝 d𝛺. 𝛺𝑝where ℎ is the virtual boundary layer thickness serving for the boundary integration (Zhang et al., 2011); 𝑚𝑝 is the mass of particle 𝑝; 𝜙𝐼𝑝 and𝜙𝐼𝑝,𝑖 are the weighting and gradient weighting functions, respectively,defined as𝜙𝐼𝑝 1𝜒 𝑁 d𝛺𝑉𝑝 𝛺𝑝 𝑝 𝐼(12)𝜙𝐼𝑝,𝑖 1𝜒 𝑁 d𝛺𝑉𝑝 𝛺𝑝 𝑝 𝐼,𝑖(13)Although several other GIMP variations have been proposed in theliterature for better accuracy, e.g., CPDI (Sadeghirad et al., 2011) andCPDI2 (Sadeghirad et al., 2013), the undeformed GIMP (uGIMP) (Bardenhagen and Kober, 2004) is adopted here as the macroscopic solverof the proposed multiscale framework for simplicity in implementationwithout losing generality.Due to the arbitrariness of the test functions, the weak forms givenby Eqs. (10) and (11) can be simplified into the following semi-discreteequations at node 𝐼: 𝑖𝑛𝑡𝐶𝐼𝐽 𝑇𝐽 ,𝑡 𝑄𝑒𝑥𝑡(14)𝐼 𝑄𝐼2.2. Spatial discretization in GIMP𝐽 Lagrangian material points and Eulerian background grid are twokey ingredients for MPM discretization of a continuum domain. TheGIMP formulation was proposed by Bardenhagen and Kober (2004) toovercome the potential cell-crossing issue in the conventional MPM.In the GIMP discretization, the continuum is spatially discretized intoa finite set of subdomains occupied by material points with a certaindomain for each. The integration over the entire domain 𝛺 governedby the weak forms in Eqs. (5) and (6) can be converted into thesummation of integration over each subdomain 𝛺𝑝 of particles. A𝑀𝐼𝐽 𝑎𝑖𝐽 𝑓𝑖𝐼𝑒𝑥𝑡 𝑓𝑖𝐼𝑖𝑛𝑡(15)𝐽with𝐶𝐼𝐽 𝑝𝑄𝑒𝑥𝑡𝐼 𝑉𝑝 𝜙𝐼𝑝 𝑄𝑡𝑝 𝑝𝑓𝑖𝐼𝑒𝑥𝑡 𝑝3𝑀𝐼𝐽 𝑚𝑝 𝑐𝑝 𝜙𝐼𝑝 𝑁𝐽 𝑝 𝑚𝑝 𝜙𝐼𝑝 𝑏𝑖𝑝 𝑝𝑚𝑝 𝜙𝐼𝑝 𝜙𝐽 𝑝(16)𝑉𝑝 𝜙𝐼𝑝,𝑖 𝑘𝑝𝑖𝑗 𝑇𝑝,𝑗(17)𝑝𝑞̄𝑝 ℎ 1 𝜙𝐼𝑝𝑝 𝑄𝑖𝑛𝑡𝐼 𝑝𝑚𝑝 𝜙𝐼𝑝 𝑡̄𝑖𝑝 ℎ 1𝑓𝑖𝐼𝑖𝑛𝑡 𝑝𝑉𝑝 𝜙𝐼𝑝,𝑗 𝜎𝑝𝑖𝑗(18)

Computers and Geotechnics 149 (2022) 104855S. Zhao et al.where 𝐶 and 𝑀 are the consistent heat capacity matrix and mass matrix, respectively; 𝑄 and 𝑓𝑖 are heat and mechanical loads, respectively(the superscripts ‘‘int’’ and ‘‘ext’’ denote internal and external loads,respectively). To facilitate the computation, the scheme of lumpedmatrix is employed, i.e., 𝑚𝑝 𝜙𝐼𝑝(19)𝑚𝑝 𝑐𝑝 𝜙𝐼𝑝𝑀𝐼 𝐶𝐼 temperature given by Eq. (29) is the so-called remapped temperature,which is smoother than that given by Eq. (28) (Tao et al., 2016). Inthis work, the fluctuating particle temperature is tracked to update thenodal temperature, while the smooth particle temperature is taken asan input for the representative volume elements (RVEs). Detail will beintroduced in Section 4.In addition, the timestep should be chosen sufficiently small tomaintain numerical stability and convergence in the coupled thermomechanical system. In the proposed staggered coupling algorithm, boththe mechanical and thermal solvers adopt the same timestep satisfyingthe following constrains (Tao et al., 2016; Lei et al., 2021):𝑝𝑝such that,𝑖𝑛𝑡𝐶𝐼 𝑇𝐼,𝑡 𝑄𝑒𝑥𝑡𝐼 𝑄𝐼𝑀𝐼 𝑎𝑖𝐼 𝑓𝑖𝐼𝑒𝑥𝑡 𝑓𝑖𝐼𝑖𝑛𝑡(20)𝑙𝛥𝑡 𝑒𝐸 𝜌2.3. Temporal discretizationWe employ the Euler forward method to discretize the thermomechanical coupled system (i.e., the semi-discrete equations inEqs. (20)) in time such that the MPM solver works in an explicitmanner. Indeed, this explicit MPM in conjunction with the schemeof lumped matrix has been prevailingly adopted in the literature forits straightforward numerical implementation (de Vaucorbeil et al.,2020). In this scheme, the discrete equations for nodal temperature andvelocity are written as:(𝑡)𝛥𝑡𝑇𝐼(𝑡 𝛥𝑡) 𝑇𝐼(𝑡) 𝑇𝐼,𝑡(21)(𝑡)(𝑡)𝑣(𝑡 𝛥𝑡) 𝑣𝑖𝐼 𝑎𝑖𝐼𝛥𝑡𝑖𝐼(22) 𝐼(𝑡) (𝑡)𝜙𝐼𝑝 𝑣𝑖𝐼 𝜙(𝑡) 𝛥𝑡𝑎(𝑡)𝑖𝐼 𝐼𝑝for PIC3. Thermo-mechanical coupling in DEM3.1. Thermal response of an RVE assembly3.1.1. Heat generated by frictionHeat can be generated due to friction dissipation at inter-particlecontacts during various granular flow processes. Note that heat canbe generated through different dissipation processes and mechanisms.Only sliding friction is considered here for simplicity as it has been reported experimentally to dominate the friction dissipation (Zhai et al.,2019). For a single timestep 𝛥𝑡, the total amount of heat within anassembly due to friction dissipation is quantified by (𝑐) (𝑐)𝐸𝑓 𝑓𝑡 𝑣𝑡 𝛥𝑡(32)(24)𝑐One major difference between the FLIP and PIC methods is that theFLIP method uses the grid accelerations only while the PIC methoduses the updated grid velocities. Stomakhin et al. (2013) suggesteda combination of the FLIP and PIC methods by regrading PIC as adamping term as follows: (𝑡) (𝑡) ) (𝑡) (𝑡)((𝑡)(25)𝑣𝑖𝐼 𝜙𝐼𝑝𝑎𝑖𝐼 𝜙𝐼𝑝 𝛥𝑡 𝛼𝑃 𝐼𝐶 𝑣(𝑡)𝑣(𝑡 𝛥𝑡) 𝑣𝑖𝑝 𝑖𝑝 𝑖𝑝𝐼where 𝑓𝑡(𝑐) and 𝑣(𝑐)𝑡 are friction (tangential force) and relative tangentialvelocity at contact 𝑐 during the present timestep. Therefore, the heatsource is given by 𝑁 (𝑖)𝑖 1 𝐸𝑓(33)𝑄 𝑉 𝑁𝛥𝑡where 𝑁 is the DEM iteration number for a given strain loading; 𝑉 isthe volume of the assembly.𝐼where 𝛼𝑃 𝐼𝐶 [0, 1] is the PIC fraction. With grid damping 𝛼𝑔 andparticle damping 𝛼𝑝 , the damped particle velocity reads (𝑡) (𝑡) (𝑡) (𝑡)(𝑡)𝑣(𝑡 𝛥𝑡) 𝑣𝑖𝑝 𝑎𝑖𝐼 𝜙𝐼𝑝 𝛥𝑡 (𝛼𝑃 𝐼𝐶 𝛼𝑝 )𝑣(𝑡)𝑣𝑖𝐼 𝜙𝐼𝑝 (26)𝑖𝑝𝑖𝑝 (𝛼𝑃 𝐼𝐶 𝛼𝑔 )𝐼3.1.2. Thermal conductivity tensorWe follow our previous study (Zhao et al., 2020) to constructthe fabric-based conductivity tensor for an assembly of non-sphericalparticles. The derivation is based on such a fundamental assumptionthat heat flows through an imaginary heat pipe joining the center ofeach contacted particle (heat reservoir) through the contact (anotherheat reservoir) as schematically illustrated in Fig. 1(a). Accordingly, athermal conductivity tensor 𝑘𝑖𝑗 can be given as𝐼and the particle position is updated with a general second-order FLIPformulation (Nairn, 2015) (𝑡 𝛥𝑡) (𝑡)𝑥(𝑡 𝛥𝑡) 𝑥(𝑡)𝑣𝑖𝐼 𝜙𝐼𝑝 𝛥𝑡𝑖𝑝𝑖𝑝 𝐼 ( (𝑡) (𝑡) )𝛥𝑡 (𝑡) (𝑡)𝑎𝑖𝐼 𝜙𝐼𝑝 𝛥𝑡 (𝛼𝑃 𝐼𝐶 𝛼𝑝 )𝑣(𝑡) (𝛼𝑃 𝐼𝐶 𝛼𝑔 )𝑣𝑖𝐼 𝜙𝐼𝑝𝑖𝑝2 𝐼𝐼(27)𝑘𝑖𝑗 Similarly to the FLIP and PIC methods on velocity update in themechanical part, two possible schemes can be formulated to update theparticle temperature, i.e., (𝑡) (𝑡)𝑇𝑝(𝑡 𝛥𝑡) 𝑇𝑝(𝑡) 𝜙𝐼𝑝 𝑇𝐼,𝑡 𝛥𝑡(28)𝑇𝑝(𝑡 𝛥𝑡) 𝑙𝑒2 𝜌𝑐𝐼𝐼𝐼(𝑡) (𝑡 𝛥𝑡)𝜙𝐼𝑝 𝑇𝐼(30)for thermal(31)𝑏𝑘where 𝐸 is Young’s modulus, which can be estimated from odometertests for an RVE packing (see Section 5.2.2 for detail); 𝑏 is an adjustablefactor, with a default value of 1; 𝑙𝑒 is the minimum element size of amesh; 𝑘 is the mean thermal conductivity, i.e., average of the diagonalterms of the thermal conductivity tensor 𝑘𝑖𝑗 .𝛥𝑡 where the superscripts 𝑡 and 𝑡 𝛥𝑡 indicate the variables at the startand end of the processing timestep hereafter, respectively. For themechanical part, the movement of material points can be solved according to the kinematic fields of the background grid. Two directcandidate methods are considered here, i.e., the so-called FLIP (FLuidImplicit Particle Brackbill and Ruppel, 1986) method and PIC (ParticleIn Cell Harlow, 1964) method, where particle velocity is updated as (𝑡) (𝑡)(𝑡)𝑣(𝑡 𝛥𝑡) 𝑣𝑖𝑝 𝑎𝑖𝐼 𝜙𝐼𝑝 𝛥𝑡 for FLIP(23)𝑖𝑝 𝑣(𝑡 𝛥𝑡)𝑖𝑝for mechanical𝑀 (𝑝) (𝑝)1 𝑟𝑖 𝑟𝑗𝑉 𝑝 1 𝛼 (𝑝) 𝑙(𝑝)(34)where 𝑉 is the volume of the assembly; 𝑟𝑖 is the vector along the 𝑝thheat pipe; 𝑙 is the branch vector joining the two adjacent particles atthe contact of the heat pipe, and 𝛼 is the thermal resistance per lengthof the heat pipe.We note that it is non-trivial to derive an elegant theoretical solutionfor the thermal conductivity tensor 𝑘𝑖𝑗 of granular media (Choo et al.,2013). Porosity has been regarded as one of the key factors influencing𝑘𝑖𝑗 . The effective thermal conductivity has been frequently assumed to(29)𝐼where the former uses only the incremental temperature of nodes,whereas the latter is a direct remapping of nodal temperatures. The4

Computers and Geotechnics 149 (2022) 104855S. Zhao et al.Fig. 1. (a) Heat conduction through an inter-particle contact; (b) thermally induced contact force due to particle expansion.where 𝑅(0)𝑖 is the initial principal length of a particle at the referencetemperature; 𝑅𝑖 is the current principal length of a particle at a temperature change 𝛥𝑇 with respect to the reference temperature; 𝛽 is thelinear thermal expansion coefficient. As a consequence, the thermallyinduced change in particle profile yields an additional penetration 𝛿𝑑𝑖 ,referring to Fig. 1(b).follow an empirical relation with porosity (Zhang and Wang, 2017),without consideration of the anisotropic nature of materials in anengineering setting, e.g., in geotechnical engineering. Further appraisalis needed to rigorously consider the various physical attributes in deriving a robust expression for effective thermal conductivity of granularmedia.Remarks: In the conventional TM modeling with pure DEM, e.g., (Nguyenet al., 2009; Choo et al., 2013), the inter-particle heat transfer has to bemodeled at the particle scale for heat flows with remarkable increase in computational cost. By contrast, in the present proposed TM coupling scheme,heat transfer is captured at the mesoscopic/RVE scale (the effect of fabricchange can be considered, e.g., on thermal conductivity tensor) and solved inMPM at the continuum scale, which helps to develop the homogenized modelwith TM coupling and significantly boost the computational efficiency.3.2.3. Thermally sensitive bond contact modelFor cemented granular media such as frozen soils in permafrostareas and gas hydrate-bearing soils in submarine settings, particlesare bonded together through the ice or ice-like cementation. Such icybonds are sensitive to temperature change which may cause phasetransformation. Such process of inter-particle bonds can be readilyconsidered via proper modification of the contact models (Brown et al.,2014; Shen et al., 2016). For simplicity yet without losing generality,we introduce the following thermally sensitive bond contact model thatgoverns the bond breakage upon heating:3.2. Mechanical response of an RVE assembly3.2.1. Basic contact modelDEM modeling of granular media commonly postulates a force–displacement law in conjunction with Coulomb’s friction conditiongoverning inter-particle contacts (Cundall and Strack, 1979). A simplelinear spring contact model is employed here to establish the relation between contact force and relative displacement for each pair ofcontacting grains as follows𝑓𝑖 𝑛 𝑘𝑛 𝑢𝑛𝑖𝛥𝑓𝑖 𝑘𝑡 𝛿𝑢𝑡𝑖𝑛(35b)𝑡where 𝑓𝑖 and 𝑓𝑖 are normal and tangential contact forces along thenormal 𝑛𝑐𝑖 and the tangential 𝑡𝑐𝑖 directions at contact 𝑐, respectively,referring to Fig. 1(b); 𝛥𝑓𝑖 𝑡 is the incremental tangential contact forceat the current timestep; 𝑘𝑛 and 𝑘𝑡 are the normal and tangent contactstiffness, respectively; 𝑢𝑛𝑖 is the penetration depth along contact normal, 𝛿𝑢𝑡𝑖 is the relative tangential displacement of the two contactingparticles at the current timestep. The Coulomb condition of friction isapplied to constraining the tangential contact force according to, 𝒇 𝑡 𝜇 𝒇 𝑛 (38) 𝐹𝑡 𝐹𝑡𝑏 ,𝐹𝑡𝑏 𝜇𝑏 𝐹𝑛𝑏(39)where 𝛽𝑎 is a tuned parameter, and 𝑇𝑠 is a threshold for temperaturechange. Both 𝛽𝑎 and 𝑇𝑠 are assumed to be constant in the numericalexamples in this work. More elegant models will be investigated in ourfuture work.(36)where 𝜇 is the coefficient of friction.3.2.2. Thermal expansion contact modelFor a dry cohesionless granular assembly, the thermally inducedstress within an assembly depends on expansion or contraction ofparticles caused by temperature change. For simplicity, the change ofparticle size is assumed to be governed by the following relation (Vargas and McCarthy, 2007; Zhao et al., 2017; Zhao and Feng, 2019):𝑅𝑖 𝑅(0)𝑖 (1 𝛽𝛥𝑇 )𝐹𝑛𝑏 𝛽𝑏 𝑘𝑛 𝑅where 𝛽𝑏 is defined as the breakage strength coefficient with respectto normal contact stiffness 𝑘𝑛 ; 𝑅 is the average radius of the twocontacting particles; 𝜇𝑏 is a threshold coefficient of the maximum shearforce at breakage with respect to the normal breakage force, whichis assumed independent of temperature. Once either condition of theabove inequalities is met, the bond breaks and can no longer sustaintension. Note that herein we only consider bond breakage subjectedto tension and shear, while neglecting bond failure caused by twistingand/or compression. The breakage strength coefficient is a function oftemperature 𝑇 , which can be simply defined as{ 𝛽 (𝑇 𝑇 )𝑎 𝑠, if 𝑇 𝑇𝑠 ,𝑇𝑠𝛽𝑏 (40)0,otherwise(35a)𝑡𝐹𝑛 𝐹𝑛𝑏 ,3.2.4. Periodic cell, deformation and stressAn RVE assembly is assumed to be sufficiently small to capturethe response of a material point from the macroscopic perspective andmeanwhile sufficiently large to yield a representative response withrespect to discrete grains. To reduce the boundary effect from rigidboundaries (e.g., rigid confining walls), periodic boundary conditionsare employed in the DEM simulations (Thornton, 2000; Yang et al.,2014; Radjai, 2018). Typically,

chemical and mining engineering (Nguyen et al.,2009). A better under-standing of the interplay between thermal and mechanical behaviors of granular media could provide crucial theoretical bases for the design, operation, and risk assessment of relevant applications and practices. The thermo-mechanical coupling not only manifests by direct be-