Stress recovery for the particle-in-cell finite element method

Abstract

The interelement stress in the Finite Element Method is not continuous in nature, and stress projections from quadrature points to mesh nodes often causes oscillations. The widely used particle-in-cell method cannot avoid this issue and produces worse results when there are mixing materials of large strength (e.g., viscosity in Stokes problems) contrast in one element.The post-processing methods including (1) distance weighted average from surrounding particles to the centroid mesh node (Post-local), (2) global projection with least square fit (Post-global), and (3) superconvergent point recovery method (SPR), cannot effectively eliminate the stress fluctuations. We propose three pre-processing methods to reduce the interface contrast in mixing elements: (1) global method with harmonic-mean averaging (GHM), (2) unification of properties at mixed-material elements (UnE), and (3) averaging particle properties within a specified distance to gauss quadrature points (AGP). For tests of Q1 elements, the results processed by combining either pre-processing method with the Post-local projection can increase the precision. The GHM pre-processing method is the least computationally expensive application and the easiest to implement, the AGP pre-processing method is the most expensive and the UnE in-between. However, for Q2 elements, the GHM pre-processing method fails in stress recovery, and produces worse results than those without any pre-processing procedures. For general cases (both Q1 and Q2 elements), the AGP preprocessing method is recommended. The optimal sampling radius used in the AGP method is close to that size of one element, beyond which it increases computational time, but does not significantly increases the accuracy of recovered stresses. In terms of the averaging approaches used in the AGP method, the harmonic mean is suitable for simple-sheardominated processes and the arithmetic mean is better for the pure-shear-dominated models. For complex models, the AGP method of harmonic mean combined with the SPR post-procedure is recommended. The AGP method is found to be able to efficiently reduce stress perturbations in a synthetic model of complex fault geometries like the San Andreas Fault system.

Intermediate-depth seismicity and intraslab stress changes due to outer-rise faulting

“…3. In order to minimize the effects of local spurious stress oscillations characteristic to particle in cell methods with strong local viscosity variations 36,37 , we averaged deviatoric stress values over sufficiently large (50 × 50 km, i.e., 2500 grid cells, Fig. 2b-j) domains.…”

Numerical Modeling of Earthquake Cycles Based On Navier‐Stokes Equations With Viscoelastic‐Plasticity Rheology

“…In this benchmark, the vertical fault is conformed to one element, and the fault material is not mixed with wall rocks in any elements. We further test models with a fault of one and a half element width and do not find intensive stress fluctuations as is common for particle‐in‐cell methods when elements are filled with materials of high viscosity contrast (Yang, Moresi, & Mansour, 2021). We find that the stress field is also comparable with that of community models, but the maximum slip rate is almost twice of that from community models.…”

“…Although the fault is designed to be of 1.5 element width, only central parts of the fault zone occupy one element, and both fault‐wall rock interfaces are in elements containing two materials. Better ways to estimate the effective strain rate for the entire fault zone (Yang, Moresi, & Mansour, 2021) may address this issue but is beyond the topic of this study.…”

“…The finite width of the fault ( w f ) is of one element size in x direction, starting from x = 0 and ending in x = Δ x , where Δ x is the one mesh size in x direction. This setup makes the fault conform to the mesh, and no mixed material in one element affects the stress estimation (Yang, Moresi, & Mansour, 2021). The fault‐normal, along‐strike and along‐dip dimensions of the computational model are marked as L 1, L 2, and L 3, respectively.…”

“…The displacement of faults in long‐term tectonic models are averaged over one time step of several millenniums, which approximates the cyclic accumulation of earthquake events. Faults designed in long‐term tectonic models have a finite thickness and are free to deform in geometries, especially for modeling with the Particle‐in‐cell method (Harlow, 1964; Moresi et al., 2003; Yang, Moresi, & Mansour, 2021). The long‐term tectonic models are driven by speculated external boundary conditions and internal body forces and can evolve in a self‐consistent fashion.…”