# A Lattice Boltzmann Method for nonlinear solid mechanics in the reference
configuration

Erik Faust
Lehrstuhl für Technische Mechanik
Technische Universität Kaiserslautern
Postfach 3049, D-67653 Kaiserslautern
<EMAIL_ADDRESS>
& Alexander Schlüter
Lehrstuhl für Technische Mechanik
Technische Universität Kaiserslautern
Postfach 3049, D-67653 Kaiserslautern
<EMAIL_ADDRESS>
& Henning Müller
Institut für Mechanik
Technische Universität Darmstadt
Franziska-Braun-Straße 7, D-64287, Darmstadt
<EMAIL_ADDRESS>
& Ralf Müller
Institut für Mechanik
Technische Universität Darmstadt
Franziska-Braun-Straße 7, D-64287, Darmstadt
<EMAIL_ADDRESS>
Corresponding author.

(August 25, 2022)

###### Abstract

With a sufficiently fine discretisation, the Lattice Boltzmann Method (LBM)
mimics a second order Crank-Nicolson scheme for certain types of balance laws
Farag et al. (2021). This allows the explicit, highly parallelisable LBM to
efficiently solve the fundamental equations of solid mechanics: the
conservation of mass, the balance of linear momentum, and constitutive
relations.

To date, all LBM algorithms for solid simulation – see e.g. Murthy et al.
(2017), Escande et al. (2020), Schlüter et al. (2021) – have been limited to
the small strain case. Furthermore, the typical interpretation of the LBM in
the current (Eulerian) configuration is not easily extensible to large
strains, as large topological changes complicate the treatment of boundary
conditions.

In this publication, we propose a large deformation Lattice Boltzmann Method
for geometrically and constitutively nonlinear solid mechanics. To facilitate
versatile boundary modelling, the algorithm is defined in the reference
(Lagrangian) configuration.

_Keywords_ Lattice Boltzmann Method $\cdot$ Nonlinear Solid Mechanics $\cdot$
Large Deformation $\cdot$ Reference Configuration $\cdot$ Transient Simulation

##  1 Introduction

The Lattice Boltzmann Method (LBM) is a staple solver for simulations in fluid
dynamics. It is particularly attractive on account of its computational
efficiency and ease of discretisation (Krüger et al., 2017, p.55). Recently,
this efficiency and the promise of coupled fluid-structure simulations using
the LBM has also driven a rising research interest into Lattice Boltzmann (LB)
schemes for solid simulation Murthy et al. (2017); Escande et al. (2020).

Marconi and Chopard (2003); O’Brien et al. (2012); Schlüter et al. (2021), for
example, proposed LB algorithms of increasing generality based on wave
equation formulations for the Lamé-Navier equation. Marconi and Chopard
(2003); Schlüter et al. (2018); Müller et al. (2022) used similar schemes for
simulations in fracture mechanics. Finally, Murthy et al. (2017) and Escande
et al. (2020) developed a solid-mechanical LBM based on balance laws in the
small strain setting. These LB algorithms allow for efficient simulations of
transient phenomena in linear elastic continua. Additional developments in
boundary rules – see e.g. Schlüter et al. (2018, 2022); Faust et al. (2022) –
have also made it possible to model the behaviour of solid bodies with more
complicated geometries.

Until now, all such LB algorithms for solid simulation have been limited to
small deformations. Furthermore, the typical interpretation of the LBM in the
current (Eulerian) configuration is not easily extensible to large strains, as
large topological changes complicate the treatment of boundary conditions.

An extension to large deformations might open up new fields of application for
the LBM in solid mechanics: impact simulations and biomechanical tissue
modelling, for example, may eventually become feasible. Additionally,
nonlinear material models could allow for simulations of tires, seals, and
other polymer components under transient load. A formulation in the reference
configuration, meanwhile, could greatly simplify the modelling of boundaries
in certain situations with large deformations. Furthermore, it may also become
feasible to use the LBM to model both the fluid and the solid behaviour in
coupled fluid-structure simulations.

This work proposes a large displacement Lattice Boltzmann Method for
geometrically and materially nonlinear solid mechanics in the reference
configuration. The publication is organised as follows: section 2 briefly
motivates why the LBM can be used as a solver for balance laws in the
reference configuration, based on a recent numerical investigation of the LBM
by Farag et al. (2021). Section 3 then outlines an LB scheme from the
Lagrangian point of view, and section 4 presents an appropriate moment chain
for nonlinear solid mechanics. In section 5, an LB algorithm is proposed to
solve this moment chain, and section 6 covers the boundary rules used.
Finally, results obtained from the LBM are validated against FEM reference
solutions for two simple examples and one more challenging benchmark.

##  2 The LBM as a solver for moment chains

As demonstrated recently Farag et al. (2021), the Lattice Boltzmann Method
(LBM) mimics a second-order Crank-Nicolson scheme for balance laws of the form

$\frac{d}{dt}\Big{(}{}^{N}\Pi_{[\alpha]}\Big{)}+\frac{d}{dx_{\beta}}\Big{(}{}^{N+1}\Pi_{[\alpha]\beta}\Big{)}={}^{N}S_{[\alpha]}\,,$
(1)

where $\bm{{}^{N}\Pi}$ represents some $N$-th order tensor field and
$[\alpha]$ designates a set of indices $\alpha_{1},...,\alpha_{N}$, while
$\bm{{}^{N}S}$ is an $N$-th order source term. This is a profound realisation:
for sufficiently small time steps, the explicit, highly parallelisable LBM
acts as an implicit Runge-Kutta method in moment space. As it uses purely
numerical arguments, the derivation in Farag et al. (2021) further generalises
the LBM beyond the spatial particle transport context of the Boltzmann
equation (Krüger et al., 2017, p.21) in which it is usually viewed.
Consequently, we can interpret the temporal and spatial derivatives in (1) as
being defined in the reference (Lagrangian) configuration, and use the LBM to
solve conservation laws of the type

$d_{t}{}^{N}\Pi_{[\alpha]}+d_{\beta}{}^{N+1}\Pi_{[\alpha]\beta}={}^{N}S_{[\alpha]}\,,$
(2)

where $d_{t}$ and $d_{\alpha}$ are the material time derivative and material
gradient/divergence operators, respectively. When unravelled, the first three
‘chain links’ of this moment chain could also be written as

