# An explicit multi-time stepping algorithm for multi-time scale coupling
problems in SPH

Xiaojing Tang<EMAIL_ADDRESS>Dong Wu<EMAIL_ADDRESS>Zhengtong Wang
<EMAIL_ADDRESS>Oskar Haidn<EMAIL_ADDRESS>Xiangyu Hu
<EMAIL_ADDRESS>TUM School of Engineering and Design, Technical University
of Munich, 85748 Garching, Germany

###### Abstract

Simulating physical problems involving multi-time scale coupling is
challenging due to the need of solving these multi-time scale processes
simultaneously. In response to this challenge, this paper proposed an explicit
multi-time step algorithm coupled with a solid dynamic relaxation scheme. The
explicit scheme simplifies the equation system in contrast to the implicit
scheme, while the multi-time step algorithm allows the equations of different
physical processes to be solved under different time step sizes. Furthermore,
an implicit viscous damping relaxation technique is applied to significantly
reduce computational iterations required to achieve equilibrium in the
comparatively fast solid response process. To validate the accuracy and
efficiency of the proposed algorithm, two distinct scenarios, i.e., a
nonlinear hardening bar stretching and a fluid diffusion coupled with Nafion
membrane flexure, are simulated. The results show good agreement with
experimental data and results from other numerical methods, and the simulation
time is reduced firstly by independently addressing different processes with
the multi-time step algorithm and secondly decreasing solid dynamic relaxation
time through the incorporation of damping techniques.

###### keywords:

Smoothed particle hydrodynamics, Multi-time scale coupling , Multi-time step
algorithm , Dynamic damping , Multi-physics problem

##  1 introduction

Smoothed Particle Hydrodynamics (SPH), a typically mesh-free method, which is
originally introduced by Lucy [1], Gigold and Monaghan [2] for studying
astrophysical problems, has been widely applied to simulate fluid-flows [3, 4,
5, 6], solid mechanics [7, 8, 9, 10, 11], fluid-structure interaction [12, 13,
14] in recent years. Comprehensive reviews can be found in Refs. [15, 16, 17,
18, 19]. Even with wide applications, SPH has some limitations when it comes
to simulating multi-scale coupling problems existing in various engineering
fields, particularly those involving solid dynamic response which is a
typically fast process [20]. The disparity in the time scales of fast and slow
processes presents a continuing challenge to numerical simulations [21].

To solve multi-time scale problems, either an implicit or explicit scheme can
be applied. The implicit scheme allows for a larger time step in the time
integration [22, 23], enabling the monolithic scheme to solve the equations
for all fast and slow processes simultaneously. For instance, Zhao [24] used
an implicit Newmark scheme to model the flow through a porous elastic solid,
where solid dynamics and fluid diffusion occur at different time scales.
Gaston [25] employed an implicit scheme to analyze the fluid, chemistry, and
structure coupling behavior in a reactor, which is a common phenomenon in the
engineering field. However, since the inversion of the stiffness matrix used
for solving equations is required for each time step [26, 27], this approach
is quite expensive concerning both computation time and memory consumption
[28].

The explicit approach is more favorable for solving multi-time scale coupling
problems due to its direct time integration and simple numerical formulation
[29, 30, 31, 32]. Some researchers have used this approach to simulate
material stretching and necking, where the load is applied during a long time
period while the material’s dynamic response is instant and fast [33, 34, 35].
Since the realistic load is applied in a long time scale, a long physical
simulation time is expected. However, with a quite small stable time step size
allowed in explicit scheme for the fast process, usually millions of time
steps are required to simulate the entire process, which is very often not
feasible. To reduce the overall simulation time, loading rate is usually
increased artificially [34]. However, high non-realistic loading rate may lead
to certain limitations and inaccuracies in the simulation results [36].

This paper presents a multi-time stepping algorithm in SPH, where a large and
a small time steps are chosen according to the slow and fast processes in the
simulation, respectively. Two, i.e., an outer and an inner loops are arranged
with these two time steps for time integration. Specifically, the slow process
is integrated with a large time step in the outer loop, while the fast solid
dynamic process with a much smaller time step in the inner loop. Since the
time step size of the fast process is small, many iterations of the solid
stress relaxation may occur within one outer loop and lead to low
computational efficiency. To address this issue, a dynamic relaxation method
based on implicit operator splitting scheme [37] is adopted to accelerate the
convergence rate of the fast dynamic process to an elastic equilibrium state.
To assess the performance and computational efficiency of the proposed
algorithm, the simulations of tensile tests, including two dimensional and
three dimensional cases, are firstly carried out; and then the evolution of
fluid diffusion in porous media coupling with elastic deformation is
simulated. The latter fluid-structure coupling process occurs in chemical
reactors, e.g. in the fuel cell of battery, where fluid mixture diffuses
through a Nafion membrane, affecting the battery performance due to the
varying fluid concentration and membrane deformation. The obtained results
demonstrate that the proposed algorithm performs better both in accuracy and
efficiency compared to previous numerical methods.

The reminder of this paper is organized as follows. Section 2 summaries the
theories and governing equations for nonlinear hardening plastic solid
mechanics and fluid-structure interaction. Section 3 describes the
corresponding SPH discretization. In Section 4, the proposed multi-time
stepping algorithm coupling with the dynamic relaxation are detailed. Section
5 states the physical problems and the results obtained using the proposed
algorithm are compared with those from previous methods and experiments.
Finally, Section 6 presents brief concluding remarks. The source code and data
needed for this numerical simulation work can be found in SPHinXsys, an open-
source multi-physics SPH library, available at https://www.sphinxsys.org.

##  2 Governing equations

###  2.1 Total Lagrangian solid dynamics

In this section, we provide a concise introduction to solid dynamics within
the framework of total Lagrange formulation, along with the relevant notations
and and symbols that will be utilized in the subsequent models. The analysis
focuses on a solid body $\mathcal{B}$, which occupies two regions:
$\mathcal{R}_{0}$ and $\mathcal{R}$, representing the body’s configurations at
time $t_{0}$ ($t=0$) and $t$ respectively. In the initial configuration
$\mathcal{R}_{0}$, the position vector of a material point is represented by
$\mathbf{X}\in\mathcal{R}_{0}$, while in the current configuration, it is
denoted as $\mathbf{x}\in\mathcal{R}$. The motion of the solid body is
described by the invertible mapping $\phi$, which transforms a material point
$\mathbf{X}$ to its corresponding vector $\mathbf{x}=\phi(\mathbf{X},t)$, as
illustrated in Figure. 2.1. Based on this definition, the Lagrangian velocity
of a material point is defined as
$\mathbf{v}(\mathbf{X},t)=\frac{d\phi(\mathbf{X},t)}{dt}$. The deformation
gradient $\mathbf{F}$, which characterizes the deviation of a material point
from its initially undeformed position to its deformed position, can be
computed from the displacement vector $\mathbf{u}=\mathbf{x}-\mathbf{X}$ using
the following equation:

$\mathbf{F}=\frac{d\mathbf{x}}{d\mathbf{X}}=\nabla^{0}{\mathbf{u}}+\mathbf{I},$
(1)

where $\mathbf{I}$ is the unit matrix, and the superscript
$\left({\bullet}\right)^{0}$ accounts for quantities in the initial reference
configuration. The corresponding Jacobian determinant term $J=$
det($\mathbf{F}$) indicates the local volume gain $J>1$ or loss $J<1$.

Figure 2.1: Finite deformation process on a body $\mathcal{B}$.

The governing equations of solid deformation within the total Lagrange
framework are derived as

$\begin{cases}\rho={\rho^{0}}\frac{1}{J}\\\
\rho^{0}\frac{\text{d}\mathbf{v}}{\text{d}t}=\nabla^{0}\cdot\mathbf{P}^{T}\end{cases},$
(2)

where $\rho$ and $\rho_{0}$ are the densities in the current configuration
$\mathcal{R}$ and the initial configuration $\mathcal{R}_{0}$ respectively,
$\mathbf{v}$ the velocity and $\mathbf{P}$ the first Piola-Kirchhoff stress
tensor. Different from the Cauchy stress $\bm{\sigma}$, which points to the
force measured in the deformed configuration, $\mathbf{P}$ relates to stress
within the initial configuration, and the two stresses are related by

$\mathbf{P}=J\bm{\sigma}\mathbf{F}^{-T}=\bm{\tau}\mathbf{F}^{-T},$ (3)

where $\bm{\tau}$ is the Kirchhoff stress tensor, which is obtained from the
constitutive relation as given in Appendix A. Also, using a multiplicative
decomposition technique [38, 39], a hardening plastic model is also given in
Appendix A.

###  2.2 Fluid-structure interaction

For the fluid diffusion in porous media coupling with elastic deformation of
the porous membrane, we propose a fluid-structure interaction model, where the
fluid diffuses in the porous solid, leading to an increased fluid pressure and
solid deformation.

Figure 2.2: Partially saturated porous medium.

In this model, the heterogeneous body is considered as a continuous solid
medium containing uniformly distributed small voids with a homogeneous
porosity ${a}$. When this medium comes into contact with a fluid, fluid flows
into these small pores and diffuses inside this medium due to the presence of
the fluid concentration gradient, resulting in the formation of a mixture
comprising solid and fluid components, as illustrated in Figure 2.2. To
simplify this model, we adopt the methodology proposed by Zhao [24] to to
present a mixture momentum equation while fluid behaviors follow the diffusion
law.

####  2.2.1 Mass and momentum equations

With a porosity ${a}$ and fluid saturation level $\widetilde{a}$ (see Appendix
B.1), the locally effective fluid density $\rho^{l}$ can be expressed as

$\rho^{l}=\rho^{l}_{0}\widetilde{a},$ (4)

where ${\rho}^{l}_{0}$ is the initial density of the fluid. The governing
equations for the solid body involving the density conservation is described
as

