5.1. Obtained Results and Discussion
In addition to the global curves matching, which gives a qualitative validation, a quantitative analysis is provided. Thus,
Table 1 shows the Root Mean Square Error (RMSE) and Normalized Root Mean Square Error (NRMSE) between the computed results with the LBMS and FEA, presented through
Figure 4,
Figure 5 and
Figure 7. The results show a great precision of the LBMS compared to the FEA method, and this is for a large number of study cases. This precision illustrates the fidelity of the LBMS. The numerical errors obtained are far lower than the amplitude of the computed fields, which yields only a few percents of error. These errors values tend to quantitatively validate the presented approach.
The numerical tests highlight the stability and accuracy of the suggested methods in a very broad range of mechanical parameters. This large range allows us to conclude that the current method could be applied to numerous constitutive laws, and not only to the isotropic linear elasticity. This wide range tends to support the theoretical development presented in this study, especially the incorporation of the stress divergence as a mean-field force.
Concerning the introduced boundary conditions for solids, some work remains in order to enhance adaptation for solids. Currently, it has been pointed out that fixed conditions are far from perfect and suitable for accurate computation near concerned areas: impreciseness on fields can arise due to this modeling issue.
Light differences with reference plots may have several explanations. First and foremost, the LBMS method detains an iterative process so as to arrive at solid static equilibrium; thus, the finite difference scheme accuracy could play an important role concerning quality of the calculated results. Moreover, we point out some post-processing issues which may be at the origin of the errors. Indeed, along cutting curves, the FEA displacement field is mainly interpolated whereas the LBMS field hardly is. All of these considerations can cause additional deviations from the reference study and may explain some of illustrated differences.
Other slight errors may be generated by the non-differentiation between Eulerian and Lagrangian frameworks. Several iterations are needed to reach the equilibrium state (i.e., when stress divergence tends towards a null vector), yielding completely different pipelines between FEA and LBMS. Furthermore, as mentioned above in
Section 4.1.4, the used plane stress assumption for two-dimensional reference studies might trigger some of the shown deviations in curve plots, especially for non-zero Poisson’s ratios. Indeed, thin plates are usually modeled with a plane stress hypothesis, but this approach is softly inconsistent with the one used to produce illustrated LBMS runs, i.e., a plane strain assignment. For all of these reasons, extra investigations seem to be mandatory to fit calculated reference fields with better accuracy.
All these achievements seem promising in terms of the use of the Boltzmann–Vlasov approach for solid dynamics. It is also worth noting that the mass of the system is always conserved. Indeed, the intrinsic advantages of the LBM are preserved by the LBMS (macroscopic equations recovered from the statistical moments, locality, simplicity, etc.). In contrast to previous attempts to cope with solid behavior with the LBM, our approach does not require statistical distribution of forces [
18,
45] nor springs [
24]. Both bring more complexity, heavier computations, are further away from the Boltzmann theory, and are less prone to FSI coupling.
For our validation case, we do not have an analytical solution due to boundary conditions issues. A simpler test case is not easily achievable in a two dimensional domain. Concerning numerical correlations, the standard FEA solution is completely approved by the mechanical community and hence serves as a reference solution. We further point out that both methods (i.e., LBMS and FEA) are truly different and must be compared carefully, see
Table 2 for more details. Our proposed LBMS method is not very competitive in terms of computational time, because it was developed to evolve to solve complex FSI dynamic problems. In these mechanical cases, FEA runs are time-consuming calculations and the previously explained trend could be drastically reversed. On the other hand, the LBM approach can trade its “locality” for an improved computational time in steady-state cases as it has been illustrated in [
40]. However, because our focus is on FSI dynamic problems, the latter approach was not used here. However, it contributes to the explanation of the limited computational time competitiveness of the suggested method.
5.2. Considerations about the Vlasov Equation
Despite the previous attempts of LBM for solids mentioned in
Section 1.1, to the best of the authors’ knowledge, none of them use the physical representation of the particles of the Boltzmann equation. Indeed, the method previously used to derive an LBM–like method is considering the LBM equation as a pure mathematical framework to solve a set of Partial Differential Equations (PDE)s, without relying on the original fluid-based historical method. These previous mathematical methods found a set of advection–relaxation equations for new variables that are equivalent to a given macroscopic model. Such mathematical approaches contrast with the suggested method. We propose keeping this link between particles and the Boltzmann method. To do so, we find the similarities between plasma and solids. In both cases, the motion of particles is almost not subject to collisions: in the plasma case, the Knudsen number is large
and particles are collisionless in the free-flight regime; whereas in solids
and particles are collisionless, because a particle does not move freely without interacting with the surrounding particles. Even if these two regimes should not be confused, they share some physical and mathematical concepts resulting in the fact that in both cases the collision term in the Boltzmann equation vanishes. Such a collisionless Boltzmann equation is called a Vlasov equation [
46].
In order to find a suitable framework for the use of the Boltzmann equation for solids, it is possible to note that, at rest (at thermodynamic equilibrium), a particle distribution having a low temperature could be described by a Dirac distribution (omitting all quantum effects). Then, the solid matter at rest condition can be seen as a combination of motionless point mass, and without a large energy perturbation, the molecular interactions can be linearized. Starting from this state, a small vibration of one edge of the network will propagate freely. This situation can be described by the Vlasov equation of Dirac distribution centered on the velocity of the the surrounding particles, . Then, the evolution of the particles described by this equation is a propagation of the particles without further modifications induced by the collisions.
The interesting idea behind the Vlasov equation for our purpose is that complex systems with few collisions and strong interactions can be modeled by a Vlasov equation. Plus, these strong interactions are incorporated through the mean-field external forces term, which is the concept used in Equation (8). Therefore, considering the Vlasov equation of Dirac distribution with a mean-field divergence of Cauchy’s stress, we obtained the following mesoscopic model:
Then, using the passage formula for the two firsts moments (
), the desired macroscopic system of the undefined equation of motion for solids is obtained from a theoretical point of view:
Thus, the Boltzmann–Vlasov equation leading to the previous system of equations seems suitable for rigid matter simulations.
Several remarks arise from these previous equations. This approach remains general for any constitutive law. Thus, the use of the Dirac distribution is somehow closer to the conventional representations of solid lattices: particles are not moving freely through the solid. In addition, the current discussion about the Vlasov equation for solids, is certainly presenting a working framework; however, the founding of this Vlasov equation requires a mathematical proof which is left for future work. Instead, we test here numerically some of the assumptions relative to this framework that are restricted to static equilibrium.
5.3. Global Outlook and Proposed Improvements
C++ programmed classes and methods integrated into LBMS library can be drastically enhanced, i.e., mostly for efficiency purposes. Although the library is robust, computation time might be substantially reduced if other guidelines are adopted.
For presented cases, tens or hundreds of seconds are necessary to reach the illustrated static equilibrium assuming sequential runs on an Intel(R) Xeon(R) CPU E5-1650 v3 @ 3.50 GHz, thus, involving a few tens of thousands of LBMS iterations. In contrast, a couple of seconds are enough to COMSOL even if the mesh is modeled very fine, because a standard static resolution is requested. Effectively, the involved computation times are hardly comparable because the methods are very different. In fact, LBMS is naturally dynamic due to its roots coming from the LBM and thus it requires an evaluation of the solid state at each time iteration. Concerning lattice discretization, additional experiments should be conducted in order to discover space convergence properties and effects.
Thanks to the pointed reproducibility of the depicted numerical results, complements to existing objects can be made freely in order to improve calculation efficiency. However, here stands one of the great strengths of the LBM: parallel computing can be seen as one of the main advantages of LBM. The LBM is intrinsically efficient and adapted to massive parallel computing, especially on Graphics Processing Unit (GPU) because of its architecture. In fact, a GPU detains hundreds or thousands of cores capable of doing very fast executions despite management of few instructions and small caches.
Finite difference schemes are playing an important role; in fact, they define external forces to be applied at each iteration for each cell. Upgrading abilities to treat mechanical singularities and enhancing accuracy under high field variation are part of the future work. For example, corners demand a special therapy because cell vicinity is restricted. Another possible way is to use robust mean-field forces as already implemented into the Shan-Chen approach [
10]. Moreover, as mentioned above, WENO schemes could help us.
Another advantage of the LBM compared to FEA, is this ease of moving from lower to higher dimensions. The transition from the 2D modeling to the 3D and 4D modeling is implied in the method. Then, LBMS extension to 3D should be straightforward and large-scale studies can be led from now on. Such computations permit a better comparison of the complex geometries with current methods usually dedicated to solving mechanical problems. The door is thus wide open in terms of further intensive analysis or significant studies. In addition, the parallelized strategy by blocks remains unchanged for solids. However, depending on the used finite difference schemes, some complementary properties are needed for block envelopes.
Furthermore, for the broad range of Young’s modulus and Poisson’s ratios simulated, the presented method must be compared with non-linear constitutive laws, to confirm its potentiality and generality. Based on the strength of the presented results, a very straightforward extension is to investigate quasi-static dynamic problems. Of course, this eventuality could only be viable if the computation time of the LBMS is more deeply characterized.
Only the static equilibrium has been compared to existing methods. Because this static equilibrium has been reached for different relaxation time, in particular for ; solid dynamics with the Boltzmann–Vlasov equation seems a natural next step. The numerical stability for the relaxation parameter is already unconventional LBM and supports the theoretical developments introduced here.
Validation of the dynamic behavior requires additional work, especially for the reference mechanical study. Indeed, the transitional state before reaching static equilibrium is hard to confirm taking into account the actual state of our research and LBM groundwork for solids.
As we mentioned before, the goal of the solid behavior with the LBM, might also be investigated through the research of an appropriate equilibrium distribution function. To go further with that idea, if
, we also observe that (see Equation (15)):
This latter formulation yields a new form of solid equilibrium where solid body forces are integrated into the equilibrium function; thus, the stress divergence term may be integrated into the BGK operator if the relaxation parameter is equal to 1. Moreover, such an approach may looks like an artifact because with , the LBM is equivalent to a particular Finite Difference Method (FDM), with more complex boundary conditions. Thus, the previous equilibrium distribution deserves further dedicated investigations, beyond the work.
Of course, many applications are plausible from the next development of the LBMS, especially regarding the interface and boundary conditions. By solving the solid interface, the motion and deformation of solids under the action of forces and contact problems between two solids could be achieved. Then, the solid–liquid interface leads to FSI in order to study the behavior of a solid product immersed in a liquid. A recent publication [
47], illustrates the fact that a unified theory is needed to solve complex FSI phenomena, taking into consideration both fluid and solid characteristics with LBM.