$\displaystyle d_{t}{}^{0}\Pi+d_{\alpha}{}^{1}\Pi_{\alpha}={}^{0}S\,,$
$\displaystyle
d_{t}{}^{1}\Pi_{\alpha}+d_{\beta}{}^{2}\Pi_{\alpha\beta}={}^{1}S_{\alpha}\,,$
$\displaystyle
d_{t}{}^{2}\Pi_{\alpha\beta}+d_{\gamma}{}^{3}\Pi_{\alpha\beta\gamma}={}^{2}S_{\alpha\beta}\,.$
(3)

Equation (2) describes balance laws with respect to material rather than
spatial coordinates, which is more convenient for solid simulations,
especially at large deformations. Solids do generally not assume the shape of
vessels containing them. Thus, boundary conditions must be formulated with
respect to parts of the material rather than regions in space. Furthermore,
the material time derivatives often appearing in solid-mechanical balance laws
yield convective terms in the Eulerian framework, but not in the Lagrangian
setting (Holzapfel, 2000, p.64-69). This fact simplifies derivations
considerably.

##  3 The LBM (in the reference configuration)

This section briefly discusses the central steps of the LBM algorithm, which
is formulated in the reference configuration here. The interested reader is
referred to Krüger et al. (2017) for an excellent textbook covering the LBM in
detail.

We introduce a lattice consisting of lattice sites $\bm{X}\in
L\subset\Omega_{0}$ in the reference configuration. The sites are spaced
regularly at distances of $\Delta X$ in each coordinate direction. For the
sake of brevity, only two-dimensional D2Q9 lattices are considered in this
work. In such a lattice, each lattice site is connected to its $8$ nearest
(direct and diagonal) neighbours via $8$ lattice links. Furthermore, a set of
$9$ lattice velocities111Note that this includes one zero velocity
$\bm{C_{i}}=\bm{0}$. $\bm{C_{i}}$ are defined to cover the distance between
one lattice site $\bm{X}$ and its neighbours $\bm{X}+\bm{C_{i}}\Delta t$ in
one time step $\Delta t$ (Krüger et al., 2017, p.86-88).

The LBM performs explicit operations on distribution functions $\bar{f}_{i}$
(one per lattice velocity $\bm{C_{i}}$) on this lattice: in each iteration,
these are first relaxed towards the value of the equilibrium distribution
function $f_{i}^{eq}$ with relaxation time $\bar{\tau}$,

$\bar{f}_{i}^{col}=\bar{f}_{i}-\frac{\Delta
t}{\bar{\tau}}(\bar{f}_{i}-f_{i}^{eq})+\Delta t\Big{(}1-\frac{\Delta
t}{2\bar{\tau}}\Big{)}\psi_{i}\,,$ (4)

and a source term contribution $\psi_{i}$ is added. The choice of equilibrium
distribution function
$f_{i}^{\text{eq}}({}^{0}\Pi,\bm{{}^{1}\Pi},\bm{{}^{2}\Pi},...)$ and source
term $\psi_{i}$ determines the equations modelled in moment space, as
discussed below. Here and in the following, we assume a simple BKGW collision
term (see (Krüger et al., 2017, p.98-101)), and a second-order He forcing
scheme He et al. (1998) (see (Krüger et al., 2017, p.233-239)).

The results of this local operation are then propagated to neighbouring
lattice sites in the streaming step

$\bar{f}_{i}(\bm{X}+\bm{C}_{i}\Delta t,t+\Delta t)=\bar{f}_{i}^{col}\,.$ (5)

The link between the distribution functions $\bar{f}_{i}$ and the macroscopic
moments $\bm{{}^{N}\Pi}$, meanwhile, is given via (Krüger et al., 2017, p.239)

${}^{N}\Pi_{[\alpha]}=\sum_{i}C_{i[\alpha]}^{\otimes^{N}}\bar{f}_{i}+\frac{1}{2}{}^{N}S_{[\alpha]}\Delta
t\,.$ (6)

Here, the dyadic product of the velocity vector $C_{i}$ with itself, taken $N$
times, is denoted via a dyadic power operator $\otimes^{N}$ in the
superscript222This dyadic power operator is defined such that
$\bm{C}_{i}^{\otimes^{1}}=\bm{C}_{i}$,
$\bm{C}_{i}^{\otimes^{2}}=\bm{C}_{i}\otimes\bm{C}_{i}$,
$\bm{C}_{i}^{\otimes^{2}}=\bm{C}_{i}\otimes\bm{C}_{i}\otimes\bm{C}_{i}$,
etc... For example, the first three moments are computed via

$\displaystyle{}^{0}\Pi=\sum_{i}\bar{f}_{i}+\frac{1}{2}{}^{0}S\Delta t\,,$
$\displaystyle{}^{1}\Pi_{\alpha}=\sum_{i}C_{i\alpha}\bar{f}_{i}+\frac{1}{2}{}^{1}S_{\alpha}\Delta
t\,,$
$\displaystyle{}^{2}\Pi_{\alpha\beta}=\sum_{i}C_{i\alpha}C_{i\beta}\bar{f}_{i}+\frac{1}{2}{}^{2}S_{\alpha\beta}\Delta
t\,.$

The moment computation by (6) is performed after the collision and streaming
steps in each iteration. The moments are then used to compute the equilibrium
distribution function
$f_{i}^{eq}({}^{0}\Pi,\bm{{}^{1}\Pi},\bm{{}^{2}\Pi},...)$ for use in the
collision step (4) in the next iteration (Krüger et al., 2017, p.67).
Similarly, the source term $\psi_{i}$ in (4) is computed as
$\psi_{i}({}^{0}S,\bm{{}^{1}S},\bm{{}^{2}S},...)$ from the source terms
$\bm{{}^{N}S}$ found on the right hand sides of the moment chain in (2) Farag
et al. (2021). Which conservation laws can be modelled depends on the
particular choice of equilibrium distribution function $f_{i}^{eq}$ and source
term $\psi_{i}$ Farag et al. (2021). Generally speaking,
$f_{i}^{eq}({}^{0}\Pi,\bm{{}^{1}\Pi},\bm{{}^{2}\Pi},...)$ and
$\psi_{i}({}^{0}S,\bm{{}^{1}S},\bm{{}^{2}S},...)$ are functions in the moments
$\bm{{}^{N}\Pi}$ and source terms $\bm{{}^{N}S}$, respectively.

Finally, boundary rules must be defined to determine the value of the
distribution functions $\bar{f}_{i}$ in accordance with the macroscopic
boundary conditions on the moments $\bm{{}^{N}\Pi}$. In a recent preprint
Faust et al. (2022), we adapted the popular bounce-back (see (Krüger et al.,
2017, p.175-181)) and anti-bounce-back boundary rules (see (Krüger et al.,
2017, p.200)) to model the behaviour of Dirichlet and Neumann boundaries with
an LBM for linear elasticity by Murthy et al. (2017) and Escande et al.
(2020). In this work, we use the same methodology for Dirichlet boundaries,
and adjust the modified anti-bounce-back boundary rule to handle boundary
conditions in the first Piola-Kirchhoff stress in the reference configuration
(see section 6).