$\rho^{s}={\rho^{s}_{0}}\frac{1}{J},$ (5)

where $\rho^{s}$ and $\rho_{0}^{s}$ are the solid density defined in current
configuration $\mathcal{R}$ and initial configuration $\mathcal{R}_{0}$
respectively, For a porous solid partially-saturated by fluid, the total
linear momentum $\mathbf{M}$ in the region $\mathcal{R}$ is the sum of fluid
momentum and solid momentum

$\mathbf{M}=\rho\mathbf{v}=\rho^{l}\mathbf{v}^{l}+\rho^{s}\mathbf{v}^{s},$ (6)

where $\rho$, $\mathbf{v}$ is the total density and velocity, $\mathbf{v}^{l}$
the velocity of fluid, $\mathbf{v}^{s}$ the velocity of dry porous solid. Due
to the difference between $\mathbf{v}^{l}$ and $\mathbf{v}^{s}$, the fluid
flux $\mathbf{q}$ on the element boundary $\partial V$ can then be expressed
as

$\mathbf{q}={\rho^{l}}(\mathbf{v}^{l}-\mathbf{v}^{s}).$ (7)

Obviously, if there is no fluid passing through the boundary, $\mathbf{q}=0$,
the fluid mass in an element is conserved. The transfer of fluid mass and
momentum between micro-scale solid constituents happens when fluid flows from
regions with higher fluid saturation to those with lower saturation.
Therefore, within an element $dV$ of the mixture, the balance of linear
momentum implies that the time derivative of momentum $\mathbf{M}$ is
determined by two factors. One is the stress exerting on the element and the
other one is the fluid flux of linear momentum
$\mathbf{v}^{l}\otimes\mathbf{q}$ on the boundary $\partial V$, where the
symbol $\otimes$ means an outer product of two vectors or tensors. It follows
that the conservation of total linear momentum of the mixture can be expressed
as

$\frac{D\mathbf{M}}{Dt}=\nabla\cdot\bm{\sigma}-\nabla\cdot\left(\mathbf{v}^{l}\otimes\mathbf{q}\right),$
(8)

where $\bm{\sigma}$ represents the cumulative Cauchy stress in the mixture
acting on the solid. $\bm{\sigma}$ is determined by Cauchy stress
$\bm{\sigma}^{s}$ and the pressure stress due to the presence of the fluid
phase $\bm{\sigma}^{l}$, which is detailed in Appendix B.2.

####  2.2.2 Fick’s law

In a partially saturated solid, the fluid saturation difference leads to the
motion of fluid from higher fluid fraction to lower parts and the flux follows
the Fick’s law

$\mathbf{q}=-K\rho^{l}\nabla\widetilde{a},$ (9)

indicating that the fluid flux is proportional to the diffusivity $K$, the
effective fluid density $\rho^{l}$ as well as the gradient of the fluid
saturation $\widetilde{a}$. Consequently, the time derivative of fluid mass in
an element $dV$ is due to the fluid flux $\mathbf{q}$ on the element boundary
$\partial V$, written as

$\frac{D\rho^{l}}{Dt}=-\nabla\cdot\mathbf{q}.$ (10)

##  3 SPH implementation

In SPH, the continuum is represented by a set of Lagrangian particles that
carry various properties, such as mass, position, velocity, and other
attributes. A variable field is approximated using a kernel function that
represents the influence of neighboring particles and the mechanics of the
continuum are approximated by modeling the interactions between these
particles. In this section, we transform the governing equations of two
previously discussed models into SPH discretization.

###  3.1 SPH discretization for solid dynamics

To discretize the solid mechanics, we employ the initial undeformed
configuration as the reference. First, aiming to restore 1st order
consistency, a correction matrix $\mathbf{B}^{0}$ [40, 10] of particle $a$ is
adopted as

$\mathbf{B}^{0}_{a}=\left(\sum_{b}V_{b}\left(\mathbf{r}^{0}_{b}-\mathbf{r}^{0}_{a}\right)\otimes\nabla^{0}_{a}W_{ab}\right)^{-1},$
(11)

where $V_{b}$ represents the volume of the neighboring particle $b$,
$\mathbf{r}^{0}_{a}$ and $\mathbf{r}^{0}_{b}$ denote the positions of
particles $a$ and $b$ in the reference configuration, and
$\nabla^{0}_{a}W_{ab}$ is the gradient of the kernel function given by

$\nabla^{0}_{a}W_{ab}=\frac{\partial
W\left(|\mathbf{r}^{0}_{ab}|,h\right)}{\partial|\mathbf{r}^{0}_{ab}|}\mathbf{e}^{0}_{ab},$
(12)

where $\mathbf{e}^{0}_{ab}$ is a unit vector pointing from particle $a$ to
$b$. In total Lagrangian formulation, the neighborhood of particle $a$ is
defined in the initial configuration, and this set of neighboring particles
remains fixed throughout the entire simulation. Additionally,
$\mathbf{B}^{0}_{a}$ is computed only once under the initial reference
configuration. The momentum conservation in Eq. (2) can be approximated in the
strong form as

$\frac{\text{d}\mathbf{v}_{a}}{\text{d}t}=\frac{2}{\rho_{a}}\sum_{b}V^{0}_{b}\tilde{\mathbf{P}}_{ab}\nabla^{0}_{a}W_{ab},$
(13)

where $\rho_{a}$ represents the density of particle $a$,
$\tilde{\mathbf{P}}_{ab}$ is the averaged first Piola-Kirchhoff stress of the
particle pair $(a,b)$, stated as

$\tilde{\mathbf{P}}_{ab}=\frac{1}{2}\left(\mathbf{P}_{a}\mathbf{B}^{0}_{a}+\mathbf{P}_{b}\mathbf{B}^{0}_{b}\right).$
(14)

Note that the first Piola-Kirchhoff stress tensor is dependent on the
deformation tensor $\mathbf{F}$, the time derivative of which is computed from

$\frac{d\mathbf{F}_{a}}{dt}=\left(\sum_{b}V_{b}\left(\mathbf{v}_{b}-\mathbf{v}_{a}\right)\otimes\nabla^{0}_{a}W_{ab}\right)\mathbf{B}^{0}_{a},$
(15)

where $\mathbf{v}_{a}$ and $\mathbf{v}_{b}$ denote the velocities of particles
$a$ and $b$. Considering the plastic response which may exist in the solid
deformation, a return mapping algorithm is used to obtain the stress-strain
evolution.

###  3.2 SPH discretization for fluid-structure interaction

In the fluid-structure interaction model discretization, each particle carries
the location $\mathbf{x}_{n}=\phi(\mathbf{X},t_{n})$ at time $t_{n}$, along
with an initial representative volume $V^{0}$ that partitions the initial
domain of the macroscopic solid. The deformation gradient $\mathbf{F}_{n}$ of
the solid phase is stored to update the solid current volume $V_{n}$ and
density $\rho^{s}_{n}$. Additionally, the fluid mass $m_{n}^{l}$, saturation
$\widetilde{a}_{n}$, and density-weighted velocity of the fluid relative to
solid $\mathbf{q}_{n}$ are stored. The fluid mass equation Eq. (10) of
particle $i$ is discretized as

$\frac{\text{D}m_{i}^{l}}{\text{D}t}=2V_{i}\sum_{j}\frac{m_{j}}{\rho_{j}}(\mathbf{q}_{i}-\mathbf{q}_{j})\nabla_{i}W_{ij}.$
(16)

Note that with the equation 1, we have the relation of gradient kernel
function in the total Lagrangian and updated Lagrangian
$\nabla_{i}W_{ij}=\mathbf{F}^{-1}\nabla_{i}^{0}W_{ij}$. Once fluid mass is
updated, the locally effective fluid density $\rho^{l}$ is obtained
subsequently. According to Eq. (39) and Eq. (9), we update the fluid
saturation $\widetilde{a}$ and the fluid flux $\mathbf{q}$ in the particle
form

$\mathbf{q}=-K\rho^{l}V_{i}\sum_{j}\frac{m_{j}}{\rho_{j}}(\widetilde{a}_{i}-\widetilde{a}_{j})\nabla_{i}W_{ij}.$
(17)

With the fluid flux and the stress in hand, we obtain discrete formulations
for the momentum balance equation Eq. (8) as

$\frac{D\mathbf{M}_{i}}{Dt}=2\sum_{j}V_{j}(\mathbf{T}_{i}+\mathbf{T}_{j})\nabla_{i}W_{ij}-2\sum_{j}V_{j}({\mathbf{v}_{i}^{l}\otimes\mathbf{q}_{i}}+{\mathbf{v}_{j}^{l}\otimes\mathbf{q}_{j}})\nabla_{i}W_{ij},$
(18)

where $\mathbf{T}_{i}$ and $\mathbf{T}_{j}$ are the stress tensors between
particles $i$ and $j$. We then compute the updated solid velocity
$\mathbf{v}^{s}$ using the total momentum definition Eq. (6), where the total
density of the mixture is the sum of the solid and fluid densities
$\rho=\rho^{s}+\rho^{l}$, written as

$\mathbf{v}^{s}=\frac{\mathbf{M}-\mathbf{q}}{\rho}=\frac{\mathbf{M}-\mathbf{q}}{\rho^{s}+\rho^{l}}.$
(19)

Subsequently, the fluid velocity $\mathbf{v}^{l}$ is calculated using Eq. (7)
as

$\mathbf{v}^{l}=\mathbf{v}^{s}-\frac{\mathbf{q}}{\rho^{l}}.$ (20)

##  4 Multi-time step algorithm

