## 3 Numerical Implementations and Algorithms
### Additional Solid Routines
For instance, around internal lattice nodes, a central difference scheme with a \(\mathcal{O}(\Delta x^{4})\) truncation error is built, see Equation (13). The first derivative along the first spatial direction X of any function \(\mathcal{Q}\) with discrete values over the lattice can be expressed as:

\[\frac{\partial Q(G_{x},G_{y})}{\partial\texttt{X}}=\frac{Q(G_{x}-2,G_{y})-8Q( G_{x}-1,G_{y})+8Q(G_{x}+1,G_{y})-Q(G_{x}+2,G_{y})}{12\Delta x}+\mathcal{O} \Big{(}\Delta x^{4}\Big{)}\,, \tag{13}\]

where \((G_{x},G_{y})\) denotes a specific lattice node in lattice global coordinates. Since this previous scheme uses a two-neighborhood environment, application near boundaries cannot be applied in this study. ### Proposed Solid Equilibrium and Colliding Step

As we aim to study static equilibrium state, we suggest using Equation (9), in order to relax towards a solid equilibrium. Moreover, the use of the Bhatnagar, Gross and Krook (BGK) collision operator, enhances the numerical stability and avoid some unwanted mechanical waves or issues near boundaries. This last point is further developed in the next Section 3.6. However, our BGK collision operator for solids differs by its equilibrium distribution from the usual one used in CFD. We tried to reach the static equilibrium through the BGK operator considering an equilibrium distribution without macroscopic speed. A simple distribution that does not depend on temperature, and ideally where matter does not move at the microscopic scale level, is targeted. However, for arithmetic reasons, mass is not concentrated on the cell's first discretized distribution function. The following numerical tests proved the functionality of the proposed distribution functions at solid equilibrium. From there, for each cell composing the lattice, distribution at solid static equilibrium for a classic \(D2\eta\) scheme can be defined by (see Equation (9)):

\[f_{i}^{(eq),sld}(\rho(x,t),\
u)=f_{i}^{(eq)}(\rho(x,t),0)=w_{i}\rho(x,t)= \begin{cases}\frac{4}{3}\rho(x,t)&\text{if }i=0,\\ \frac{1}{3}\rho(x,t)&\text{if }i\in[1,\ldots,4],\\ \frac{1}{36}\rho(x,t)&\text{if }i\in[5,\ldots,8],\end{cases}  \tag{14}\]

where \(f_{i}^{(eq),sld}\) represents a discrete value \(i\) of \(f^{(eq),sld}\) for a cell node \((G_{x},G_{y})\). The presented distribution is applied both for initial conditions and equilibrium during time iterations. Hence, the proposed equilibrium only varies in space and time through density, and does not depend on the solid constitutive law chosen. Thus, mass conservation is assured at solid equilibrium thanks to the summation of Gauss-Hermite weights \(w_{i}\). One can also note that in Equation (14) the common choice of \({c_{s}}^{2}=\frac{1}{3}\) is chosen. However, in this special case of solid static equilibrium, the equations are of diffusion type. In such a case, the sound celerity can be chosen freely since the BGK or Two Relaxation Times (TRT) operators are unconditionally stable [40]. The colliding step (see Equation (2)) is then modified to include the above suggested solid equilibrium and related force calculated from stress divergence. For instance, considering an internal node, not on boundaries, the discretized collision operator is expressed by:

\[\Omega\Big{(}f_{i}(G_{x},G_{y},t),f_{i}^{(eq),sld}\Big{)}=-\omega\Big{(}f_{i}(G_{ x},G_{y},t)-f_{i}^{(eq),sld}\Big{)}+s_{i}(G_{x},G_{y},t),\forall i\in[0,\ldots,8]. \tag{15}\]

Observe that \(f_{i}(G_{x},G_{y},t)\) and \(f_{i}^{(eq),sld}\) depends on space and time whereas \(f_{i}^{(eq),sld}\) is also dependent on lattice scheme properties. Moreover, \(s_{i}(G_{x},G_{y},t)\) represents the lattice source term related to a body force for a node \((G_{x},G_{y})\) at time \(t\) applied on distribution, see Algorithm 3. \(s_{i}(G_{x},G_{y},t)\) is also scheme dependent because its computation is based, for instance, on Gauss-Hermite weights \(w_{i}\). Different force implementations exist for LBM. All the following illustrated results are made using the Guo approach [35], see Equations (2) and (7).