##  4 A moment chain for nonlinear solid mechanics

To facilitate their solution via the LBM, we would like to rewrite the
fundamental equations of solid mechanics – the conservation of mass, the
balance of linear momentum, and the material law – in the moment chain form
given in (2).

From the Lagrangian point of view, mass conservation becomes trivial: the
reference density $\rho_{0}$ remains constant, and can be treated as a
material parameter, i.e. $d_{t}\rho_{0}=0$ (Holzapfel, 2000, p.135). Thus, we
do not need to account for mass conservation with the zeroth-order equation in
(2), and remain flexible as to what to model with this equation and with the
(scalar) zeroth moment ${}^{0}\Pi$.

The balance of linear momentum

$\rho_{0}a_{\alpha}=d_{\beta}P_{\alpha\beta}+\rho_{0}b_{\alpha}\,,$ (7)

meanwhile, relates the acceleration $\bm{a}$ to the divergence of the first
Piola-Kirchhoff stress tensor $\bm{P}$ as well as the body force per unit mass
$\bm{b}$ (Holzapfel, 2000, p.143). The acceleration $\bm{a}$ is defined as the
material time derivative of the velocity $\bm{v}$, meaning that the inertial
term in (7) can be rewritten as

$\rho_{0}a_{\alpha}=\rho_{0}d_{t}v_{\alpha}=d_{t}(\rho_{0}v_{\alpha})=d_{t}\bar{j}_{\alpha}\,,$
(8)

with $\bm{\bar{j}}=\rho_{0}\bm{v}$ denoting the linear momentum density in the
Lagrangian configuration. The balance of linear momentum thus becomes

$d_{t}\bar{j}_{\alpha}-d_{\beta}P_{\alpha\beta}=\rho_{0}b_{\alpha}\,,$ (9)

which is more convenient for implementation in an LB algorithm: the time
derivative of the linear momentum density $\bm{\bar{j}}$ appearing here fits
neatly into the first-order equation in (2), meaning that $\bm{\bar{j}}$ might
be used as a first moment $\bm{{}^{1}\Pi}$.

Finally, a material law is required to account for the evolution of the first
Piola-Kirchhoff stress tensor $\bm{P}$. In this publication, we assume
hyperelastic behaviour of the form333Note that material law formulations in
terms of alternative strain measures are also possible. The deformation
gradient is chosen here because it can easily be calculated from the
displacement $\bm{u}$. $\bm{P}(\bm{H})$, with the displacement gradient
$\bm{H}$

$H_{\alpha\beta}=d_{\beta}u_{\alpha}\,.$ (10)

Here, $\bm{u}$ denotes the displacement field. Accounting for such a material
law in an LB algorithm is not trivial: ideally, we would like to use the
second-order equation in (2), i.e.

$d_{t}{}^{2}\Pi_{\alpha\beta}+d_{\gamma}{}^{3}\Pi_{\alpha\beta\gamma}=0\,,$
(11)

to model the evolution of the first Piola-Kirchhoff stress tensor $\bm{P}$ as
a second moment $\bm{{}^{2}\Pi}$ via the divergence of a third-order tensor
$\bm{Q}$ to be used as a third moment $\bm{{}^{3}\Pi}$, i.e.

$-d_{t}P_{\alpha\beta}+d_{\gamma}Q_{\alpha\beta\gamma}=0\,.$ (12)

However, an analytical solution for a third order tensor $\bm{Q}$ satisfying
this relation for an arbitrary hyperelastic material law $\bm{P}(\bm{H})$ can
– to the best of our knowledge – not be identified easily444This equation
might also be solved numerically. If discretised with finite difference
approximations to the derivatives, (12) yields a very sparse linear system. In
each iteration of an LB algorithm, the value of the unknown third order tensor
$\bm{Q}$ from the previous iteration might be used as an initial value for an
iterative solver. Due to the small time step, the change in the values of
$\bm{P}$ and thus $\bm{Q}$ from one iteration to the next is limited,
promising quick convergence. Note that (12) can be treated separately for
coordinate indices $\alpha$ and $\beta$..

We instead use the second-order equation in (2) to model a part of the
material law $\bm{\bar{P}}(\bm{H})$ which is linear in the displacement
gradient $\bm{H}$, i.e.

$d_{t}\bar{P}_{\alpha\beta}+d_{\gamma}\bar{Q}_{\alpha\beta\gamma}=0\,.$ (13)

It turns out that for the linear expression

$\bar{Q}_{\alpha\beta\gamma}=C_{s}^{2}(\bar{j}_{\alpha}\delta_{\beta\gamma}+\bar{j}_{\beta}\delta_{\alpha\gamma}+\bar{j}_{\gamma}\delta_{\alpha\beta})\,,$
(14)

which can readily be replicated by the LBM Murthy et al. (2017),
$\bm{\bar{P}}$ turns out to be

$\bar{P}_{\alpha\beta}=-\mu(d_{\beta}u_{\alpha}+d_{\alpha}u_{\beta}+d_{\gamma}u_{\gamma}\delta_{\alpha\beta})\,.$
(15)

Here, $C_{s}=\sqrt{\nicefrac{{\mu}}{{\rho_{0}}}}$ is the shear wave speed. In
a previous publication Faust et al. (2022), we referred to $\bm{\bar{P}}$ as
the Poisson stress tensor due to its similarity to the Cauchy stress tensor
for Poisson solids555For small strains and Poisson solids with $\lambda=\mu$
($\nu=\nicefrac{{1}}{{4}}$), $\bm{\bar{P}}=-\bm{\sigma}$, hence the name.. For
the sake of convenience, the same name will be used in the following.

The (nonlinear) remainder of the material law can then be modelled using an
additional artificial source term in the balance of linear momentum, i.e.

$d_{t}\bar{j}_{\alpha}+d_{\beta}\bar{P}_{\alpha\beta}=\rho_{0}b_{\alpha}+d_{\beta}\Big{(}P_{\alpha\beta}(\bm{H})+\bar{P}_{\alpha\beta}(\bm{H})\Big{)}\,,$
(16)

with the modified source term
$\bm{{}^{1}S}=\rho_{0}\bm{b}+\text{div}(\bm{P}(\bm{H})+\bm{\bar{P}}(\bm{H}))$.