In multi-time scale coupling involving solid dynamic problems, different time
scales simultaneously exist. A multi-time step algorithm using explicit scheme
to match different time scale processes is introduced in this section. In this
paper, the slow process, e.g., fluid diffusion is integrated with larger time
step sizes, while the fast solid dynamics with smaller ones. With small time
step size, the solid dynamics evolves to a quasi-equilibrium state to update
velocity, position and other solid information. Further, in order to reduce
the stress relaxation time of solid dynamics, a damping scheme is applied to
accelerate the equilibrium process. For the following numerical simulations,
stretch loading or fluid diffusion is performed with a larger time step size,
while the dynamic stress relaxation coupled with a damping term is executed
with a smaller time step size.

###  4.1 Multi-time criteria

Since the explicit integration operator is conditionally stable, a time step
criterion ${\Delta t_{s}}$ in solid simulation is required when using explicit
scheme, stated as

$\Delta
t_{s}=0.6\min\left(\frac{h}{c_{s}+|\mathbf{v}_{s}|_{max}},\sqrt{\frac{h}{|\frac{\text{d}\mathbf{v}_{s}}{\text{d}t}|_{max}}}\right),$
(21)

where the artificial speed of sound of a solid structure
$c_{s}=\sqrt{K/\rho_{s}}$. In multi-time scale coupling problems, considering
that the solid dynamic relaxation process is comparatively fast, $\Delta
t_{s}$ is usually limited under a small value. In comparison, the time step
for internal diffusion evolution or stretching is allowed to be much larger.
For the tensile test simulation, we divide the stretching process into $N_{S}$
steps and the time step is

$\Delta t_{l}=\frac{T_{t}}{N_{S}},$ (22)

where $T_{t}$ is the entire process time of the tensile test, $\Delta t_{l}$
accordingly the time step for stretch loading. Similarly, for the fluid-
structure interaction, according to the Fick’s law, the maximum time step
allowed for explicit time stepping is characterized as [41]

$\Delta t_{d}=0.5\frac{h^{2}}{D},$ (23)

stating that the time step is mainly limited by the diffusivity constant $D$
and the kernel smoothing length $h$. To address the difference between these
time step sizes of different time scale processes, we present a multi-time
step algorithm to simulate these processes respectively with a iterative
scheme.

###  4.2 Iterative scheme

Figure 4.1: Flowchart of the iterative scheme in multi-time step algorithm.

Figure. 4.1 shows the iterative scheme of the proposed multi-time step
algorithm schematically. It can be seen this algorithm consists of two loops,
where the outer loop indicates that the entire dynamic progress is controlled
by the prescribed displacements or diffusion relaxation, which are executed
incrementally with a subscript $l$ denoting each increment. The inner loop
describes the solid dynamics evolution with a subscript $k$ signifying each
stress relaxation step. The loading or diffusion criterion $\Delta t_{l}$ or
$\Delta t_{d}$ controls the external force exerting or the fluid diffusion
process and ${\Delta t_{s}}$ determines the frequency of solid stress
relaxation. However, within one external loading time step $\Delta t_{l}$ or
diffusion time step $\Delta t_{d}$, the time integration of structure should
be computed as $k_{0}=[\frac{{\Delta t_{l/d}}}{\Delta t_{s}}]+1$ times. With a
limited ${\Delta t_{s}}$ and much larger $\Delta t_{l}$ and $\Delta t_{d}$,
$k_{0}$ is supposed to be very large and the computation of solid dynamics
will be trapped into a meaningless iteration, increasing the unnecessary
computation time.

Since once solid dynamics achieves the static state, the inner loop can be
finished to begin another outer loop. Therefore, in order to save computation
time, the inner loop is executed with a damping term to dissipate the kinetic
energy and accelerate the relaxation of the transient response. Solid
governing equations with extra damping can be solved a small number of times
$k$ until the kinetic energy is reduced to a sufficient small value $E_{k}$.
Specific criteria values of the kinetic energy are given in different cases.
After the equilibrium state of the solid deformation is achieved in the inner
steps, a new outer step begins and this procedure is performed once again
until the physical computation time ends.

###  4.3 Damping scheme

As we mentioned before, obtaining equilibrium for a dynamic system is
excessively time-consuming in SPH method with explicit time-stepping. To
address this issue, we apply a damping term into the stress relaxation to
dissipate the extra kinetic energy inside the system and accelerate the
convergence of stress relaxation process. Following Zhu et al. work [37], a
viscous damping term $\mathbf{f}^{v}$ is added in the solid momentum equation
as

$\frac{d\mathbf{v}}{dt}=\mathbf{f}^{s}+\mathbf{g}+\mathbf{f}^{v},$ (24)

Where $\mathbf{f}^{s}$ and $\mathbf{g}$ represents the surface and body
forces, the added damping term $\mathbf{f}^{v}$ can be discreized in the total
Lagrangian form as

$\mathbf{f}^{v}_{a}=\frac{\eta}{\rho_{a}}\nabla_{a}^{2}\mathbf{v}=\frac{2\eta}{m_{a}}\sum_{a}V_{a}V_{b}\mathbf{v}_{ab}\nabla_{a}^{0}W_{ab},$
(25)

where $\eta$ is the dynamic viscosity, given separately in different cases,
and usually it depends on the characteristic length scale of the problem and
materials parameters. $\mathbf{v}_{ab}=\mathbf{v}_{a}-\mathbf{v}_{b}$ denotes
the velocity difference between a particle pair $(a,b)$. This viscous force
can deduce the system oscillation caused by large velocity gradient and
eliminates the extra kinetic energy. Therefore, the solid stress is relaxed
much faster to a equilibrium state where the kinetic energy decreases below a
criterion value. Also, a pairwise splitting scheme is adopted to update the
velocity implicitly and locally, keeping the conservation of momentum in each
particle pair. More detailed information can be referred to Zhu’s work [37].

##  5 Numerical examples

In this section, several tests including the stretching-necking and the fluid
diffusion coupled solid deformation in two and three dimensions, are simulated
using the present method to show its accuracy and efficiency.

###  5.1 Necking of a two-dimensional bar

The standard tensile necking test simulation has been previously studied in
several papers [38, 42, 43, 44] with experimental and numerical results to
compare.

Table 5.1: Necking test simulation: physical material parameters. Parameters | Value
---|---
Shear modulus | 80.1938 Gpa
Bulk modulus | 164.21 Gpa
Initial flow stress | 450 MPa
Saturation flow stress | 715 MPa
Saturation exponent | 16.93
Linear hardening coefficient | 129.24 MPa

With a length of 53.334 mm and a width of 12.826 mm, the test sample is
stretched from the surface under an increasing (uniaxial) load. A reduction in
the width and thickness happens consistent with the elongation of this
specimen. A slight imperfection of this sample (1.8% reduction) is imposed
initially in the center part as shown in Figure 5.1 to trigger the necking
phenomenon.

Figure 5.1: 2D tensile necking: geometry and initial and boundary condition
setup.

The specimen is composed of a elastic deformation depicted by the Neo-Hookean
law and a plastic response by the nonlinear isotropic hardening law. The
material parameters are given in Table 5.1. A total stretching of 10 mm is
realized via a symmetric displacement boundary conditions. Here,
$dp=PH/50$=0.25652 mm. Three layers of particle are imposed with the
aforementioned boundary condition. Consistent with the experimental time
around 2 minutes, the physical time in this simulation is set to $t=100$s,
with stretching steps $N_{S}=$ 10000 the corresponding velocity $v=0.5\times
10^{-4}$ m/s. This is different from that in reference papers where the
velocity usually is increased to about 1 m/s to reduce the physical time to
$1.5\times 10^{-3}$ s. After each step of stretch loading, stress relaxation
coupled with damping is performed. The damping ratio is set to an experienced
value of $\eta=1.0\times 10^{4}$ based on the work of Zhu [37].

Figure 5.2 shows the deformation evolution colored by von Mise Strain at
different time instants. A clear necking pattern is observed in the center of
the specimen, which is consistent with that observed in both experimental and
other numerical works [42, 44, 43]. The specimen undergoes three distinct
stages: elastic strain, followed by uniform plastic strain, and finally
necking strain.

(a) 0s

(b) 50s

(c) 80s

(d) 100s

Figure 5.2: 2D tensile necking: the deformation colored by von Mise Strain at
different time instants.

Figure 5.3 plots the radius evolution of the central part where necking occurs
as a function of the imposed stretching displacement. It is compared with the
results from the reference Elguedj and Hughes [44] where different mesh
discretization and element types Q1, mixed Q1/P0, etc. are used to model this
test. As time progresses and the sample elongates, the radius displacement of
the central part increases linearly, while after necking occurs, it
experiences a rapid increase.

Figure 5.3: 2D tensile necking: the evolution of the radial displacement as a
function of the imposed vertical displacement of the central part.

Figure 5.4 depicts the evolution of the reaction force as time progresses.
After a short elastic response, represented by the initial straight line, the
specimen enters the stage of uniform plastic deformation with a smooth
increase of reaction force. During this stage, plastic deformation spreads
slowly and shows a homogeneous state throughout the specimen. Eventually, when
the boundary displacement reaches a certain value, necking occurs in the
central part, and the reaction force reaches its peak value. Subsequently, the
deformation changes to a mode where the plastic effect is concentrated in the
central zone, resulting in a decreasing reaction force, which is more obvious
in the following three dimensional case.

To determine when equilibrium is achieved, we monitor the kinetic energy
$E_{k}$ until it is damped below a threshold value derived from the elastic
energy $E_{e}$. Here, $E_{e}$ is calculated using the formula
$E_{e}=\frac{1}{2}F\Delta x$, where $F$ is the load force of $8000$ N deduced
from Figure 5.4, and $\Delta x$ is the stretching length of $10$ mm. To
investigate the effect of the kinetic energy threshold on the simulation
results, we conducted a series of stretching simulations with varying
criteria. Figure 5.5 plots the variation of the radius displacement and
reaction force for different kinetic energy criteria. Initially, we chose a
larger criterion value of $E_{k}=$ 5%$E_{e}$ and gradually decreased the
criterion. The results reveal that when $E_{k}$ is set to be 5%$E_{e}$, either
the radius displacement and loading force evolution is not smooth enough,
indicating that equilibrium is not achieved. This suggests that 5%$E_{e}$ is
too large as a criterion value. On the other hand, with too small criterion
values, unnecessary calculation steps are performed, increasing computation
time. The results demonstrate that for this 2D case, the appropriate kinetic
energy criterion value is 0.5%$E_{e}$.

