SPH modelling of hydrodynamic lubrication: laminar fluid flow-structure interaction with no-slip conditions for slider bearings

The FOSS CFD-SPH code SPHERA v.9.0.0 (RSE SpA) is empowered to deal with fluid-solid body interactions under no-slip conditions and laminar regimes for the simulation of hydrodynamic lubrication. The code is herein validated in relation to a uniform slider bearing (i.e., for a constant lubricant film depth) and a linear slider bearing (i.e., for a film depth with a linear profile variation along the main flow direction). Validations refer to comparisons with analytical solutions, herein generalized to consider any Dirichlet boundary condition. Further, this study allows a first code validation of the fluid-fixed frontier interactions under no-slip conditions. With respect to the most state-of-the-art models (2D codes based on Reynolds' equation for fluid films), the following distinctive features are highlighted: (i) 3D formulation on all the terms of the Navier-Stokes equations for incompressible fluids with uniform viscosity; (ii) validations on both local and global quantities (pressure and velocity profiles; load-bearing capacity); (iii) possibility to simulate any 3D topology. This study also shows the advantages of using a CFD-SPH code in simulating the inertia and 3D effects close to the slider edges, and it opens new research directions overcoming the limitations of the codes for hydrodynamic lubrication based on the Reynolds' equation for fluid films. This study finally allows SPHERA to deal with hydrodynamic lubrication and empowers the code for other relevant application fields involving fluid-structure interactions (e.g., transport of solid bodies by floods and earth landslides; rock landslides). SPHERA is developed and distributed on a GitHub public repository.


Introduction
In tribology, the science and technology of interacting solid surfaces under relative motion, three lubrication regimes are defined: the full-film lubrication regime (i.e. the solid surfaces are completely separated by the fluid), the boundary lubrication regime (i.e. the solid surfaces are in direct contact by means of their asperities and possible additives) and the mixed lubrication regime  [3] found that the assumption of no advection/inertia terms is responsible for up-stream under-estimations (around the leading edge) and down-stream overpredictions (around the trailing edge) in the pressure field for step bearings.
Almqvist et al. [4] provided inter-comparisons between a FD (Finite Difference) model based on Reynolds' equation for fluid films and a commercial CFD-FVM code based on Navier-Stokes equations.In a series of follow-up articles, Almqvist et al. [5,6] reported analytical solutions on pressure longitudinal profiles, velocity vertical profiles, friction force and load-bearing capacity (the frictional coefficient is the ratio between these two forces) for both a linear slider and the Rayleigh 3 step slider, with null Dirichlet boundary conditions for pressure. They also derived the optimum geometric configuration for a linear slider to maximize the load-bearing capacity (k=1.189, where k is the distance increase of the surface gap per unit length). Further, they analysed the effects of surface roughness by means of a homogenization technique.
Rahmani et al. [7] presented an analytical approach based on Reynolds' equation for asymmetric partially textured slider bearings with surface discontinuities, to optimize the choice of the textures parameters with respect to the load-bearing capacity and the friction force. Papadopoulos et al. [8] use a 2D CFD-FVM code to optimize micro-thrust bearings with surface texturing by means of numerical inter-comparisons.Fouflias et al. [9] used a commercial CFD-FVM code to simulate bearings with pockets/dimples and surface texturing, providing model inter-comparisons on steadyloads for different designs.
Regarding modelling of complex surface topologies, Paggi and He [10] investigated the effect of roughness on the evolution of the channel network influencing the fluid flow in the mixed lubrication regime.
Gropper et al. [11] proposed a detailed review on hydrodynamic lubrication of textured surfaces, included (multi-scale) roughness effects and cavitation. Hajishaflee et al. [12] adopted a 2D CFD-FVM model to reproduce elasto-hydrodynamic lubrication problems for rolling element bearings, including cavitation effects. Snyder and Braun [13] provided analytical solutions and 2D CFD -FVM (Finite Volume Method) results to quantify the squeezing effects on a sliding bearing load in terms of dynamic coefficients (dimensionless stiffness and damping). They also used a PRE (Perturbed Reynolds Equation) -FD model, which is based on the representation of perturbed quantities (film thickness and pressure) within Reynolds' equation, and provided three separated differential equations for static pressure, dynamic pressure associated with stiffness, and dynamic pressure associated with damping.
With respect to the state-of-the-art on CFD modelling for bearings, mostly based on 2D codes or Reynolds' simplified equation, the present study adopts a 3D CFD discretization of all the terms of 4 the Navier-Stokes equations for incompressible fluids with uniform viscosity. It also provides validations on local quantities (pressure and velocity profiles) and it is able to take into account any 3D surface in input. In particular, the present study empowers, validates and applies the FOSS (Free/Libre and Open Source Software) CFD-SPH (Smoothed Particle Hydrodynamics) code SPHERA [14] to "fluidsolid body" interactions under no-slip conditions, simulating uniform and linear sliders. Validations are provided by comparisons with analytical solutions, here generalized to deal with non-null Dirichlet boundary conditions for pressure. Validations refer to pressure longitudinal profiles, velocity vertical profiles and load-bearing capacity. Furthermore, a demonstrative test case is simulated to represent a linear slider moving over a complex 3D surface, to show its applicability to any 3D surface input data. Beyond the new numerical developments, other code features are herein first validated, in particular the "fluidfixed frontier" interactions under no-slip conditions. This study represents one the first applications of the SPH method to bearings. The basic features of this numerical method are briefly recalled hereafter.
Smoothed Particle Hydrodynamics (SPH) is a mesh-less CFD method, whose computational nodes are represented by numerical fluid particles. In the continuum, the functions and derivatives in the fluid dynamics balance equations are approximated by convolution integrals, which are weighted by interpolating (or smoothing functions), called kernel functions.
The integral SPH approximation (<>I) of a generic function (f) is defined as: where W (m -3 ) is the kernel function [15], x0 (m) the position of a generic computational point and Vh (m 3 ) the integration volume, which is called kernel support. This is represented by a sphere of radius 2hSPH (m), possibly truncated by the frontiers of the fluid domain.
Any first derivative of a generic function, calculated along i-axis, can be computed as in (1), after replacing f with the targeted derivative. After integration by parts, one obtains: The integration also involves the surface Ah (m 2 ) of the kernel support. The associated surface integral is non-zero in case of a truncated kernel support. The representation of this term noticeably differentiates SPH codes among each other [16][17][18][19][20].
Far from boundaries, the SPH particle approximation of (2) reads: where a summation on particle volumes  (m  ) replaces the volume integral. The subscripts "0" and "b" refer to the computational particle and its "neighbouring particles" (fluid particles within the kernel support of the computational particle), respectively.
Usually, the approximation (3) is replaced by more complicated and accurate formulas. Further, the SPH method can also approximate a generic n-th derivative, analogously to (3).
Among the various numerical methods, Smoothed Particle Hydrodynamics (SPH) has several advantages: a direct estimation of free surface and phase/fluid interfaces; effective simulations of multiple moving bodies and particulate matter within fluid flows; direct estimation of Lagrangian derivatives (absence of non-linear advective terms in the balance equations); effective numerical simulations of fast transient phenomena; no meshing; simple non-iterative algorithms (in case the "Weakly Compressible" approach is adopted). On the other hand, SPH models are affected by the following drawbacks, if compared with mesh-based CFD tools: computational costs are slightly higher due to a larger stencil (around each computational particle), which causes a high number of interacting elements (neighbouring particles) at a fixed time step (nonetheless SPH codes are more suitable to parallelization); local refining of spatial resolution represents a current issue and is only addressed by few, advanced and complex SPH algorithms; accuracy is relatively low for classical CFD applications where mesh-based methods are well established (e.g., confined mono-phase flows). Detailed reviews on SPH assets and drawbacks are reported in [21][22][23][24]. Nevertheless, SPH 6 models are effective in several, but peculiar, application fields. Some of them are here briefly recalled: flood propagation (e.g., [25,26]); sloshing tanks (e.g., [20]); gravitational surface waves (e.g. [27]); hydraulic turbines (e.g., [28]); liquid jets (e.g., [28]); astrophysics and magnetohydrodynamics (e.g., [29]; body dynamics in free surface flows (e.g., [30]); multi-phase and multifluid flows; sediment removal from water reservoirs (e.g., [31]); landslides (e.g., [32,33]).
The paper is organized as follows. Sec. 2 reports the state-of-the-art mathematical models and analytical solutions relevant for this study. Sec. 3 derives a generalization of the analytical solutions for the linear slider to take into account non-null Dirichlet boundary conditions for pressure. Sec. 4 reports those basic features of the reference CFD-SPH code relevant for this study, whereas Sec. 5 introduces the new numerical developments of the code. Validations refer to a uniform slider (Sec. 6) and a linear slider (Sec. 7). A demonstrative test case for a linear slider with a 3D complex surface is briefly outlined in Sec. 8. The overall conclusions are synthesized in Sec. 9.

Benchmark analytical solutions
One considers a linearization of Navier-Stokes' equations for incompressible Newtonian fluids, with a fluid film flowing between two solid plates: where p is pressure, u is the velocity vector,  is the dynamic viscosity, ij is Kronecker's delta and x represents a generic position. Advective and gravity terms are herein neglected, and the viscous shear stress terms only affect the horizontal projection of momentum. Combining the above expressions, defining the fluid depth h (x,y,t) and assuming the following hypotheses (h0 represents the minimum value of h and L the upper plate length): one obtains Reynolds' equation for fluid films: where the subscripts "s1" and "s2" denote the upper and lower plate, respectively. The last formula represents a 2D time-dependent equation to describe the dynamics of fluid films between two solid plates.
In case of stationary regime and uniformity along the y-axis, the above equation reduces to 1D Reynolds' equation for fluid films: In case of uniform viscosity, one obtains: Under these assumptions, the analytical solution of any velocity profile for generic plate geometries assumes the following form: where velocity is proportional to the pressure derivative along x-axis and the mass flow rate is uniform [6]. The analytical solutions in [6] for both the uniform slider and the linear slider are reported in Sec. 2.1Errore. L'origine riferimento non è stata trovata. and Sec. 2.2, respectively.

Uniform slider
The uniform slider is featured by a uniform value for the fluid depth h. This geometrical configuration implies a 1D Laplace equation for pressure, which assumes a linear horizontal profile: where the subscripts "L" and "0" denote the plate edges.
Considering Eq. (9), the analytical solution for the velocity profiles within a uniform slider is: The load-bearing capacity lc represents the hydrodynamic thrust exerted on the plate. On the uniform slider, the analytical solution for the load-bearing capacity (per unit width) reads: The following non-dimensional quantities are defined for pressure, velocity magnitude and time: with Cp denoting the pressure coefficient.
The non-dimensional load-bearing capacity LC is defined as follows: For a uniform slider, the analytical solution for LC is equal to the average of the pressure coefficient values provided as boundary conditions:

Linear slider
Within a linear slider, the depth rate of change of the lubricant along the x-axis is constant: Almqvist [6] assumed null pressure at the edges of the linear slider and provided analytical solutions for the longitudinal profile of pressure (which is uniform along the vertical direction), the vertical profiles of the horizontal velocity, and the load-bearing capacity: (17) Assuming a vanishing pressure at the edges of a linear slider is a theoretical simplification, which does not take into account the over-pressure and the under-pressure zones determined by the interaction between the fluid flow and the solid plate at its leading and trailing edges, respectively.

Generalized analytical solutions for a linear slider
The analytical solutions of Almqvist et al. [6] are herein generalized to impose non-null Dirichlet's boundary conditions for pressure at the edges of a linear slider, to cope with more general and practical configurations. No matter about the particular value of pressure at the edges of the plate, the longitudinal pressure profile reads [6]: where C1 and C2 are integration constants.
The velocity scale U is the x-component of the vector summation of the velocities of the solid plates: One imposes generic now non-null Dirichlet's boundary conditions for pressure at the slider edges: Solving the first formula of (20) for C2 and replacing this expression in the second formula of (20), one obtains the following expression for the integration constant C1: Replacing C1 into the second expression of (20), one obtains the integration constant C2: Replacing the expressions for the integration constants (21)- (22) and the fluid depth (16) into Eq.
(18), the pressure profile assumes the following form: After minor arrangements, one obtains the following analytical solution for the longitudinal profile of pressure within a linear slider under generic non-null Dirichlet boundary conditions: It is remarkable to note that the first term of the Right Hand Side (RHS) of Eq. Under normal conditions (pressure at the leading edge "0" is higher than pressure at the trailing edge "L"), the second and third terms of the RHS of Eq. (24) raise pressure levels. The second term also shifts the pressure peak towards the body edge with higher pressure.
The integration over the plate length of the analytical solution for pressure (Eq. (24)) provides the expression for the load-bearing capacity: The second to last term in Eq. (25) assumes the following expression: And, once it is introduced back into Eq. (25), the load-bearing capacity reads: and its non-dimensional formulation renders:

12
A generic analytical form for the 2D field of the x-component of velocity (u) is expressed by Almqvist et al. [14]: Considering the integration constant C1 as provided by Eq. (21), the generalized analytical solution for the 2D velocity field assumes the following expression: The configuration of the linear slider used to validate the code SPHERA (Sec. 6.

Description of the SPH computational framework
The reference code SPHERA used in this study [14] has been applied to floods (with transport of solid bodies, bed-load transport and a domain spatial coverage up to some hundredths of squared kilometres), fast landslides and wave motion, sediment removal from water reservoirs, sloshing tanks. SPHERA is featured by: a scheme for dense granular flows [34]; a scheme for the transport of solid bodies in free surface flows [30]; a scheme for a boundary treatment ("DB-SPH") based on discrete surface and volume elements, and on a 1D Linearized Partial Riemann Solver coupled with a MUSCL (Monotonic Upstream-Centered Scheme for Conservation Laws) spatial reconstruction scheme [20]; a scheme for a 2D erosion criterion [31]; a scheme for a boundary treatment ("semianalytic approach or SA-SPH" for simplicity of notation) based on volume integrals, numerically computed outside of the fluid domain [35].

SPH approximation of the balance equations for fluid dynamics and the boundary treatment scheme called "semi-analytic approach"
The numerical scheme for the main flow is a Weakly-Compressible (WC) SPH model, which takes benefit from a boundary treatment for fixed boundaries based on the semi-analytic approach of Vila [36], as developed by Di Monaco et al. [35].
One considers Navier-Stokes' momentum and continuity equations for incompressible fluids with uniform viscosity: 14 One needs to compute Eq. (32) at each fluid particle position by using the SPH formalism and taking into account the boundary terms (fluid-frontier and fluid-body interactions), as described in the following.
One considers the discretization of Eq. (32), as provided by the SPH approximation of the first derivative of a generic function (f), according to the semi-analytic approach [36] ("SA"): The inner fluid domain here involved is filled with numerical particles. At boundaries, the kernel support is (formally) not truncated because it can partially lie outside the fluid domain. In other words, the summation in Eq. (33) is performed over all the fluid particles "b" (neighbouring particles with volume ) in the kernel support of the computational fluid particle ("0"). At the same time, the volume integral in Eq. (33) represents the boundary term, which is a convolution integral on the truncated portion of the kernel support. In this fictitious and outer volume Vh' (m 3 ), one needs to define the generic function f (pressure, velocity or density alternatively).
The semi-analytic approach ("SA") (as developed in [35]) introduces the following linearization and assumptions to compute f in Vh': The peculiar "SA" values of the functions and derivatives within Vh' are assigned to represent a null normal gradient of reduced pressure at the frontier interface (while considering uniform density): ; , represented as a linear extrapolation from the computational fluid particle velocity. The latter is 15 equal to its analogous vector of the same computational fluid particle (the subscript "w" refers to a generic frontier), in case of no-slip conditions: The velocity field within the truncated portion of the kernel support assumes is defined as follows: where n is the normal vector of the wall surface, according to its local orientation.
At this point, one can write the continuity equation for a Weakly Compressible SPH model (Einstein's notation works for the subscript "j"), using the semi-analytic approach for the boundary integral term (second term on the Right Hand Side): where Cs (kg×m -3 ×s -1 ) represents a "fluid-body" interaction term.
On the other hand, one can analogously derive the approximation of the momentum equation (the notation " " indicates the SPH particle -discrete-approximation): where as (m×s -2 ) is a the acceleration term due to the fluid-body interactions, M (m 2 ×s -1 ) is the artificial viscosity [15], m (kg) the particle mass and r (m) the relative distance between the neighbouring and the computational particle. 16 The last two terms of the RHS of Eq. (39) represent the inner and boundary components of the shear stress gradient term. They are presented in [35] but had never been validated on a test case where their role was non-negligible. This study allowed fixing several bugs in the 3D implementation of these terms which could not be used with the release SPHERA v.8.0 [14].
Finally, a barotropic equation of state (EOS) is linearized as follows: The artificial sound speed c (m/s, "ref" stands for a reference state) is usually 10 times higher than the maximum fluid velocity (WC approach) to obtain a density relative error within 1%. However, Monaghan [15] demonstrated this statement by assuming that the maximum pressure coefficient is 2. Each test case of this study considers an external force continuously applied to the fluid so that the maximum value of Cp is noticeably higher than 2. A generalization for the definition of the sound speed value is here obtained assuming that c is 5Cp,max times higher than the maximum fluid velocity.
The reader is referred to [30] and in [35] for more details.

SPH balance equations for rigid body transport
Body dynamics is ruled by Euler-Newton equations, whose discretization takes advantage from the SPH formalism and the coupling terms derived in the following sections: Here the subscript "B" refers to a generic computational body and "CM" to its centre of mass. In this sub-section, r implicitly represents the relative distance from the body centre of mass.
In order to solve the system (41), we need to model the global force and torque, as described in the following.
The resultant force is composed of several terms: where G (kg×m×s -2 ) represents the gravity force, whereas PF (kg×m×s -2 ) and TF (kg×m×s -2 ) the vector sums of the pressure and shear forces provided by the fluid. Analogously, PS (kg×m×s -2 ) and TS (kg×m×s -2 ) are the vector sums of the normal and the shear forces provided by other bodies or boundaries (solid-solid interactions). In case of inertial and quasi-inertial fluid flows, turbulence schemes and tangential stresses are not mandatory (simplifying hypothesis).
The fluid-solid interaction is expressed by the hydrodynamic thrust: The computational body is numerically represented by solid volume elements, here called (solid) "body particles" ("s"). Some of them describe the body surface and are referred to as "surface body particles". These particular elements are also characterized by an area and a vector n of norm 1, is perpendicular to the body face of the particle (it belongs to) and pointing outward the fluid domain (inward the solid body). The pressure of a body particle is computed as described in Sec.4.3.
The torque in (41) is discretized as the summation of each vector product between the relative position rs, of a surface body particle with respect to the body centre of mass, and the corresponding total particle force: Time integration of (41) is performed using a Leapfrog scheme synchronized with the fluid dynamics balance equations. This means that the body particle pressure is computed simultaneously to the fluid pressure, so that this parameter is staggered of around dt/2 with respect to all the other body particle parameters.
After time integration, the model provides the velocity of a body particle as the vector sum of the velocity of the corresponding body barycentre and the relative velocity:

Fluid-body interaction terms
The fluid-body interaction terms rely on the boundary technique introduced by Adami et al. [16], as implemented and adapted for free-slip conditions by Amicarelli et al. [30]. If the boundary is fixed, then this method can be interpreted as a discretization of the semi-analytic approach used to treat fluid-boundary interactions (Sec.4.1). The outer domain of () is herein represented by all the body particles inside the kernel support of the computational fluid particle. Further, Adami et al. [16] introduced a new term, related to the acceleration of the fluid-solid interface, which influences the estimation of body particle pressure.
The fluid-body interaction term in the momentum equation (39) assumes the form: More details are available in Amicarelli et al. [30]. The new formulation for the pressure of a generic body particle under no-slip conditions is presented in Sec.5.1.

Time integration schemes
Time integration is ruled by a second-order Leapfrog scheme (stability analysis and time integration schemes in SPH modelling are discussed in [37], as described in [35] and [30]: where CFL is the Courant-Friedrichs-Lewy number. Following [16], the viscous term stability parameter is set to C=0.05. 20

New developments for the mathematical and the numerical models: shear stress gradient term and no-slip conditions for fluid-body interactions
This section describes the new numerical developments implemented into the code SPHERA [14] by RSE SpA to describe the shear stress gradient term and no-slip conditions at the interface between the liquid domain and a mobile solid body (i.e. fluid-body interactions). The new coupling terms are introduced in the momentum (Sec. 5.1) and the continuity (Sec. 5.2) equations.

Fluid-body coupling terms for the momentum equation
Modelling the shear stress gradient term in the Navier-Stokes equation requires the introduction of an additional fluid-body coupling term in the RHS of the momentum equation. The discretization of this coupling term, which takes into account the shear stress exchanged at the fluid-body interface, needs to be coherent with the SPH particle approximation in the inner domain and assumes the following form: where the subscript "s0" denotes a generic solid-fluid inter-particle interaction and the solid particle volume is computed as follows: The inter-particle velocity us0 represents the field of the fluid velocity virtually reconstructed within the portion of the kernel support, which is truncated by the solid body.
The component of the inter-particle velocity, which is normal to the interface, guarantees no mass penetration (symmetric conditions) at the interface: Under no-slip conditions, the component of the inter-particle velocity, which is tangential (subscript "T") to the interface guarantees a uniform velocity gradient around the interface (if its position is assumed to be the average of the positions of the interacting particles): where the unit vector t is tangential to the interface.
Under no-slip conditions, the difference between the inter-particle velocity and the fluid particle velocity in Eq. (51) is expressed as follows: No-slip conditions also affect the formulation of the pressure gradient coupling term. The pressure value of the generic neighbouring surface body particle "s" in Eq. (48) depends on the particular computational fluid particle "0" we are considering, so that we can refer to the interaction subscript "s,0".
The pressure value of the generic neighbouring surface body particle "s" is derived as follows.
Consider a generic point at a generic fluid-body interface. In case of free-slip conditions, the normal projection of the acceleration on the fluid side ("f") and on the solid side ("w") are equal (in-built motion in the direction normal to the interface): w w w f f w f n a n g p n dt The "wall" acceleration at the position of a generic body particle can then be derived by linearizing Eq. (56). This depends on the particular computational fluid particle "0" we are considering, so that we can refer to the interaction subscript "s,0": where dln is a vectorial length element along the centreline of the two particles, projected along the wall element normal: One applies a SPH interpolation over all the pressure values estimated according to [16] to derive a unique pressure value for a body particle. In case of no-slip conditions, Eq. (58) assumes the following form [16]: As h is the gap between the upper "s,1" and the lower "s,2" solid plates of the bearing, then the length scale in the Reynolds' number is half the fluid depth. According to [38], the 1D infinitive linear slider shows a laminar regime for Re≤290. As the current validation refers to a finite slider, a lower 23 value (Re=100) is used here to guarantee a laminar regime along the 95% of the slider length (far enough from the leading and the trailing edges of the mobile plate).
The initial fluid velocity is the average velocity of the solid plates (uf,0=us,1/2). The reference velocity is U * =50m/s. A high viscosity liquid allows representing a motor oil and makes the presence of the artificial viscosity (Sec.4.1) irrelevant.
Omitting gravity is equivalent and alternative to replacing pressure with the reduced pressure, which is defined as the difference between pressure and its hydrostatic component. Under the The spatial resolution is defined by dx=2.1×10 -6 m, hSPH/dx=1.3 and dx/dxs=2. Since the spatial discretization of the numerical pressure profile is uniform, the numerical non-dimensional load bearing capacity is estimated as the average of the pressure coefficient values over the mobile plate bottom:  Although the configuration of this uniform slider is 2D, the model is able to represent an equivalent 3D configuration, with a proper representation of the fluid volume, the solid volume, the fixed bottom and the symmetry planes. This feature is relevant for its applicability to complex topological configurations, see Sec. 8. equations for fluid films are not fully respected, as also observed by Vakilian et al. [3] and Dobrica and Fillon [2]. Contrarily to the analytical solution, the code can simulate the z-dependent pressure field, the velocity vertical component and the inertial effects due to the fluid-structure interactions at the leading and the trailing edges. Further, the inter-comparisons between the SPH profiles at different heights show the accuracy of the code in reproducing the uniformity of the pressure field along z-axis, far enough from the plate edges (Fig. 3, centre panel).

Applicability of the computational framework to complex rough surfaces
The linear slider bearing of Sec. 7 is herein modified by imposing a complex 3D surface as the bottom fixed frontier, to demonstrate the versatility of the approach to take into account complex rough surfaces as input data. Its shape is taken from a real natural geometry of a strawberry leaf that was acquired using the Leica DCM 3D confocal profilometer in the MUSAM-Lab of the IMT School for Advanced Studies Lucca. The spatial resolution needs be finer than the previous test cases due to roughness, with dx=1.66×10 -6 m. This study allowed SPHERA to deal with hydrodynamic lubrication and empowered the code on other application fields with fluid-structure interactions (e.g., transport of solid bodies by floods and earth landslides; rock landslides). SPHERA is developed and distributed on a GitHub public repository.