As mentioned previously, we remain flexible as to what to model with the first
equation in (2). We are committed to $\bm{{}^{1}\Pi}=\bm{\bar{j}}$ by the
considerations above. For simplicity, we could set ${}^{0}S=0$ and obtain

$d_{t}r+d_{\alpha}\bar{j}_{\alpha}=0\,,$ (17)

i.e.,

${}^{0}\Pi=r=-\rho_{0}d_{\alpha}u_{\alpha}\,.$ (18)

With all the above, we have identified a moment chain formulation for
nonlinear solid mechanics in the reference configuration, which is given by

$\displaystyle d_{t}r+d_{\alpha}\bar{j}_{\alpha}=0\,,$ $\displaystyle
d_{t}\bar{j}_{\alpha}+d_{\beta}\bar{P}_{\alpha\beta}=\rho_{0}b_{\alpha}+d_{\beta}\Big{(}P_{\alpha\beta}(\bm{H})+\bar{P}_{\alpha\beta}(\bm{H})\Big{)}\,,$
$\displaystyle
d_{t}\bar{P}_{\alpha\beta}+d_{\gamma}\bar{Q}_{\alpha\beta\gamma}=0\,.$ (19)

The zeroth moment $r=-\rho_{0}d_{\alpha}u_{\alpha}$ is proportional to the
trace of the deformation gradient $\bm{H}$ and thus encodes first-order
information about tensile and compressive strains. Its evolution is determined
by the linear momentum density $\bm{\bar{j}}$, which is used as a first
moment. The evolution of $\bm{\bar{j}}$ is determined by the divergence of the
second moment $\bm{\bar{P}}$ and by the first order source terms on the right
hand side of the second equation in (4).

The second moment $\bm{\bar{P}}$ encodes a part of the material law which is
linear in $\bm{H}$. The nonlinear remainder of the material law is modelled
via an additional source term – note the opposing signs of $\bm{P}$ and
$\bm{\bar{P}}$, which mean that the linear terms accounted for by
$\bm{\bar{P}}$ are ’subtracted from’ the first Piola-Kirchhoff tensor for the
source term.

As shown by Farag et al. (2021), the LBM mimics a second-order Crank-Nicolson
scheme for moment chains like this666For the LBM to mimic a fully implicit
Crank-Nicolson solver, the source terms must be evaluated – for some iteration
starting at time $t$ – at time $t$ for the collision operation in (4) and at
time $t+\Delta t$ for the moment computation via (6).. The relaxation and
streaming operations in an LBM algorithm account for the left hand sides of
(4), while the source terms – both the bulk force and the nonlinear parts of
the material law – are modelled via the forcing contribution in the collision
step.

##  5 An LBM for nonlinear solid mechanics

In the preceding section, we introduced a moment chain for nonlinear, large-
deformation solid mechanics in the reference configuration. As moments, we
selected ${}^{0}\Pi=r$, $\bm{{}^{1}\Pi}=\bm{\bar{j}}$,
$\bm{{}^{2}\Pi}=\bm{\bar{P}}$, and $\bm{{}^{3}\Pi}=\bm{\bar{Q}}$, with

$\displaystyle\bar{P}_{\alpha\beta}=-\mu(d_{\beta}u_{\alpha}+d_{\alpha}u_{\beta}+d_{\gamma}u_{\gamma}\delta_{\alpha\beta})\,,$
$\displaystyle\bar{Q}_{\alpha\beta\gamma}=C_{s}^{2}(\bar{j}_{\alpha}\delta_{\beta\gamma}+\bar{j}_{\beta}\delta_{\alpha\gamma}+\bar{j}_{\gamma}\delta_{\alpha\beta})\,,$

from (15) and (14). These moments may readily be replicated via the LBM: the
equilibrium distribution function Murthy et al. (2017)

$f_{i}^{eq}=w_{i}\Big{(}r+\frac{1}{c_{s}^{2}}C_{i\alpha}\bar{j}_{\alpha}+\frac{1}{2C_{s}^{4}}(\bar{P}_{\alpha\beta}-rC_{s}^{2}\delta_{\alpha\beta})(C_{i\alpha}C_{i\beta}-C_{s}^{2}\delta_{\alpha\beta})\Big{)}\,,$
(20)

yields precisely the tensors $r$, $\bm{\bar{j}}$, $\bm{\bar{P}}$, and
$\bm{\bar{Q}}$ via the moments,

$\displaystyle r$ $\displaystyle=\sum_{i}f_{i}^{eq}\,,$
$\displaystyle\bar{j}_{\alpha}$
$\displaystyle=\sum_{i}f_{i}^{eq}C_{i\alpha}\,,$
$\displaystyle\bar{P}_{\alpha\beta}$
$\displaystyle=\sum_{i}f_{i}^{eq}C_{i\alpha}C_{i\beta}\,,$
$\displaystyle\bar{Q}_{\alpha\beta\gamma}$
$\displaystyle=\sum_{i}f_{i}^{eq}C_{i\alpha}C_{i\beta}C_{i\gamma}\,,$ (21)

which guarantees that the moments of $\bar{f}_{i}$ behave as in (6) Farag et
al. (2021). The source term Murthy et al. (2017),

$\psi_{i}=w_{i}\frac{1}{c_{s}^{2}}C_{i\alpha}{}^{1}S_{\alpha}\,,$ (22)

with

${}^{1}S_{\alpha}=\rho_{0}b_{\alpha}+d_{\beta}(P_{\alpha\beta}(\bm{H})+\bar{P}_{\alpha\beta}(\bm{H}))\,,$
(23)

accounts for the remainder of the material law. Since the displacement
$\bm{u}$ is not directly available in the LBM algorithm, it must be computed
from the linear momentum density via integration. Here, a trapezoidal rule

$u_{\alpha}(x,t+\Delta
t)=u_{\alpha}(x,t)+\frac{1}{2\rho_{0}}\Big{(}\bar{j}_{\alpha}(x,t+\Delta
t)+\bar{j}_{\alpha}(x,t)\Big{)}\,,$ (24)

is used. The displacement gradient $H_{\alpha\beta}=d_{\beta}u_{\alpha}$ can
then be determined from the displacement and be used to compute the first
Piola-Kirchhoff and Poisson stress tensors $\bm{P}$ and $\bm{\bar{P}}$ for use
in the source term $\bm{{}^{1}S}$. We utilise central finite difference
stencils for the gradient and divergence computations.