Figure 5.4: 2D tensile necking: the evolution of the reaction force versus the
imposed vertical displacement.

(a)

(b)

Figure 5.5: 2D tensile necking: radius displacement (a) and the loading force
(b) convergence with different kinetic energy criteria.

During the simulation, the evolution of the kinetic energy after one
stretching at four different time instants, as evaluated by the elastic energy
$E_{e}$, is shown in Figure 5.6. As expected, due to the stretching force,
there is a kinetic energy fluctuation. After each stretching event, the
kinetic energy first increases, followed by a decrease to a certain criterion
value of 0.5%$E_{e}$, which is due to the damping effects. Throughout the
simulation process, stress relaxation occurs with viscous damping immediately
after each stretching. The relative kinetic energy at the end of each
stretching step approaches 0.5%$E_{e}$, showing that the equilibrium is
achieved.

Figure 5.6: 2D tensile necking: evolution of kinetic energy evaluated by the
elastic energy after one stretching at different time.

With a physical time in simulation $t=100$s, due to the time step size
limitation in explicit scheme, the performed stretching times $N_{S}$ and
stress relaxation times $N_{s}$ are supposed to be $N_{S}=N_{s}=t/\Delta
t_{s}=2.58\times 10^{9}$. With this multi-time criteria algorithm, we firstly
decrease the number of stretching time steps from $2.58\times 10^{9}$ to
$N_{S}=1.0\times 10^{4}$. Secondly, we decrease the stress relaxation times
from $2.58\times 10^{9}$ to $N_{s}=3.26\times 10^{5}$ by coupling the damping
term to accelerate the equilibrium obtaining. Table 5.2 lists the stress
relaxation iterations performed in straightforward and multi-time step
algorithms respectively and gives the quantitative efficiency of the present
algorithm compared against the straightforward one in terms of stretching
$N_{S}$ and stress relaxation iterations $N_{s}$ with the same total particle
number $N_{p}$. It is obvious that the proposed algorithm yields a drastic
reduction in computation time.

Table 5.2: 2D tensile necking: quantitative validation of the efficiency of this multi-time step algorithm. algorithm | $N_{p}$ | $N_{S}$ | $N_{s}$ | $N_{damping}$
---|---|---|---|---
straightforward algorithm | 10788 | $2.58e^{9}$ | $2.58e^{9}$ | -
multi-time step algorithm | 10788 | $1.0e^{4}$ | $3.26e^{5}$ | $3.26e^{5}$

###  5.2 Necking of a three-dimensional bar

Further, a three-dimensional necking analysis of a cylindrical bar is carried
out, which has been studied by Simo and Armero [45, 38], de Souza Neto et al.
[42], Elguedj and Hughes [44]. The same geometry of radius 6.413 mm and length
53.334 mm with a slight reduction (1.8%) in the center of the bar as in the
previous 2D case is considered. Loading is imposed using displacement control,
with a total vertical displacement of 7 mm applied on both the top and bottom
surface of the bar. The same material properties in Table 5.1 and elastic-
plastic response as that applied in previous two-dimensional case are employed
herein. In this work, initial particle spacing $dp$ = 0.3 mm with a total
particle number almost $N_{p}=2.5e^{5}$. With physical time $t=100$s and
stretching steps $N_{S}=10000$, the corresponding velocity is $0.7\times
10^{-4}$ m/s, which allows problem to be simulated in a real stretching rate.
The damping ratio used here is $\eta=1.0\times 10^{4}$.

(a) 0s

(b) 50s

(c) 80s

(d) 100s

Figure 5.7: 3D tensile necking: the deformation colored by von Mise Strain at
different time instants (top side view).

(a) 0s

(b) 50s

(c) 80s

(d) 100s

Figure 5.8: 3D tensile necking: the deformation colored by von Mise Strain at
different time instants (front view with half the specimen).

Contour plots of the von Mise Strain at different time instants from different
views are shown in Figure 5.7-5.9. The last plots depict the deformed shape of
the specimen at the final stage of the simulation, indicating the occurrence
of a necking in the center of the specimen. Based on these figures, we can
deduce the deformation evolution of this specimen: initially, the boundary
conditions enabled the specimen to maintain an uniform elastic response in the
short stage of loading history; subsequently, in the post-peak regime, a
diffuse necking mode emerged, which eventually led to the formation of shear
bands at high strain levels. These bands accumulated plastic deformations,
ultimately leading to the final necking even failure of the specimen. The
evolution of this pattern is well-reproduced by the force and deformation data
presented in Figures 5.10 and 5.11, which agrees well with experimental
findings.

Specially, Figure 5.10 collected the numerical data of the radius
displacement, normalized by the initial radius, versus the length
displacement, normalized by the initial bar length, and the results are
compared with experimental data and solutions reported in other works [46, 47,
34, 38]. The results obtained with the present method seems to be in good
agreement with the experimental data of material 2501R. The load displacement
curve obtained from this simulation is presented in Figure. 5.11 along with
the results from other numerical works. Some differences can be observed in
the later stages of the deformation and force variation curves, which may be
attributed to the hourglass mode and the limitations of particle
discretization. From the findings of de Saracibar [48], different mesh
discretization and element types can result in distinct stabilization and
locking patterns, as evidenced by the disparate peak and final values of force
in load-displacement curves obtained using different mesh discretization
strategies, including uniform and non-uniform mesh distributions in the
necking area. A comprehensive comparison of these results, as well as a clear
explanation, is provided in de Saracibar’s work. Although this hourglass
phenomenon depicted in Figure 5.10 cannot be fully removed in this study, the
same pattern is observed in prior research that utilized standard 3D elements
in FEM [42]. The current paper does not aim to address this issue
comprehensively, but rather suggests that further improvements will be
implemented in the future research.

(a) 0s

(b) 50s

(c) 80s

(d) 100s

Figure 5.9: 3D tensile necking: the deformation colored by von Mise Strain at
different time instants (top view with quarter the specimen).  Figure 5.10: 3D
tensile necking: the evolution of radial displacement of the central part
compared with the reference [44, 49, 42].  Figure 5.11: 3D tensile necking:
the overall evolution of the reaction force versus the imposed vertical
displacement compared with the reference [44, 49, 42].

As for determining the static state achieving, similar with previous 2D case,
the kinetic energy criterion is derived from the elastic energy
$E_{e}=\frac{1}{2}F\Delta x$, while the load force $F=80000N$, consistent with
the value in Figure 5.11, and the stretching length $\Delta x=0.014$ m. Using
the same method applied in the 2D case, we find that the results are converged
when the kinetic energy is reduced to 0.5% of the elastic energy. Finally, to
check the efficiency of the present algorithm, we computed the relaxation
iterations and the results are summarized in Table 5.3. Evident reduction is
obtained in stress relaxation times, saving the computation time
significantly.

Table 5.3: 3D tensile necking: quantitative validation of the efficiency of this multi-time step algorithm. algorithm | $N_{p}$ | $N_{S}$ | $N_{s}$ | $N_{d}$
---|---|---|---|---
straightforward algorithm | 250852 | $2.21e^{9}$ | $2.21e^{9}$ | -
multi-time step algorithm | 250852 | $1.0e^{4}$ | $3.21e^{5}$ | $3.21e^{5}$

###  5.3 Two-dimensional fluid-structure interaction

In this section, we perform a two-dimensional simulation of fluid diffusion
coupling with porous solid deformation and the model is described in Section
2.2, to verify the efficiency of the presented method. As Figure 5.12 shows, a
thin porous beam with a length of $L$ = 10.0 mm and width of $W$ = 0.125 mm is
considered, with the left and right sides being constrained to prevent any
curling or movement. The simulation starts with a fluid droplet contacting the
center part of the beam with a length of 0.3$L$, and this contact continues
for 10 seconds while the total physical time is 100 seconds. Given the thin
nature of the beam, we assume all pores in the upper half part are filled with
fluid initially. As we stated before, the relationship between fluid
saturation $\widetilde{a}$ and solid porosity ${a}$ is
$0\leq\widetilde{a}\leq{a}<1$. For this 2D and 3D cases discussed later, we
assume a solid porosity of ${a}=0.4$, meaning that the fluid saturation
$\widetilde{a}$ in the central part($0.5W\times 0.3L$) is constrained to
$\widetilde{a}={a}=0.4$ for the initial 10 seconds, while in other regions
$\widetilde{a}_{0}=0.0$.

Figure 5.12: 2D fluid-structure interaction: physical configuration of the
thin porous beam.

In accordance with the experimental setup, the solid material is considered as
a porous and elastic Nafion membrane, with water serving as the fluid. The
physical properties and material parameters of this membrane are listed in
Table 5.4. The pressure coefficient C has been calibrated to fit the
experimentally measured flexure curves, while other parameters are obtained
from previous research papers [50, 51].

Table 5.4: Fluid-structure interaction: physical material parameters value of
Nafion film. Data estimated from Motupally and Goswami [50, 51].

Parameters |  $\rho$ ${\rm(kg/m^{3})}$ | K ${\rm(m^{2}/s)}$ | Pressure coefficient C ${\rm(Pa)}$ | Young modulus ${\rm(Pa)}$ | Poisson ratio
---|---|---|---|---|---
Value | 2000 | $1.0e^{-10}$ | $3.0e^{6}$ | $8.242e^{6}$ | 0.2631