It is important to note that the source term $\bm{{}^{1}S}$ is required at two
points in each iteration: firstly, to compute the forcing contribution to the
collision step (4), and secondly, in the moment computation (6). Note that,
for the LB algorithm to mimic a fully implicit Crank-Nicolson solver,
$\bm{{}^{1}S}$ must be evaluated at time $t$ in the collision step, but at
time $t+\Delta t$ in the moment computation, for some iteration starting at
time $t$. Since the displacement at time $t+\Delta t$ is unknown prior to the
moment computation, we instead evaluate $\bm{{}^{1}S}$ at time $t$ in (6).
Thus, we do not mimic a fully explicit second-order Crank-Nicolson scheme for
the nonlinear component of the material law (but do for the part contained in
$\bm{\bar{P}}$). The approximation error resulting from this operation is of
order $\mathcal{O}(\Delta t^{2})$. If we approximate $\bm{\bar{j}}$ at
$t+\Delta t$ by linear extrapolation from $t$, and subsequently use the
trapezoidal rule to obtain $\bm{u}$ at $t+\Delta t$, the error reduces to
$\mathcal{O}(\Delta t^{3})$. We therefore obtain a partially explicit scheme
for the part of the material law not contained in $\bm{\bar{P}}$, albeit with
a rather small error777If this resulted in issues with the stability of the
method, we could implement a predictor-corrector scheme for the computation of
$\bm{u}$ at $t+\Delta t$. In the numerical experiments performed by the
authors, however, this was not necessary..

Note that this is a ’quick and dirty’ strategy: it would be preferable, from
the viewpoint of efficiency and stability, to do away with this source term,
but we have not been able to do so yet. Encouragingly, the LB scheme presented
here seems to perform as desired in spite of this. We expect further
performance gains from a ’cleaner’ algorithm.

The complete LB algorithm – used to solve the solid-mechanical moment chain in
(2) – is outlined in algorithm 1. A brief discussion of the boundary handling
scheme follows in the next section.

Set _$\Delta t$ , $\Delta X$, $C_{s}$, $\bm{C}_{i}$, $w_{i}$, $\bar{\tau}$,
$\lambda$, $\mu$, $\rho_{0}$_;

Initialise _Lattice, Boundary, Moments: $r$, $\bm{\bar{j}}$, $\bm{\bar{P}}$,
$\bm{\bar{Q}}$_;

// initialise distribution functions

$\bar{f}_{i}\leftarrow
w_{i}\Big{(}r+\frac{1}{c_{s}^{2}}C_{i\alpha}\bar{j}_{\alpha}+\frac{1}{2C_{s}^{4}}(\bar{P}_{\alpha\beta}-rc_{s}^{2}\delta_{\alpha\beta})(C_{i\alpha}C_{i\beta}-C_{s}^{2}\delta_{\alpha\beta})\Big{)}$;

$t\leftarrow 0$;

// main simulation loop

while  _$t\leq t_{max}$_ do

for  _$X\in L$_ do

${}^{1}S_{\alpha}\leftarrow\rho_{0}b_{\alpha}+d_{\beta}(\bar{P}_{\alpha\beta}+P_{\alpha\beta})$;

for  _$i\leftarrow 0$ to $9$_ do

// calculate equilibrium distribution function

$f_{i}^{eq}\leftarrow
w_{i}\Big{(}r+\frac{1}{C_{s}^{2}}C_{i\alpha}\bar{j}_{\alpha}+\frac{1}{2C_{s}^{4}}(\bar{P}_{\alpha\beta}-rC_{s}^{2}\delta_{\alpha\beta})(C_{i\alpha}C_{i\beta}-C_{s}^{2}\delta_{\alpha\beta})\Big{)}$;

// compute source

$\psi_{i}\leftarrow w_{i}\frac{1}{C_{s}^{2}}C_{i\alpha}{}^{1}S_{\alpha}$;

// collide

$\bar{f}_{i}^{col}\leftarrow\bar{f}_{i}-\frac{\Delta
t}{\bar{\tau}}\big{(}\bar{f}_{i}-\bar{f}_{i}^{eq}\big{)}+\Big{(}1-\frac{\Delta
t}{2\bar{\tau}}\Big{)}\psi_{i}\Delta t$;

// stream (except for $X+C_{i}\Delta t$ crossing boundary)

$\bar{f}_{i}(X+C_{i}\Delta t,t+\Delta t)\leftarrow\bar{f}_{i}^{col}(X,t)$;

end for

end for

// calculate missing $\bar{f}_{i}$ at the boundary

for  _$X\in\partial L_{u}$_ do

for  _$i$ crossing $\partial\Omega_{u}$_ do

$\bar{f}_{\bar{i}}(X,t+\Delta
t)\leftarrow\bar{f}_{i}^{col}(X,t)-\frac{2}{C_{s}^{2}}w_{i}C_{i\alpha}\bar{j}_{\alpha}^{bd}$;

end for

end for

for  _$X\in\partial L_{t}$_ do

for  _$i$ crossing $\partial\Omega_{t}$_ do

$\bar{f}_{\bar{i}}(X,t+\Delta
t)\leftarrow-\bar{f}_{i}^{col}(X,t)+2w_{i}\Big{(}r^{bd}+\frac{1}{2C_{s}^{4}}(\bar{P}_{\alpha\beta}^{bd}-rC_{s}^{2}\delta_{\alpha\beta})(C_{i\alpha}C_{i\beta}-C_{s}^{2}\delta_{\alpha\beta})\Big{)}$;

end for

end for

// calculate moments (explicit here)

for  _$X\in L$_ do

${}^{1}S_{\alpha}\leftarrow\rho_{0}b_{\alpha}+d_{\beta}(\bar{P}_{\alpha\beta}+P_{\alpha\beta})$;

$r\leftarrow\sum_{i}\bar{f}_{i}$,
$\bar{j}_{\alpha}\leftarrow\sum_{i}C_{i\alpha}\bar{f}_{i}+\frac{\Delta
t}{2}{}^{1}S_{\alpha}$,

end for

// update displacement and stress

for  _$X\in L$_ do

$u_{\alpha}(X,t+\Delta t)\leftarrow
u_{\alpha}(X,t)+\nicefrac{{1}}{{2\rho_{0}}}\Big{(}\bar{j}_{\alpha}(X,t+\Delta
t)+\bar{j}_{\alpha}(X,t)\Big{)}$;

$H_{\alpha\beta}\leftarrow d_{\beta}u_{\alpha}$, $P_{\alpha\beta}\leftarrow
P_{\alpha\beta}(\bm{H})$,
$\bar{P}_{\alpha\beta}\leftarrow\bar{P}_{\alpha\beta}(\bm{H})$;

end for

$t\leftarrow t+\Delta t$;

end while

Algorithm 1 An LB algorithm for nonlinear solid mechanics.

##  6 Boundary conditions

In a recent publication Faust et al. (2022), we adapted the the popular
bounce-back (see (Krüger et al., 2017, p.175-181)) and anti-bounce-back
boundary rules (see (Krüger et al., 2017, p.200)) to approximate the behaviour
of Dirichlet and Neumann boundaries with an LBM for linear elasticity by
Murthy et al. (2017) and Escande et al. (2020). Here, we use an analogous
methodology to impose boundary values in the linear momentum density $\bm{j}$

$\bar{f}_{\bar{i}}(\bm{x},t+\Delta
t)=\bar{f}_{i}^{\text{col}}(\bm{x},t)-\frac{2}{C_{s}^{2}}w_{i}C_{i\alpha}\bar{j}_{\alpha}^{*}\,,$
(25)

and on the Poisson stress tensor $\bm{\bar{P}}$

$\bar{f}_{\bar{i}}(\bm{x},t+\Delta
t)=-\bar{f}_{i}^{\text{col}}(\bm{x},t)+2w_{i}\Big{(}r^{\text{bd}}+\frac{1}{2C_{s}^{4}}(\bar{P}_{\alpha\beta}^{*}-r^{\text{bd}}C_{s}^{2}\delta_{\alpha\beta})(C_{i\alpha}C_{i\beta}-C_{s}^{2}\delta_{\alpha\beta})\Big{)}\,.$
(26)

Boundary values in the Poisson tensor may be obtained from the boundary
stresses $\bm{T^{*}}$ via

$\displaystyle P_{\alpha\beta}N_{\beta}=T_{\alpha}^{*}\,,$
$\displaystyle-\bar{P}_{\alpha\beta}N_{\beta}+(P_{\alpha\beta}+\bar{P}_{\alpha\beta})N_{\beta}=T_{\alpha}^{*}\,,$
$\displaystyle\bar{P}_{\alpha\beta}N_{\beta}=-T_{\alpha}^{*}+(P_{\alpha\beta}+\bar{P}_{\alpha\beta})N_{\beta}\,.$

As in Faust et al. (2022), we transform these equations into a coordinate
system normal to the boundary, in which the quasi-boundary tractions
$-\bm{T^{*}}+(\bm{P}+\bm{\bar{P}})\bm{N}$ may simply be substituted into the
first row and column of $\bm{\bar{P}}$, while the remaining entries are
extrapolated to the boundary. The Poisson and first Piola-Kirchhoff tensors
appearing in the second term are also extrapolated to the boundary.

##  7 Numerical experiments

Having proposed an LB algorithm for nonlinear solid mechanics in the previous
sections, we now consider example problems for validation. In the absence of
analytical solutions for nonlinear dynamic field problems, we instead resort
to comparisons with Finite Element (FE) solutions. Firstly, two simple tension
and shear load cases are investigated, and a discussion of a more challenging
example follows. For the sake of simplicity, only the two-dimensional plane
strain case is treated in this publication.

The LBM is used to model the response of a compressible neo-Hooke solid to
these load cases. For the two-dimensional case, the strain energy density
function could be written as

$W=\frac{\mu}{2}(I_{C}-2)+\frac{\lambda}{4}(J^{2}-1)-\frac{\lambda}{2}\log
J-\mu\log J\,,$ (27)

and Lamé parameters $\lambda=\mu=1$ as well as density $\rho_{0}=1$ are
assumed throughout. Here, $\bm{C}$ is the right Cauchy-Green tensor
$\bm{C}=(\bm{H}+\bm{I})^{T}(\bm{H}+\bm{I})$ – of which $I_{C}$ is the first
invariant – and $J$ the volume ratio $J=\text{det}(\bm{H}+\bm{I})$. For non-
dimensionalisation, we use $l_{\text{ref}}=u_{\text{ref}}=l$.

###  7.1 Simple tension and simple shear

As a simple benchmark geometry, we consider a square block with side length
$l$. This block – pictured in figure 2 – is first subjected to a pure tensile
load $T^{*}$ in the $x_{2}$ direction. $T^{*}$ is increased linearly from
$T^{*}=0\nicefrac{{\mu l}}{{l}}$ at $t=0l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$
to $T^{*}=1\nicefrac{{\mu l}}{{l}}$ at
$t=1l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$, and then held at this level until
the end of the simulation at $t=3l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$.

In the second example, a pure shear load $T^{*}$ is applied to the same block,
as displayed in figure 4. $T^{*}$ is increased linearly from
$T^{*}=0\nicefrac{{\mu l}}{{l}}$ at $t=0l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$
to $T^{*}=0.05\nicefrac{{\mu l}}{{l}}$ at
$t=1l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$, and then held at this level until
the end of the simulation888A reduced load is used for the shear load case due
to the reduced stiffness of the block in the tangential direction. Due to the
lower speed of shear waves, the simulation is run for a longer time. at
$t=8l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$.

A lattice spacing of $\Delta X=0.025l$ is used in both cases, yielding a time
step of $\Delta t\approx 1.443\times
10^{-2}l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$ for the LBM. A FE simulation
using the same material parameters is performed using FEAP Taylor (2014). The
displacements at points $P_{1}=(\nicefrac{{\Delta
x}}{{2}},0.5-\nicefrac{{\Delta x}}{{2}})l$, $P_{2}=(-0.5+\nicefrac{{\Delta
x}}{{2}},0.5-\nicefrac{{\Delta x}}{{2}})l$, and $P_{3}=(-0.5+\nicefrac{{\Delta
x}}{{2}},\nicefrac{{\Delta x}}{{2}})l$ are then postprocessed from the
simulation results produced by the LBM and the FEM.

Figure 1: Square plate subjected to an in-plane tension load in the
$x_{2}$-direction.

Figure 2: Displacements postprocessed from LBM (blue) and FE (red) simulations
of a neo-Hooke solid under tension loading.

For the simple tension simulations, figure 2 shows the $x_{2}$ component of
the displacement $u_{2}^{(1)}$ at $P_{1}$, the $x_{2}$ component of the
displacement $u_{2}^{(2)}$ at $P_{2}$, and the $x_{1}$ component of the
displacement $u_{1}^{(3)}$ at $P_{3}$. While $u_{2}^{(1)}$ and $u_{2}^{(2)}$
are measures of the tensile response of the block in the loaded direction,
$u_{1}^{(3)}$ quantifies the lateral contraction. As is apparent from the
figure, the LBM and FEM are in agreement with regard to the computed
displacements throughout the simulation, with the LBM yielding slightly lower
displacements in general. Furthermore, the LBM predicts slightly later peaks
in the displacement curves. While the match between LBM and FEM is not
perfect, the results suggest that the LB scheme proposed above captures the
large-strain tensile response of the square block in figure 2 roughly as well
as the FEM.