In the simulation, eight particles are placed in the vertical direction, with
a particle spacing of $dp_{y}=W/8=1.5625\times 10^{-2}$ mm. However, due to
the high aspect ratio of the beam, using the same particle spacing $dp$ in
horizontal $x$ and vertical $y$ directions would require a large number of
particles, thus increasing the computation time. To address this issue, an
anisotropic kernel algorithm is employed, with an anisotropic ratio of 4.0,
meaning $dp_{x}=4dp_{y}=0.0625$ mm. In this simulation, an experienced damping
ratio of $\eta=1.0e^{3}$ in the damping term is utilized.

(a) 0s

(b) 10s

(c) 50s

(d) 100s

Figure 5.13: 2D fluid-structure interaction: the deformation colored by fluid
saturation at different time instants.  Figure 5.14: 2D fluid-structure
interaction: bending amplitude of the beam at different time instants.

With the conditions given above, the simulation produces a deformed
configuration colored by fluid saturation, as shown in Figure 5.13. Initially,
the presence of a water droplet in the upper central region generates a fluid
pressure, as explained in Eq. 43, leading to a localized bending in the
central region. As time progresses, the saturation difference drives water
diffusing continuously, and the total water amount within the porous solid
increases, causing a rising flexure. This is also depicted in Figure 5.14,
which records the vertical position $y$ versus the horizontal $x$ position of
the beam at different time instants. After the contact finishes, no more water
is added into the beam, and the central water flows slowly into the side
areas. Clearly, the fluid saturation shows a smooth transition from the center
to the surrounding area in Figure 5.13. Accordingly, a more uniform pressure
distribution is developing, resulting in a more smooth flexure of the beam as
shown in Figure 5.14 in the later period.

Figure 5.15: 2D fluid-structure interaction: bending amplitude convergence
with different density kinetic energy criteria. Figure 5.16: 2D fluid-
structure interaction: the density kinetic energy variation within the
diffusion period when $t=20$s valuated by the water pressure $p^{l}$.

For determining the density kinetic energy criterion $E_{k}$, we use the
pressure from water, $p^{l}$, stated in Eq. 43, as the reference since the
fluid pressure induces the beam swelling. To evaluate the effect of the
relative density kinetic energy threshold on the simulation results, a series
of simulations are conducted using various criteria $E_{k}$. The time
evolution of the bending amplitude with different kinetic energy criteria is
presented in Figure 5.15. With a relatively large criterion value of $E_{k}$ =
5%$p^{l}$, it is observed that the equilibrium state is not achieved and the
energy is not fully eliminated with a relatively light deformation. On the
other hand, using a very small criterion value leads to unnecessary
calculation steps, increasing computation time. Therefore, it can be concluded
that the appropriate density kinetic energy criterion value for this 2D case
is 0.05%$p^{l}$.

Referring to Figure 5.16, the evolution of the density kinetic energy within
the diffusion period when $t=20$s, evaluated by the water pressure $p^{l}$, is
presented. Due to the water pressure, the density kinetic energy firstly
experiences a peak after one diffusion performance, then followed by a
decrease to a certain criterion value of 0.05% $p^{l}$ we set before, which is
attributed to the damping effects. Throughout the simulation process, the
stress relaxation takes place accompanying with viscous damping immediately
after each diffusion relaxation event. The relative density kinetic energy at
the end of each diffusion step approaches 0.05%$p^{l}$, indicating that the
velocity almost vanishes. This signifies that equilibrium is achieved at the
end of each diffusion time step.

The efficiency of the proposed approach is demonstrated through Table 5.5,
which presents a quantitative comparison of the algorithm against the
straightforward approach in terms of diffusion and stress relaxation
iterations $N_{D}$, $N_{s}$ with a total particle number $N_{p}$. The results
reveal a great reduction in computation iterations, thus demonstrating the
significant improvement in efficiency achieved by the proposed approach.

Table 5.5: 2D fluid-structure interaction: quantitative validation of the efficiency of this multi-time step algorithm. algorithm | $N_{p}$ | $N_{D}$ | $N_{s}$ | $N_{d}$
---|---|---|---|---
straightforward algorithm | 1336 | 1.58$e^{7}$ | 1.58$e^{7}$ | -
multi-time step algorithm | 1336 | 125 | 2.76$e^{5}$ | 2.76$e^{5}$

###  5.4 Three-dimensional fluid-structure interaction

Next, we consider the fluid diffusion coupling swelling in a three-dimensional
film, specifically the diffusion of water within a porous Nafion membrane.
This system has been previously studied numerically by Zhao [24] and
experimentally by Goswami [51]. This reference thin porous body is in the form
of a polymer film with a x-y plane of dimensions $L_{x}=10.0$ mm, $L_{y}=10.0$
mm and a height of $L_{z}=0.125$ mm. Four boundary sides are constrained to
prevent any curling or movement. The physical parameters are taken to be the
same as those listed in Table 5.4. The initial conditions are similar to those
used in the two-dimensional case. The central square part of the membrane in
contact with water occupies a region of dimensions $0.3L_{x}\times
0.3L_{y}\times 0.5L_{z}$, and this contact lasts for 450 seconds, while the
total physical time is 2500 seconds. No fluid is allowed to diffuse out from
the membrane. The fluid saturation $\widetilde{a}$ in the central square part
is constrained to $\widetilde{a}={a}=0.4$ for the initial 450 seconds, while
in other regions $\widetilde{a}_{0}=0.0$. Similar with that in the previous
two-dimension case, an anisotropic kernel algorithm to reduce the total
particle number evolving in this membrane simulation is used. Specifically, 8
particles are set in the vertical $z$ direction, meaning that the particle
spacing $dp_{z}=W/8=1.5625\times 10^{-2}$ mm. Here, the anisotropic ratio is
8.0, meaning the $dp_{x}=dp_{y}=8dp_{z}=0.125$ mm. In the stress relaxation
process of the simulation, the experienced damping ratio is set to
$\eta=1.0e^{4}$. In terms of the convergence study of density kinetic energy
criteria, by using the same method with that in 2D case, the 3D case has a
converged criterion value of $E_{d}$ = 0.1% $p^{l}$.

(a) t = 450s, front view

(b) t = 1500s, front view

(c) t = 450s, top side view

(d) t = 1500s, top side view

(e) t = 450s, top view

(f) t = 1500s, top view

Figure 5.17: 3D fluid-structure interaction: the deformation colored by water
saturation at different time instants.

In order to provide a more accurate representation of the experiment, the
evaporation process is taken into consideration, i.e., the water loses as time
progresses. Deformation flexure occurs during the initial period, and later as
the mass of fluid loses from the membrane, it eventually returns back to the
original shape.

Figure 5.17 shows the membrane deformation colored by water saturation at
different time instants. In the first 450 seconds, water amount continues to
increase as time progresses, leading to a rising flexure as depicted in Figure
5.18, which records the time history of the height $z$ of the central point.
Once the contact period finishes, no further water is added into the beam, and
the central water flows slowly into the side areas. At the same time, water
evaporates from the membrane, resulting in a rapid decrease of water pressure
and a corresponding decrease of the flexure, as shown by the blue line in
Figure 5.18 beyond 450 seconds. Figure 5.18 also includes corresponding data
points measured experimentally by Goswami [51] and results from other
numerical models for the swelling degree of the very center point versus
different time instants. Clearly, the present numerical simulation results
exhibit good agreement with experimental results in terms of the deformation
amplitude pattern, reproducing the increasing flexure during the water contact
period and the subsequent decrease after the contact finishes, consistent with
the saturation variation.

Figure 5.18: 3D fluid-structure interaction: bending amplitude of the center
point compared with experimental data and results from other numerical models.
Figure 5.19: 3D fluid-structure interaction: the stress iterations history
during the whole simulation.

Drawing from the previous discussion, the optimal large outer time step is
determined by the diffusion constant and the smoothing length, while the small
inner time step is dictated by the material properties of the solid. Ideally,
the outer time step allowed in principle is hundreds or thousands of times
larger than the inner time step size allowed. However, in the standard
explicit algorithm, the time step is limited to the smaller one, resulting in
the execution of numerous stress relaxation steps and consuming a substantial
amount of time. In the presented method, first, diffusion is performed with
the larger time step, while stress relaxation is executed multiple times with
damping effects until a kinetic energy is reached. Our approach saves time in
two ways. Firstly, the number of diffusion relaxation times is reduced since
multi-time step algorithm allows diffusion to be performed with its own time
step as the outer loop. Secondly, once the kinetic energy criterion is
satisfied, we consider the equilibrium achieved, and the inner loop is halted
accordingly, avoiding unnecessary stress relaxation calculations. Figure 5.19
indicates the stress iterations $N_{s}$ during this 3D simulation. There is an
increase in the initial 450 seconds when the fluid is in contact with the
film, and then a slower increase in the later stages. Table 5.6 presents the
quantitative efficiency of our new algorithm compared to the straightforward
one, by listing the diffusion iterations $N_{D}$ and $N_{s}$ separately. As
shown in the table, both two iterations are obviously reduced, representing a
significant improvement in saving computation time.

Table 5.6: 3D fluid-structure interaction: quantitative validation of the efficiency of this multi-time step algorithm. algorithm | membrane | $N_{D}$ | $N_{s}$ | $N_{d}$
---|---|---|---|---
straightforward algorithm | 60552 | 1.5$e^{10}$ | 1.5$e^{10}$ | -
multi-time step algorithm | 60552 | 1.25$e^{5}$ | 2.89$e^{6}$ | 2.89$e^{6}$

##  6 Conclusion