Figure 3: Square plate subjected to an in-plane shear load in the
$x_{2}$-direction.

Figure 4: Displacements postprocessed from LBM (blue) and FE (red) simulations
of a neo-Hooke solid under shear loading.

Similarly, figure 4 shows the $x_{1}$ and $x_{2}$ components of the
displacement $u_{1}^{(2)}$ and $u_{2}^{(2)}$ at $P_{2}$, as well as the
$x_{1}$ component of the displacement $u_{1}^{(3)}$ at $P_{3}$. $u_{1}^{(2)}$
and $u_{2}^{(2)}$ measure the translation of the top left corner of the block
in figure 4 in the horizontal and vertical directions under shear load,
respectively. $u_{1}^{(3)}$, meanwhile, quantifies the horizontal translation
of the point $P_{3}$ in the centre of the left side of the block. Again, the
agreement between the LBM and the FEM is encouraging, though the rebound
predicted by the LB scheme is less extreme. Once more, solutions computed via
the LBM also peak slightly later. As both the FEM and the LBM produce only
approximations to the analytical solution, it is difficult to make more
precise judgements on the correctness of either result. In general, the LBM
and FEM produce solutions which agree decently qualitatively and
quantitatively, suggesting that the LB scheme captures a majority of the large
strain shear mode response of the block in figure 4.

###  7.2 Plate with hole under sudden loading and unloading

As a slightly more challenging example, we consider a quadratic plate with
side length $l$ featuring a quadratically shaped hole with side length
$e=0.4l$. At $t=0l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$, this plate is suddenly
loaded with a traction of $T^{*}=0.1\nicefrac{{\mu l}}{{l}}$ in the
$x_{2}$-direction. The traction is held constant until
$t=1l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$, at which point the load is removed
suddenly and $T^{*}=0\nicefrac{{\mu l}}{{l}}$ until the end of the simulation
at $t=3l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$. Once more, a neo-Hookean
material model is used. A lattice spacing of $\Delta X=0.0125l$ is chosen,
yielding a time step of $\Delta t\approx 7.217\times
10^{-3}l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$.

The displacements computed by the LBM and the FEM for
$Q_{1}=(0.2+\nicefrac{{\Delta X}}{{2}},\nicefrac{{\Delta X}}{{2}})l$ and
$Q_{2}=(\nicefrac{{\Delta X}}{{2}},0.5-\nicefrac{{\Delta X}}{{2}})l$, in the
$x_{1}$ and $x_{2}$ directions, respectively, are shown in figure 6. As is
clear from the figure, the normal and transverse displacements computed by the
FEM and LBM in response to the sudden loading and unloading are in excellent
agreement throughout the simulation. The plate is deformed considerably in
response to the load, being lengthened by approximately $18\%$ in the
$x_{2}$-direction.

Figure 5: Plate with hole subjected to an in-plane tension load in the
$x_{2}$-direction.

Figure 6: Displacements postprocessed from LBM (blue) and FE (red) simulations
of a neo-Hooke solid under tension loading.

Furthermore, the shear component $\sigma_{12}$ of the Cauchy stress
$\bm{\sigma}$ is computed at $Q_{3}=(0.3+\nicefrac{{\Delta
X}}{{2}},0.2+\nicefrac{{\Delta X}}{{2}})l$ and $Q_{4}=(0.2+\nicefrac{{\Delta
X}}{{2}},0.3+\nicefrac{{\Delta X}}{{2}})l$. Additionally, the
$x_{2}$-direction tensile stress component $\sigma_{22}$ is evaluated at
$Q_{1}$, $Q_{3}$, and $Q_{5}=(0.5-\nicefrac{{\Delta X}}{{2}},\nicefrac{{\Delta
X}}{{2}})l$. There is nothing fundamental about this choice of points, they
are simply chosen because they yield high peak stress values and an
interesting evolution of the stress in time.

Figure 7: Shear stresses postprocessed from LBM (blue) and FE (red)
simulations of a neo-Hooke solid under tension loading.

The evolution of the shear stresses computed by the LBM and the FEM for the
selected points is shown in figure 7. The two numerical solution match well,
yielding very similar overshoot profiles. The fluctuations induced by passing
S-waves are also very similar in both simulations. Generally, the LBM predicts
slightly lower shear stresses throughout.

Meanwhile, figure 8 outlines the evolution of the tensile stress
$\sigma_{22}$, the magnitude of which is significantly higher generally. Here,
the solutions obtained via the LBM and the FEM match exceptionally well. The
overall overshoot profiles are nearly identical, and the sudden jumps induced
by P-wave fronts in the LB and FE simulations also line up very well, with the
LBM producing slightly more pronounced fluctuations overall.

Figure 8: Tensile stresses postprocessed from LBM (blue) and FE (red)
simulations of a neo-Hooke solid under tension loading.

An animation highlighting the evolution of the deformation and the
$x_{2}$-component of the Cauchy stress throughout the simulation is available
from the authors upon request. Figure 9 additionally includes four exemplary
stills from this animation. The significance of the deformation – which is
shown to scale – is apparent. Additionally, the dilatational wave fronts
traversing the material domain are clearly visible. Along with the
considerable dynamic overshoots in the deformation, these give rise to a
complex stress field with a similarly complex time evolution.

Figure 9: Deformed configurations predicted by the LBM at $t\approx
0.7971l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$ (top left), $t\approx
1.232l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$ (top right), $t\approx
1.558l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$ (bottom left), and $t\approx
2.355l\sqrt{\nicefrac{{\rho_{0}}}{{\mu}}}$ (bottom right). A contour plot for
the $x_{2}$-component of the Cauchy stress $\sigma_{22}$ is superimposed on
the warped domain. Deformation to scale.

##  8 Discussion and conclusions

Generally, FEM and LBM are in better agreement with regard to the tensile
stress components than the shear stresses, both in the simple example and the
more challenging benchmark. Overall, the agreement is good, with very similar
peak displacement and stress values and profiles being predicted by the LBM
and the FEM. It is difficult to asses these results in the absence of
analytical benchmarks. On the other hand, the FEM is the state-of-the-art
method for nonlinear, transient, large displacement solid simulation. Because
of this, we think that it is encouraging that the LB algorithm presented above
can replicate FE results well. This is especially the case because the LBM for
nonlinear solid mechanics is still at an early stage of development and we had
to resort to artificial source terms to model part of the material law.

Additionally, the LBM seems to produce more pronounced fluctuations and jumps
in the computed stresses (and, albeit to a lesser extent, in the
displacements). This might be because the FEM tends to produce slight
numerical damping and more stable solutions, while the LBM damps solutions to
a lesser extent and is more prone to instabilities.

##  9 Summary and outlook

This work proposed a large displacement Lattice Boltzmann Method for
materially and geometrically nonlinear solid mechanics in the reference
configuration. To this end, the LBM was treated as a solver for balance laws
in the reference configuration, based on a recent numerical investigation of
the LBM by Farag et al. (2021). A moment chain for nonlinear solid mechanics
was formulated, and an LB algorithm proposed to solve it. Finally, results
produced by the LBM were validated against FEM reference solutions for two
simple examples and one slightly more challenging transient benchmark.

The agreement between the LB and FE solutions is encouraging, which suggests
that the LBM can indeed be used as a solver for nonlinear, transient, large
displacement solid mechanics in the reference configuration. This is exciting,
because the LB algorithm consists of computationally simple, highly
parallelisable operations. An efficient implementation of the solid-mechanical
LBM might have the potential to reduce simulation times for transient
simulations – for which computational efficiency is presently a bottleneck.
Furthermore, the simple spatial and temporal discretisation makes for easy
pre- and post-processing compared to FE workflows. Finally, LB schemes lend
themselves to coupled simulations for fluid-structure interaction: if LB
algorithms can be used on both sides of a fluid-solid interface, the required
coupling might be simplified considerably.

The extension to large displacements might open up new fields of application
for the LBM in sold mechanics: impact simulations and biomechanical tissue
modelling, for example, may eventually become feasible. Additionally,
nonlinear material models could allow for simulations of tires, seals, and
polymer components under transient load. The formulation in the reference
configuration, meanwhile, greatly simplifies the modelling of boundaries.

We are currently working on a three-dimensional implementation of the solid-
mechanical LBM to allow for the simulation of a wider range of interesting
phenomena. To improve the efficiency and stability of the method, we are
additionally looking for a cleaner implementation of the material law which
does not require source terms. Furthermore, we would like to extend the
algorithm presented in this work to inelastic material behaviour – such as
elasto-plasticity and viscoelasticity – and to use the LBM for forming
simulations.

## References

  * Farag et al. [2021] G. Farag, S. Zhao, G. Chiavassa, and P. Boivin.  Consistency study of lattice-boltzmann schemes macroscopic limit.  _Physics of Fluids_ , 33(3):037101, 2021.  doi:10.1063/5.0039490.
  * Murthy et al. [2017] J. Murthy, Praveen Kolluru, Vishwanathan Kumaran, and Santosh Ansumali.  Lattice boltzmann method for wave propagation in elastic solids.  _Communications in Computational Physics_ , 23, 08 2017.  doi:10.4208/cicp.OA-2016-0259.
  * Escande et al. [2020] Maxime Escande, Praveen Kumar Kolluru, Louis Marie Cléon, and Pierre Sagaut.  Lattice boltzmann method for wave propagation in elastic solids with a regular lattice: Theoretical analysis and validation.  _CoRR_ , abs/2009.06404, 2020.  URL https://arxiv.org/abs/2009.06404.
  * Schlüter et al. [2021] Alexander Schlüter, Sikang Yan, Thomas Reinirkens, Charlotte Kuhn, and Ralf Müller.  Lattice boltzmann simulation of plane strain problems.  _PAMM_ , 20(1):e202000119, 2021.  doi:https://doi.org/10.1002/pamm.202000119.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/pamm.202000119.
  * Krüger et al. [2017] Timm Krüger, Halim Kusumaatmaja, Alexandr Kuzmin, Orest Shardt, Goncalo Silva, and Erlend Magnus Viggen.  _The Lattice Boltzmann Method: Principles and Practice_.  Springer, Heidelberg, 1st edition, 2017.  doi:10.1007/978-3-319-44649-3.
  * Marconi and Chopard [2003] Stefan Marconi and Bastien Chopard.  A lattice boltzmann model for a solid body.  _International Journal of Modern Physics B_ , 17(01n02):153–156, 2003.  doi:10.1142/S0217979203017254.
  * O’Brien et al. [2012] G. S. O’Brien, T. Nissen-Meyer, and C. J. Bean.  A Lattice Boltzmann Method for Elastic Wave Propagation in a Poisson Solid.  _Bulletin of the Seismological Society of America_ , 102(3):1224–1234, 06 2012.  ISSN 0037-1106.  doi:10.1785/0120110191.
  * Schlüter et al. [2018] Alexander Schlüter, Charlotte Kuhn, and Ralf Müller.  Lattice boltzmann simulation of antiplane shear loading of a stationary crack.  _Computational Mechanics_ , 62(5):1059–1069, Nov 2018.  ISSN 1432-0924.  doi:10.1007/s00466-018-1550-4.  URL https://doi.org/10.1007/s00466-018-1550-4.
  * Müller et al. [2022] Henning Müller, Ali Touil, Alexander Schlüter, and Ralf Müller.  Dynamic propagation of mode iii cracks in a lattice boltzmann method for solids, 2022.  URL https://arxiv.org/abs/2206.12412.
  * Schlüter et al. [2022] Alexander Schlüter, Henning Müller, and Ralf Müller.  Lattice boltzmann method for antiplane shear deformation: non-lattice-conforming boundary conditions.  _Archive of Applied Mechanics_ , Aug 2022.  ISSN 1432-0681.  doi:10.1007/s00419-022-02240-z.  URL https://doi.org/10.1007/s00419-022-02240-z.
  * Faust et al. [2022] Erik Faust, Alexander Schlüter, Henning Müller, and Ralf Müller.  Dirichlet and neumann boundary conditions in a lattice boltzmann method for elastodynamics, 2022.  URL https://arxiv.org/abs/2208.04088.
  * Holzapfel [2000] A Gerhard Holzapfel.  _Nonlinear solid mechanics: A continuum approach for engineering_.  John Wiley & Sons, Inc., Chichester, 1st edition, 2000.  doi:10.1023/A:1020843529530.
  * He et al. [1998] Xiaoyi He, Shiyi Chen, and Gary D. Doolen.  A novel thermal model for the lattice boltzmann method in incompressible limit.  _Journal of Computational Physics_ , 146(1):282–300, 1998.  ISSN 0021-9991.  doi:10.1006/jcph.1998.6057.  URL https://www.sciencedirect.com/science/article/pii/S0021999198960570.
  * Taylor [2014] Robert L Taylor.  Feap-a finite element analysis program, 2014.