This paper proposed an approach employing a multi-time step algorithm to solve
multi-time coupling problem involving solid dynamics. In this algorithm, the
explicit scheme in time integration is used to simplify the equation system
solving. Inner and outer loops with different time step sizes are carried out
to match different time scale process. Another crucial feature of this
algorithm is the utilization of a kinetic energy criterion to ascertain the
attainment of equilibrium of solid dynamics and a damping term to accelerate
this equilibrium attainment process, thereby enabling the earlier termination
of the inner loop of solid stress relaxation and avoiding redundant
computations. Two types of multi-time coupling problem, including a nonlinear
hardening bar stretching and a fluid diffusion in porous media coupling solid
deformation are simulated to test the performance of this algorithm. Results
demonstrate the accuracy and a significant decrease in computation time.
Further, the application of this algorithm in practical fluid diffusion
coupling hydrogel deformation paves the way for simulating complex multi-
physics problems of multi-time scales in the field of complex chemistry
reaction.

Authorship contribution statement

Xiaojing Tang made the methodology, designed the research, developed code and
tested the present library components, performed the visualization and
validation, and wrote the original draft of the manuscript. Dong Wu
investigated the topic, made the methodology, developed code and tested the
present library components, conducted the formal analysis, modified the draft.
Zhentong Wang developed code and tested the present library components, and
revised the manuscript. Oskar Haidn and Xiangyu Hu made the conceptualization,
supervised and administered the project, and revised the manuscript.

Statements and Declarations

The authors have no known competing financial interests or personal
relationships that could have appeared to influence the work reported in this
paper.

Acknowledgments

Xiaojing Tang was partially supported by the China Scholarship Council (Grant
No. 201906120034). Dong Wu was partially supported by the China Scholarship
Council (Grant No. 20190613018). Xiangyu Hu would like to express his
gratitude to Deutsche Forschungsge meinschaft (DFG) for their sponsorship of
this research (Grant No. DFG HU1527/12-4).

## References

  * [1] L. B. Lucy, “A numerical approach to the testing of the fission hypothesis,” _The Astronomical Journal_ , vol. 82, pp. 1013–1024, 1977.
  * [2] R. A. Gingold and J. J. Monaghan, “Smoothed particle hydrodynamics: Theory and application to non-spherical stars,” _Monthly Notices of the Royal Astronomical Society_ , vol. 181, no. 3, pp. 375–389, 1977.
  * [3] J. J. Monaghan, “Simulating free surface flows with SPH,” _Journal of Computational Physics_ , vol. 110, no. 2, pp. 399–406, 1994.
  * [4] X. Y. Hu and N. A. Adams, “A multi-phase SPH method for macroscopic and mesoscopic flows,” _Journal of Computational Physics_ , vol. 213, no. 2, pp. 844–861, 2006.
  * [5] S. Shao, C. Ji, D. I. Graham, D. E. Reeve, P. W. James, and A. J. Chadwick, “Simulation of wave overtopping by an incompressible SPH model,” _Coastal engineering_ , vol. 53, no. 9, pp. 723–735, 2006.
  * [6] C. Zhang, G. Xiang, B. Wang, X. Hu, and N. Adams, “A weakly compressible SPH method with WENO reconstruction,” _Journal of Computational Physics_ , vol. 392, pp. 1–18, 2019.
  * [7] L. D. Libersky and A. G. Petschek, “Smooth particle hydrodynamics with strength of materials,” in _Advances in the Free-Lagrange Method Including Contributions on Adaptive Gridding and the Smooth Particle Hydrodynamics Method_. Springer, 1991, pp. 248–257.
  * [8] W. Benz and E. Asphaug, “Simulations of brittle solids using smooth particle hydrodynamics,” _Computer Physics Communications_ , vol. 87, no. 1-2, pp. 253–265, 1995.
  * [9] J. J. Monaghan, “SPH without a tensile instability,” _Journal of Computational Physics_ , vol. 159, no. 2, pp. 290–311, 2000.
  * [10] P. Randles and L. D. Libersky, “Smoothed particle hydrodynamics: Some recent improvements and applications,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 139, no. 1-4, pp. 375–408, 1996.
  * [11] C. Zhang, Y. Zhu, Y. Yu, D. Wu, M. Rezavand, S. Shao, and X. Hu, “An artificial damping method for total Lagrangian SPH method with application in biomechanics,” _Engineering Analysis with Boundary Elements_ , vol. 143, pp. 1–13, 2022.
  * [12] C. Antoci, M. Gallati, and S. Sibilla, “Numerical simulation of fluid–structure interaction by SPH,” _Computers and Structures_ , vol. 85, no. 11-14, pp. 879–890, 2007.
  * [13] L. Han and X. Hu, “SPH modeling of fluid-structure interaction,” _Journal of Hydrodynamics_ , vol. 30, no. 1, pp. 62–69, 2018.
  * [14] C. Zhang, M. Rezavand, and X. Hu, “A multi-resolution SPH method for fluid-structure interactions,” _Journal of Computational Physics_ , vol. 429, p. 110028, 2021.
  * [15] M. Liu and G. Liu, “Smoothed particle hydrodynamics (SPH): an overview and recent developments,” _Archives of computational methods in engineering_ , vol. 17, no. 1, pp. 25–76, 2010.
  * [16] J. J. Monaghan, “Smoothed particle hydrodynamics and its diverse applications,” _Annual Review of Fluid Mechanics_ , vol. 44, pp. 323–346, 2012.
  * [17] M. Luo, A. Khayyer, and P. Lin, “Particle methods in ocean and coastal engineering,” _Applied Ocean Research_ , vol. 114, p. 102734, 2021.
  * [18] C. Zhang, Y. Zhu, D. Wu, and X. Hu, “Review on smoothed particle hydrodynamics: Methodology development and recent achievement,” _arXiv preprint arXiv:2205.03074_ , 2022.
  * [19] H. Gotoh and A. Khayyer, “On the state-of-the-art of particle methods for coastal and ocean engineering,” _Coastal Engineering Journal_ , vol. 60, no. 1, pp. 79–103, 2018.
  * [20] J. U. Brackbill and B. I. Cohen, _Multiple time scales_. Academic Press, 2014, vol. 3.
  * [21] D. Knoll, L. Chacon, L. Margolin, and V. Mousseau, “On balanced approximations for time integration of multiple time scale systems,” _Journal of Computational Physics_ , vol. 185, no. 2, pp. 583–611, 2003.
  * [22] A. Prior, “Applications of implicit and explicit finite element techniques to metal forming,” _Journal of Materials Processing Technology_ , vol. 45, no. 1-4, pp. 649–656, 1994.
  * [23] E. Gavalas, I. Pressas, and S. Papaefthymiou, “Mesh sensitivity analysis on implicit and explicit method for rolling simulation,” _International Journal of Structural Integrity_ , 2018.
  * [24] Q. Zhao and P. Papadopoulos, “Modeling and simulation of liquid diffusion through a porous finitely elastic solid,” _Computational Mechanics_ , vol. 52, no. 3, pp. 553–562, 2013.
  * [25] D. Gaston, C. Newman, G. Hansen, and D. Lebrun-Grandie, “Moose: A parallel computational framework for coupled systems of nonlinear equations,” _Nuclear Engineering and Design_ , vol. 239, no. 10, pp. 1768–1778, 2009\.
  * [26] E. Gavalas and S. Papaefthymiou, “Brief overview on the application of finite element method on flat rolling processes,” _Technology_ , vol. 105, no. 1-2, pp. 110–118.
  * [27] J. Sun, K. Lee, and H. Lee, “Comparison of implicit and explicit finite element methods for dynamic problems,” _Journal of materials processing technology_ , vol. 105, no. 1-2, pp. 110–118, 2000.
  * [28] M. Rezaiee-Pajand and J. Alamatian, “The dynamic relaxation method using new formulation for fictitious mass and damping,” _Structural engineering and mechanics_ , vol. 34, no. 1, p. 109, 2010.
  * [29] R. Yaghmaie and S. Ghosh, “Multi-time scale based modeling of piezoelectric materials coupling transient electrical and dynamic fields with finite deformation damage,” _International Journal of Solids and Structures_ , vol. 202, pp. 338–355, 2020.
  * [30] J. C. Ragusa and V. S. Mahadevan, “Consistent and accurate schemes for coupled neutronics thermal-hydraulics reactor analysis,” _Nuclear Engineering and Design_ , vol. 239, no. 3, pp. 566–579, 2009.
  * [31] L. Beuth, “Formulation and application of a quasi-static material point method,” 2012.
  * [32] F. Harewood and P. McHugh, “Comparison of the implicit and explicit finite element methods using crystal plasticity,” _Computational Materials Science_ , vol. 39, no. 2, pp. 481–494, 2007.
  * [33] S. Doll, K. Schweizerhof, R. Hauptmann, and C. Freischläger, “On volumetric locking of low-order solid and solid-shell elements for finite elastoviscoplastic deformations and selective reduced integration,” _Engineering Computations_ , vol. 17, no. 7, pp. 874–902, 2000.
  * [34] P. R. R. de Campos, A. J. Gil, C. H. Lee, M. Giacomini, and J. Bonet, “A new updated reference lagrangian smooth particle hydrodynamics algorithm for isothermal elasticity and elasto-plasticity,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 392, p. 114680, 2022.
  * [35] Y. G. Rao, M. Nataraj, and P. Srinivas, “Explicit dynamic analysis of tensional & torsional propagations on composite material with dog bone shaped testing specimen,” _International Journal for Modern Trends in Science and Technology ISSN_ , pp. 2455–3778, 2019.
  * [36] M. Yılmaz, Ö. Anıl, B. Alyavuz, and E. Kantar, “Load displacement behavior of concrete beam under monotonic static and low velocity impact load,” _International Journal of Civil Engineering_ , vol. 12, no. 4, pp. 488–503, 2014.
  * [37] Y. Zhu, C. Zhang, and X. Hu, “A dynamic relaxation method with operator splitting and random-choice strategy for SPH,” _Journal of Computational Physics_ , vol. 458, p. 111105, 2022.
  * [38] J. C. Simo and T. J. Hughes, _Computational inelasticity_. Springer Science & Business Media, 2006, vol. 7.
  * [39] Y. Yue, B. Smith, C. Batty, C. Zheng, and E. Grinspun, “Continuum foam: A material point method for shear-dependent flows,” _ACM Transactions on Graphics (TOG)_ , vol. 34, no. 5, pp. 1–20, 2015.
  * [40] R. Vignjevic, J. R. Reveles, and J. Campbell, “SPH in a total Lagrangian formalism,” _CMC-Tech Science Press-_ , vol. 4, no. 3, p. 181, 2006.
  * [41] P. W. Cleary and J. J. Monaghan, “Conduction modelling using smoothed particle hydrodynamics,” _Journal of Computational Physics_ , vol. 148, no. 1, pp. 227–264, 1999.
  * [42] E. D. S. Neto, F. A. Pires, and D. Owen, “F-bar-based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. part I: Formulation and benchmarking,” _International Journal for Numerical Methods in Engineering_ , vol. 62, no. 3, pp. 353–383, 2005.
  * [43] C. A. de Saracibar, M. Chiumenti, Q. Valverde, and M. Cervera, “On the orthogonal subgrid scale pressure stabilization of finite deformation J2 plasticity,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 195, no. 9-12, pp. 1224–1251, 2006.
  * [44] T. Elguedj and T. J. Hughes, “Isogeometric analysis of nearly incompressible large strain plasticity,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 268, pp. 388–416, 2014.
  * [45] J.-C. Simo and F. Armero, “Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes,” _International Journal for Numerical Methods in Engineering_ , vol. 33, no. 7, pp. 1413–1449, 1992.
  * [46] F. Armero and E. Love, “An arbitrary lagrangian–eulerian finite element method for finite strain plasticity,” _International Journal for Numerical Methods in Engineering_ , vol. 57, no. 4, pp. 471–508, 2003.
  * [47] J. C. Simo, “A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. part II: Computational aspects,” _Computer methods in applied mechanics and engineering_ , vol. 68, no. 1, pp. 1–31, 1988.
  * [48] C. A. de Saracibar, M. Chiumenti, Q. Valverde, and M. Cervera, “On the orthogonal subgrid scale pressure stabilization of small and finite deformation J2 plasticity,” _Monograph Series on Computational Methods in Forming Processes_ , vol. 2, 2004.
  * [49] A. Rodríguez-Ferran, A. Pérez-Foguet, and A. Huerta, “Arbitrary lagrangian–eulerian (ALE) formulation for hyperelastoplasticity,” _International Journal for Numerical Methods in Engineering_ , vol. 53, no. 8, pp. 1831–1851, 2002.
  * [50] S. Motupally, A. J. Becker, and J. W. Weidner, “Diffusion of water in Nafion 115 membranes,” _Journal of The Electrochemical Society_ , vol. 147, no. 9, p. 3171, 2000.
  * [51] S. Goswami, S. Klaus, and J. Benziger, “Wetting and absorption of water drops on Nafion films,” _Langmuir_ , vol. 24, no. 16, pp. 8627–8633, 2008.
  * [52] R. v. Mises, “Mechanik der festen körper im plastisch-deformablen zustand,” _Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse_ , vol. 1913, pp. 582–592, 1913\.
  * [53] M. L. Wilkins, “Calculation of elastic-plastic flow,” California Univ Livermore Radiation Lab, Tech. Rep., 1963.
  * [54] J.-P. Ponthot, “Unified stress update algorithms for the numerical simulation of large deformation elasto-plastic and elasto-viscoplastic processes,” _International Journal of Plasticity_ , vol. 18, no. 1, pp. 91–126, 2002\.
  * [55] F. Dunne and N. Petrinic, _Introduction to computational plasticity_. OUP Oxford, 2005.
  * [56] M.-H. Yu, _Generalized plasticity_. Springer Science & Business Media, 2006.
  * [57] D. Gawin, P. Baggio, and B. A. Schrefler, “Coupled heat, water and gas flow in deformable porous media,” _International Journal for numerical methods in fluids_ , vol. 20, no. 8-9, pp. 969–987, 1995.
  * [58] J. Korsawe, G. Starke, W. Wang, and O. Kolditz, “Finite element analysis of poro-elastic consolidation in porous media: Standard and mixed approaches,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 195, no. 9-12, pp. 1096–1115, 2006.
  * [59] J. Ghaboussi and E. L. Wilson, “Flow of compressible fluid in porous elastic media,” _International Journal for Numerical Methods in Engineering_ , vol. 5, no. 3, pp. 419–442, 1973.
  * [60] R. J. Atkin and R. Craine, “Continuum theories of mixtures: Basic theory and historical development,” _The Quarterly Journal of Mechanics and Applied Mathematics_ , vol. 29, no. 2, pp. 209–244, 1976.

##  Appendix A Plasticity theory and nonlinear hardening plastic model

In this appendix, we present the J2 plasticity theory coupling with a
hardening elastic-plastic model to determine the plastic deformation. In A1 we
describe the multiplicative decomposition technique for the material
deformation. A2 presents the constitutive relation in this model. To describe
the strain-stress evolution, the flow rule and a hardening plastic model is
stated in A3. Then a return mapping algorithm is given in A4 to explain the
time integration.

###  A.1 Multiplicative decomposition technique

To describe the elastoplastic model, we adopted the flow plasticity theory
where total strain can be multiplicatively decomposed into an elastic part and
a plastic part [38, 39]. Using this technique, $\mathbf{F}$ can be written as
the product of its elastic volumetric part $\mathbf{F}^{e}$ and plastic
deviatoric part $\mathbf{F}^{p}$:

$\mathbf{F}=\mathbf{F}^{e}\mathbf{F}^{p}.$ (26)

Similarly, $\mathbf{b}^{e}$, the elastic part of the left Cauchy-Green tensor
$\mathbf{b}=\mathbf{F}\mathbf{F}^{T}$, is defined as
$\mathbf{b}^{e}=\mathbf{F}^{e}\mathbf{F}^{e\,T}$. When strains are within the
elastic range, $\mathbf{F}=\mathbf{F}^{e}$ and $\mathbf{b}=\mathbf{b}^{e}$.
For plasticity analysis, the plastic Lagrangian tensor $\mathbf{C}^{p}$ is
introduced as

$\mathbf{C}^{p}=\mathbf{F}^{p\,T}\mathbf{F}^{p}.$ (27)

The relationship between $\mathbf{b}^{e}$ and $\mathbf{C}^{p}$ can be
described as

$\mathbf{b}^{e}=\mathbf{F}\mathbf{C}^{p\,-1}\mathbf{F}^{T},$ (28)

which is used in the subsequent hardening plastic model. Additionally, to
adhere to the volume preserving assumption in plasticity, we assume that the
determinant of the plastic deformation part, det($\mathbf{F}^{p}$) $=1$.

###  A.2 Constitutive relation

According to the theoretical framework proposed by Simo and Hughes [38], with
an isotropic stress response assumption, the elastoplastic constitutive model
incorporates a nonlinear elastic strain energy function which is decomposed
into volumetric and deviatoric parts:

$W=W_{v}(J)+W_{s}(\overline{\mathbf{b}}^{e}),$ (29)

where $\overline{\mathbf{b}}^{e}$ is the volume-preserving left-Cauchy Green
tensor. The volumetric part weighted by the bulk modulus $K$, is given by

$W_{v}(J)=\frac{1}{2}K[\frac{1}{2}(J^{2}-1)-\ln J].$ (30)

The deviatoric part related to the shear modulus $\mu$, is obtained by

$W_{s}({\overline{\mathbf{b}}^{e}})=\frac{1}{2}\mu[\operatorname{Tr}(\overline{\mathbf{b}}^{e})-D].$
(31)

Here, $D=$ {1,2,3} depends on the dimension of the problem. With $E$ denoting
Young’s modulus and $\nu$ the Poisson ratio, the bulk and shear moduli are
interconnected through the relationship:

$E=2\mu\left(1+\nu\right)=3K\left(1-2\nu\right).$ (32)

With the energy function Eq. (29) in hand, the Kirchhoff stress tensor, which
characterizes the stress response, can be expressed as

$\bm{\tau}=\frac{\partial
W}{\partial\mathbf{F}^{e}}\mathbf{F}^{eT}=\frac{\kappa}{2}\left(J^{2}-1\right)\mathbf{I}+\mu\operatorname{dev}(\overline{\mathbf{b}}^{e}),$
(33)

where the two parts account for the volumetric and shear stresses.

Note that in the equations above, the expression

$\overline{\mathbf{T}}=[\operatorname{det}(\mathbf{T})]^{-1/3}\mathbf{T}$ (34)

indicates the volume preserving treatment of a tensor $\mathbf{T}$.
Additionally,

$\operatorname{dev}(\mathbf{T})=\mathbf{T}-\frac{\operatorname{Tr}(\mathbf{T})}{3}\mathbf{I}$
(35)

represents the trace free part of the tensor $\mathbf{T}$, i.e.,
$\operatorname{Tr}[\operatorname{dev}(\mathbf{T})]=0$ and
$\operatorname{dev}(\mathbf{I})=0$.

###  A.3 Flow rule and hardening plasticity model

With the flow plasticity theory, a flow rule is needed to determine the
orientation and magnitude of plastic deformation. In this paper, the classical
$J_{2}$ flow theory, also known as the Mises–Huber yield condition proposed by
Hube and von Mises [52], is used to model plasticity stress-strain evolution.
This theory states that the plastic behavior is governed by the deviatoric
part of the Kirchhoff stress tensor $\bm{\tau}$, represented by the second
term $\mu\operatorname{dev}(\overline{\mathbf{b}}^{e})$ in Eq. (33). For
simplicity, we define
$\mathbf{s}=\operatorname{dev}(\bm{\tau})=\mu\operatorname{dev}(\overline{\mathbf{b}}^{e})$.
With a Frobenius norm $\|.\|_{F}$, the magnitude scalar $s=\|\mathbf{s}\|_{F}$
is used to compare with the yield criterion to determine the onset of
plasticity, and the normalized tensor of ${\mathbf{s}}$ is given by
$\hat{\mathbf{s}}=\mathbf{s}/s$.

In mechanical engineering, isotropic work hardening plastic behavior is
commonly observed. To incorporate this behavior, a scalar yield function
$f(\bm{\tau},\alpha)$ that depends on the hardening function $k(\alpha)$ is
introduced, where $\alpha$ represents the equivalent plastic strain. The yield
function is formulated as

$f(\bm{\tau},\alpha)=\|\operatorname{dev}(\bm{\tau})\|_{F}-\sqrt{\frac{2}{3}}k(\alpha)=s-\sqrt{\frac{2}{3}}k(\alpha)\leq
0,$ (36)

where $k(\alpha)$ is defined by a nonlinear isotropic hardening law, as
proposed by Simo et al. [38, 44]:

$k(\alpha)=\sigma_{0}+\left(\sigma_{\infty}-\sigma_{0}\right)[1-\exp(-\delta\alpha)]+H\alpha,$
(37)

where $\sigma_{0}$ represents the initial flow stress, also called yield
stress, $\sigma_{\infty}$ the saturation flow stress, $\delta$ the saturation
exponent $\delta>0$, and $H$ the linear hardening coefficient. The yield
function $f$ defines the yield surface (when $f=0$), which classifies the
purely elastic response when $f<0$. When the yield condition is violated
($f>0$), the stress response to deformation consists of both plastic and
elastic components. Once the deformation enters the plastic regime, the
material cannot go back to its original shape, with a permanently plastic
deformation.

###  A.4 Return mapping algorithm

To integrate the material deformation over time, the return mapping algorithm,
which has been explained and widely applied in literature [53, 54, 55, 56,
47], is applied here. In this algorithm, the deviatoric part of the Kirchhoff
stress tensor
$\mathbf{s}^{pre}=\|\operatorname{dev}(\bm{\tau})\|=\mu\operatorname{dev}(\overline{\mathbf{b}}^{e})$
is obtained using the predicted update of $\mathbf{b}^{e}$, considering only
the elastic strain. The yield condition is then checked using Eq. (36) to
determine if plastic deformation occurs. If the yield condition is not
satisfied, the strain in the current step remains elastic, and the predicted
update $\mathbf{b}^{e}$ is considered acceptable. Otherwise, plastic
correction (returning map) is introduced to obtain the final $\mathbf{b}^{e}$
and $\mathbf{s}$ for the next step. The framework of this algorithm is
presented in Agorithm 1.

##  Appendix B Fluid-structure interaction model

In this appendix, referring to Zhao’s [24] algorithm, we briefly discuss the
porosity assumption and the corresponding relations, including porosity and
fluid saturation(B.1), and stress relations(B.2). In this simplified mixture
model, the solid and fluid densities $\rho_{s}$ and $\rho_{l}$, solid velocity
$\mathbf{v}^{s}$, and fluid saturation $\widetilde{a}$ are treated as state
variables, enabling the fluid velocity to be calculated referring to solid
velocity, rather than being an independent variable. This approach is
practically significant because it significantly reduces the complexity of the
system, as it eliminates the need for two sets of equations to describe the
fluid and solid separately.

###  B.1 Porosity and fluid saturation

Considering a representative volume element $dV$, the macroscopic porosity $a$
is defined as the ratio of the total volume of the pores $dV^{p}$ to $dV$,
yielding $a=\frac{dV^{p}}{dV}$. Note that $0<{a}<1$ holds for all cases.

When the porous solid is partially saturated by fluid, the fluid volume in the
representative element $dV$ is denoted by $dV^{l}$. The fluid saturation level
$\widetilde{a}$ can be defined as

$\widetilde{a}=\frac{dV^{l}}{dV}.$ (38)

Clearly, $\widetilde{a}$ is always less than or equal to the maximum possible
saturation ${a}$, i.e., $\widetilde{a}\leq{a}$. The locally effective fluid
density $\rho^{l}$, defined as the mass of the fluid per unit volume, varies
depending on the extent of fluid saturation and can be expressed as

$\rho^{l}=\frac{dm^{l}}{dV}=\frac{dm^{l}}{dV^{l}}\frac{dV^{l}}{dV}=\rho^{l}_{0}\widetilde{a},$
(39)

where $dm^{l}$ represents the mass of the fluid within a representative volume
element $dV$, ${\rho}^{l}_{0}$ the initial density of the fluid which is
assumed to be a constant for incompressible fluids.

###  B.2 Effective stress on solid

Following [57, 58, 59, 60], the total stress acting on the solid is the sum of
Cauchy stress $\bm{\sigma}^{s}$ and the pressure stress due to the presence of
the fluid phase $\bm{\sigma}^{l}$, written as:

$\bm{\sigma}=\bm{\sigma}^{s}+\bm{\sigma}^{l}=\bm{\sigma}^{s}-p^{l}\mathbf{I}.$
(40)

where $p^{l}$ is fluid pressure. For a hyper-elastic material, the
constitutive equation for the solid component is given by

$\bm{\sigma}^{s}=2\mu\mathbf{e}+\lambda\text{tr}(\mathbf{e})\mathbf{I},$ (41)

where the Eulerian-Almansi finite strain tensor $\mathbf{e}$ can be evaluated
by

$\mathbf{e}=\frac{1}{2}(\mathbf{I}-\mathbf{F}^{T}\mathbf{F}).$ (42)

The Lam$\acute{e}$ parameters $\lambda$ can be calculated via shear modulus
$\mu$ and bulk modulus K as $\lambda=K-\frac{2\mu}{3}$.

The excess fluid pressure simply depends on the fluid saturation level within
the porous solid element, with a function $p^{l}=p^{l}(\widetilde{a})$. The
relationship between fluid diffusion and the solid deformation satisfies a
fundamental principle: when fluid flows out of a given zone, the saturation
level decreases, resulting in a drop in pressure, and consequently, the
material tends to contact inwardly. Conversely, when fluid penetrates a porous
solid area, there exists a higher saturation level corresponding to a stronger
pressure force, leading to a material expansion. In the present model, this
behavior is described mathematically using a linear relation, taking the form

$p^{l}=C(\widetilde{a}-\widetilde{a}_{0}),$ (43)

where $C$ is a material constant, $\widetilde{a}_{0}$ the initial saturation.
Details can be referred to [60].

1Update deformation tensor  $\mathbf{F}_{n+1}=\Delta
t\frac{d\mathbf{F}}{dt},\quad\overline{\mathbf{F}}_{n+1}=[\operatorname{det}(\mathbf{F}_{n+1})]^{-1/3}\mathbf{F}_{n+1}$

2Predict the elastic
$\overline{\mathbf{b}}_{n+1}^{e,\mathrm{pre}}=\overline{\mathbf{F}}_{n+1}\overline{\mathbf{C}}_{p}^{n}\overline{\mathbf{F}}_{n+1}^{T},\quad\mathbf{s}_{n+1}^{\mathrm{pre}}=\mu\operatorname{dev}(\overline{\mathbf{b}}_{n+1}^{e,\mathrm{pre}})$

3Check the von Mises criterion
$f_{n+1}^{\mathrm{pre}}=s_{n+1}^{\mathrm{pre}}-\sqrt{\frac{2}{3}}k(\alpha_{n})$

4if  _$f_{n+1}^{\mathrm{pre}}\leq 0$_ then

5 Elastic state, Set $(\textbf{.})_{n+1}=(\textbf{.})_{n+1}^{\mathrm{pre}}$,
go to 1.

6 else

7 Perform 1 (the return mapping).

8 end if

9

10 Compute the normlized shear modulus
$\tilde{\mu}=\frac{1}{3}\operatorname{Tr}(\overline{\mathbf{b}}_{n+1}^{e,\mathrm{pre}})\mu$
Initiate  $\Delta\gamma=0$ Compute $\Delta\gamma$ so that
$\hat{f}(\Delta\gamma)=\left\|\mathbf{s}_{n+1}^{\text{trial
}}\right\|-\sqrt{\frac{2}{3}}k\left(\alpha_{n}+\sqrt{\frac{2}{3}}\Delta\gamma\right)-2\bar{\mu}\Delta\gamma=0$
$\hat{\mathbf{s}}=\mathbf{s}_{n+1}^{\mathrm{pre}}/s_{n+1}^{\mathrm{pre}}$
Return map
$\mathbf{s}_{n+1}=\mathbf{s}_{n+1}^{\mathrm{pre}}-2\tilde{\mu}\Delta\gamma\hat{\mathbf{s}},\quad\alpha_{n+1}=\alpha_{n}+\sqrt{\frac{2}{3}}\Delta\gamma$

11 Update stress
$J_{n+1}=det(\mathbf{F}_{n+1}),\quad\bm{\tau}_{n+1}=\frac{\kappa}{2}\left(J_{n+1}^{2}-1\right)\mathbf{I}+\mathbf{s}_{n+1},\quad\mathbf{P}_{n+1}=\bm{\tau}_{n+1}\mathbf{F}_{n+1}^{-T}$

Update local configuration
$\overline{\mathbf{b}}_{n+1}^{e}=\frac{1}{\mu}\mathbf{s}_{n+1}+\frac{1}{3}\operatorname{Tr}\left(\overline{\mathbf{b}}_{n+1}^{e,\text{
pre}}\right)\mathbf{I},\quad\overline{\mathbf{C}}_{p}^{n+1}=\overline{\mathbf{F}}_{n+1}^{-1}\overline{\mathbf{b}}_{n+1}^{e}\overline{\mathbf{F}}_{n+1}^{-T}$

Algorithm 1 Returning-mapping algorithm for $J_{2}$ flow theory with nonlinear
isotropic hardening.