An SPH framework for fluid–solid and contact interaction problems including thermo-mechanical coupling and reversible phase transitions

Abstract

The present work proposes an approach for fluid–solid and contact interaction problems including thermo-mechanical coupling and reversible phase transitions. The solid field is assumed to consist of several arbitrarily-shaped, undeformable but mobile rigid bodies, that are evolved in time individually and allowed to get into mechanical contact with each other. The fluid field generally consists of multiple liquid or gas phases. All fields are spatially discretized using the method of smoothed particle hydrodynamics (SPH). This approach is especially suitable in the context of continually changing interface topologies and dynamic phase transitions without the need for additional methodological and computational effort for interface tracking as compared to mesh- or grid-based methods. Proposing a concept for the parallelization of the computational framework, in particular concerning a computationally efficient evaluation of rigid body motion, is an essential part of this work. Finally, the accuracy and robustness of the proposed framework is demonstrated by several numerical examples in two and three dimensions, involving multiple rigid bodies, two-phase flow, and reversible phase transitions, with a focus on two potential application scenarios in the fields of engineering and biomechanics: powder bed fusion additive manufacturing (PBFAM) and disintegration of food boluses in the human stomach. The efficiency of the parallel computational framework is demonstrated by a strong scaling analysis.

Introduction

In many applications in science and engineering, like for example in some areas of biomechanics, fluid–solid and contact interaction problems characterized by a large number of solid bodies immersed in a fluid flow and undergoing reversible phase transitions, are of great interest. Often, explicitly considering the deformation of solid bodies can be neglected, which reduces the complexity of the problem to the treatment of undeformable but mobile rigid bodies, in favor of simplified modeling. Most current mesh- or grid-based methods, e.g., the finite element method (FEM), the finite difference method (FDM), or the finite volume method (FVM), require substantial methodological and computational efforts to model the motion of rigid bodies in fluid flow. To overcome those issues, several approaches, e.g., based on the particle finite element method (PFEM) [1,2,3], on the discrete element method (DEM) [4,5,6,7], or on smoothed particle hydrodynamics (SPH) [8,9,10,11,12,13,14], have been proposed. SPH as a mesh-free discretization scheme is, due to its Lagrangian nature, very well suited for flow problems involving multiple phases, dynamic and reversible phase transitions, and complex interface topologies. This makes SPH very appropriate for a wide range of applications in engineering, e.g., in metal additive manufacturing melt pool modeling [15, 16], or in biomechanics, e.g., for modeling the digestion of food in the human stomach [17]. For the former application, an SPH formulation for thermo-capillary phase transition problems involving solid, liquid, and gaseous phases has recently been proposed [18], amongst others, focusing on evaporation-induced recoil pressure forces, temperature-dependent surface tension and wetting forces, Gaussian laser beam heat sources, and evaporation-induced heat losses. However, for simplicity, this and also other current state-of-the-art approaches, e.g., [19, 20] are restricted to immobile powder grains.

All aforementioned SPH formulations for modeling rigid body motion in fluid flow have in common, that rigid bodies are fully resolved, that is spatially discretized as clusters of particles. It is generally accepted that advanced boundary particle methods, e.g., based on the extrapolation of field quantities from fluid to boundary particles [21,22,23], are beneficial, because one can model the fluid field close to the boundary with high accuracy. In many of the aforementioned applications, an exact representation of the fluid–solid interface plays an important role. Therefore, herein a formulation of this kind proposed in [23] is utilized. To the best of the authors’ knowledge none of the aforementioned SPH formulations modeling rigid body motion in fluid flow simultaneously consider thermal conduction, reversible phase transitions, and multiple (liquid and gas) phases.

To help close this gap, this contribution proposes a fully resolved smoothed particle hydrodynamics framework for fluid–solid and contact interaction problems including thermo-mechanical coupling and reversible phase transitions. The solid field is assumed to consist of several arbitrarily-shaped, undeformable but mobile rigid bodies, that are evolved in time individually. Based on a temperature field, provided by solving the heat equation, reversible phase transitions, i.e., melting and solidification, are evaluated between the fluid and the solid field. As a result, the shape and the total number of rigid bodies may vary over time. In addition, contact between rigid bodies is considered by employing a spring-dashpot model. Note that some characteristic phenomena for thermo-capillary flow [24, 25] and, especially, for metal PBFAM melt pool modeling [18,19,20, 26,27,28] are not addressed in this work, thus, referring to the literature.

While parallel implementation aspects along with detailed scalability studies are not in the focus of the aforementioned references, in this work, a concept for the parallelization of the computational framework is proposed, setting the focus in particular on an efficient evaluation of rigid body motion. The parallel behavior is demonstrated, confirming that detailed studies at a large scale become possible. It shall be noted, that the parallel implementation of such a computational framework is far from trivial but indispensable when examining numerical examples that are of practical relevance. Note that the introduced concept for the parallelization of the computational framework is applicable not only when using SPH as a discretization scheme, but also for other particle-based methods, e.g., discrete element method (DEM), or molecular dynamics (MD).

The remainder of this work is organized as follows: To begin with, the governing equations for a fluid–solid and contact interaction problem including thermo-mechanical coupling and reversible phase transition are outlined. Next, the numerical methods are presented and the details of the computational framework are discussed. Finally, the accuracy and robustness of the proposed formulation is demonstrated by several numerical examples.

Governing equations

Consider a domain $\Omega$ of a fluid–solid interaction problem that consists at each time $t \in [0,T]$ of the non-overlapping fluid domain $\Omega ^{f}$ and solid domain $\Omega ^{s}$ that share a common interface $\Gamma ^{fs}$, with $\Omega = \Omega ^{f} \cup \Omega ^{s}$ and $\Omega ^{f} \cap \Omega ^{s} = \Gamma ^{fs}$. In general, the fluid domain $\Omega ^{f}$ may consist of multiple (liquid and gas) phases. For ease of notation, in the following it will not be distinguished between the different fluid phases. The solid domain $\Omega ^{s}$ is composed of several non-overlapping sub-domains $\Omega ^{s}_{k}$, which represent rigid bodies k, such that $\Omega ^{s} = \bigcup _{k} \Omega ^{s}_{k}$. In the event of contact between two rigid bodies k and $\hat{k}$, a common interface $\Gamma ^{ss}_{k,\hat{k}} = \Omega ^{s}_{k} \cap \Omega ^{s}_{\hat{k}}$ exists, separating the two solid sub-domains $\Omega ^{s}_{k}$ and $\Omega ^{s}_{\hat{k}}$. A detailed illustration of the problem is given in Fig. 1. In the following the (standard) governing equations of the fluid and the solid field as well as the respective coupling conditions are briefly given. In addition, reversible phase transitions between the fluid and the solid field, e.g., temperature-induced melting and solidification, may occur. For this reason, the temperature field is modeled solving the heat equation.

Fig. 1

Domain $\Omega$ consisting of several disjunct domains, the fluid domain $\Omega ^{f}$ and the solid sub-domains $\Omega ^{s}_{k}$ representing rigid bodies k, with fluid–solid interface $\Gamma ^{fs}_{k}$ and solid-solid interface $\Gamma ^{ss}_{k,\hat{k}}$ in the event of contact between the rigid bodies k and $\hat{k}$

Fluid field

The fluid field is governed by the instationary Navier–Stokes equations in the domain $\Omega ^{f}$, which consist in convective form of the mass continuity equation and the momentum equation

$$\begin{aligned} \frac{\mathrm{d}\rho ^{f}}{\mathrm{d}t}= & {} -\rho ^{f} {\varvec{\nabla }}\cdot \mathbf {u}^{f} {\quad \text {in}\quad }\Omega ^{f} \, , \end{aligned}$$

(1)

$$\begin{aligned} \frac{\mathrm{d}\mathbf {u}^{f}}{\mathrm{d}t}= & {} -\frac{1}{\rho ^{f}} {\varvec{\nabla }}_{p^{f}} + \mathbf {f}_{\nu } + \mathbf {b}^{f} {\quad \text {in}\quad }\Omega ^{f} \, , \end{aligned}$$

(2)

with viscous force $\mathbf {f}_{\nu }$ and body force $\mathbf {b}^{f}$ each per unit mass. For a Newtonian fluid the viscous force is $\mathbf {f}_{\nu } = \nu ^{f} \nabla ^{2}{\mathbf {u}^{f}}$ with kinematic viscosity $\nu ^{f}$. The mass continuity equation (1) and the momentum equation (2) represent a system of $d+1$ equations with the $d+2$ unknowns, velocity $\mathbf {u}^{f}$, density $\rho ^{f}$, and pressure $p^{f}$, in d-dimensional space. The system of equations is closed with an equation of state $p^{f} = p^{f}(\rho ^{f})$ relating fluid density $\rho ^{f}$ and pressure $p^{f}$. The Navier–Stokes equations (1) and (2) are subject to the following initial conditions

$$\begin{aligned} \rho ^{f} = \rho ^{f}_{0} {{\quad }\text {and}\quad }\mathbf {u}^{f} = \mathbf {u}^{f}_{0} {\quad \text {in}\quad }\Omega ^{f} \quad \text {at}\quad t = 0 \end{aligned}$$

(3)

with initial density $\rho ^{f}_{0}$ and initial velocity $\mathbf {u}^{f}_{0}$. In addition, Dirichlet and Neumann boundary conditions are applied on the fluid boundary $\Gamma ^{f} = \partial \Omega ^{f} \setminus \Gamma ^{fs}$

$$\begin{aligned} \mathbf {u}^{f} = \hat{\mathbf {u}}^{f} \quad \text {on}\quad \Gamma ^{f}_{D} {{\quad }\text {and}\quad }\mathbf {t}^{f} = \hat{\mathbf {t}}^{f} \quad \text {on}\quad \Gamma ^{f}_{N} \, , \end{aligned}$$

(4)

with prescribed boundary velocity $\hat{\mathbf {u}}^{f}$ and boundary traction $\hat{\mathbf {t}}^{f}$, where $\Gamma ^{f} = \Gamma ^{f}_{D} \cup \Gamma ^{f}_{N}$ and $\Gamma ^{f}_{D} \cap \Gamma ^{f}_{N} = \emptyset$. Furthermore, on the fluid–solid interface $\Gamma ^{fs} = \bigcup _{k} \Gamma ^{fs}_{k}$ the so-called kinematic and dynamic coupling conditions are

$$\begin{aligned} \mathbf {u}^{f} = \mathbf {u}^{fs}_{k} {{\quad }\text {and}\quad }\mathbf {t}^{f} = \mathbf {t}^{fs}_{k} \quad \text {on}\quad \Gamma ^{fs}_{k} \quad \forall {k} \, , \end{aligned}$$

(5)

resembling a no-slip boundary condition and ensuring equilibrium of fluid and solid traction across the interface $\Gamma ^{fs}$. Herein, $\mathbf {u}^{fs}_{k}$ and $\mathbf {t}^{fs}_{k}$ denote the velocity respectively traction of a rigid body k on the fluid–solid interface $\Gamma ^{fs}_{k}$.

Remark 1

In Eqs. (1)–(2) governing the fluid field, all time derivatives follow the motion of material points, i.e., are material derivatives $\frac{\mathrm{d}(\cdot )}{\mathrm{d}t} = \frac{{\partial }(\cdot )}{{\partial }t} + \mathbf {u} \cdot {\varvec{\nabla }}{(\cdot )}$. Besides, ${\varvec{\nabla }}{(\cdot )}$ denotes derivatives with respect to spatial coordinates.

Solid field

The solid field is assumed to consist of several mobile rigid bodies k each represented by a sub-domain $\Omega ^{s}_{k}$ embedded in the fluid domain $\Omega ^{f}$. Thus, the interface of a rigid body k is $\Gamma ^{s}_{k} = \Gamma ^{fs}_{k} \cup \left( \bigcup _{\hat{k}} \Gamma ^{ss}_{k, \hat{k}} \right)$ with contacting rigid bodies $\hat{k}$, cf. Fig. 1. The kinematics of each rigid body k are uniquely defined by three respectively six degrees of freedom in two- and three-dimensional space, i.e., the position of the center of mass $\mathbf {r}^{s}_{k}$ and the orientation $\varvec{\psi }^{s}_{k}$. As a result, the equations of motion of an individual rigid body k are described by the balance of linear and angular momentum

$$\begin{aligned} m^{s}_{k} \frac{\mathrm{d}^{2}{\mathbf {r}^{s}_{k}}}{\mathrm{d}{t}^{2}}= & {} \mathbf {f}^{fs}_{k} + \sum _{\hat{k}} \mathbf {f}^{ss}_{k,\hat{k}} + m^{s}_{k} \mathbf {b}^{s}_{k} {\quad \text {in}\quad }\Omega ^{s}_{k} \, , \end{aligned}$$

(6)

$$\begin{aligned} \mathbf {I}^{s}_{k} \frac{\mathrm{d}\varvec{\omega }^{s}_{k}}{\mathrm{d}t}= & {} \mathbf {m}^{fs}_{k} + \sum _{\hat{k}} \mathbf {m}^{ss}_{k,\hat{k}} {\quad \text {in}\quad }\Omega ^{s}_{k} \, , \end{aligned}$$

(7)

with mass $m^{s}_{k}$ and mass moment of inertia $\mathbf {I}^{s}_{k}$ with respect to the center of mass position $\mathbf {r}^{s}_{k}$. Herein, $\varvec{\omega }^{s}_{k}$ denotes the angular velocity of a rigid body k, cf. Remark 2. Furthermore, $\mathbf {f}^{fs}_{k}$ and $\mathbf {m}^{fs}_{k}$ describe the resultant coupling force respectively torque acting on the fluid–solid interface $\Gamma ^{fs}_{k}$ of rigid body k. Contacting rigid bodies k and $\hat{k}$ exchange the resultant contact force $\mathbf {f}^{ss}_{k,\hat{k}}$ respectively torque $\mathbf {m}^{ss}_{k,\hat{k}}$ at the solid-solid interface $\Gamma ^{ss}_{k, \hat{k}}$. Finally, the body force $\mathbf {b}^{s}_{k}$ given per unit mass is contributing to the balance of linear momentum.

Remark 2

The orientation $\varvec{\psi }^{s}_{k}$ is expressed by a (pseudo-)vector whose direction and magnitude represent the axis and angle of rotation. Note that in general, the angular velocity $\varvec{\omega }^{s}_{k}$ of a rigid body k is different from the time derivative of the orientation $\varvec{\psi }^{s}_{k}$, i.e., $\varvec{\omega }^{s}_{k} \ne \mathrm{d}{\varvec{\psi }^{s}_{k}}/\mathrm{d}{t}$, due to the non-additivity of large rotations [29,30,31,32]. Direct evolution of the orientation $\varvec{\psi }^{s}_{k}$ of a rigid body k requires a special class of time integration schemes, so-called Lie group time integrators [33, 34].

Thermal conduction

Thermal conduction in the combined fluid and solid domain $\Omega = \Omega ^{f} \cup \Omega ^{s}$ in the absence of heat sources or heat sinks (which are neglected herein for simplicity) is governed by the heat equation

$$\begin{aligned} c_{p}^{\phi } \frac{\mathrm{d}T}{\mathrm{d}t} = \frac{1}{\rho ^{\phi }} {\varvec{\nabla }} ( \kappa ^{\phi } {\varvec{\nabla }}{T} ) {\quad \text {in}\quad }\Omega \, , \end{aligned}$$

(8)

with temperature T and heat flux $\mathbf {q}=-\kappa ^{\phi } {\varvec{\nabla }}{T}$. The material parameters heat capacity $c_{p}^{\phi }$ and thermal conductivity $\kappa ^{\phi }$ are in general different for fluid and solid field, and hence for clarity are denoted by the index $(\cdot )^{\phi }$ with $\phi \in \{f,s\}$. The heat equation (8) is subject to the following initial condition

$$\begin{aligned} T = T_{0} {\quad \text {in}\quad }\Omega \quad \text {at}\quad t = 0 \end{aligned}$$

(9)

with initial temperature $T_{0}$. In addition, Dirichlet and Neumann boundary conditions are required on the domain boundary $\Gamma = \partial \Omega$

$$\begin{aligned} T = \hat{T} \quad \text {on}\quad \Gamma _{D} {{\quad }\text {and}\quad }\mathbf {q} = \hat{\mathbf {q}} \quad \text {on}\quad \Gamma _{N} \, , \end{aligned}$$

(10)

with prescribed boundary temperature $\hat{T}$ and boundary heat flux $\hat{\mathbf {q}}$, where $\Gamma = \Gamma _{D} \cup \Gamma _{N}$ and $\Gamma _{D} \cap \Gamma _{N} = \emptyset$.

Reversible phase transition

Reversible phase transitions, i.e., melting and solidification, are considered between the solid and the fluid field. Within this publication, solid material points that exceed a transition temperature $T_{t}$ undergo a phase transition to a fluid material point and vice versa, cf. Remark 3. Consequently, the shape of a rigid body k, i.e., its sub-domain $\Omega ^{s}_{k}$, is changing due to a loss or gain of material points resulting in a varying mass $m_{k}$, center of mass position $\mathbf {r}_{k}$, and mass moment of inertia $\mathbf {I}_{k}$.

Remark 3

For the sake of simplicity, only temperature-independent parameters are considered herein, and latent heat is neglected. Latent heat could be included by an apparent capacity scheme relying on an increased heat capacity $c_{p}$ within a finite temperature interval [35] in a straightforward manner.

Remark 4

The proposed framework is general enough to model also chemically-induced phase transitions based on a concentration field. For this purpose, the diffusion equation $\mathrm{d}{C}/\mathrm{d}{t} = {1}/{\rho ^{\phi }} {\varvec{\nabla }} ( D^{\phi } {\varvec{\nabla }}{C} )$ with diffusivity $D^{\phi }$ modeling the transport of a concentration C is solved. Considering the similarity between the heat equation (8) and the diffusion equation, the latter can be discretized following a similar SPH discretization [36, 37] as applied for the heat equation. Similarly, modeling phase transitions a transition concentration $C_{t}$ is defined.

Numerical methods and parallel computational framework

This section presents the methods applied for discretization and numerical solution of fluid–solid and contact interaction problems including thermo-mechanical coupling and reversible phase transitions. The presented parallel computational framework is implemented in the in-house parallel multiphysics research code BACI (Bavarian Advanced Computational Initiative) [38] using the Message Passing Interface (MPI) for distributed-memory parallel programming.

Spatial discretization via smoothed particle hydrodynamics

For the spatial discretization smoothed particle hydrodynamics (SPH) is used, allowing for a straightforward particle-based evaluation of fluid–solid coupling conditions. In the following, the basics of this method are recapitulated briefly.

Approximation of field quantities applying a smoothing kernel

The fundamental concept of SPH is based on the approximation of a field quantity f via a smoothing operation and on the discretization of the domain $\Omega$ with discretization points, so-called particles j, each occupying a volume $V_{j}$. Introducing a smoothing kernel $W(r,h)$ that fulfills certain consistency properties [36, 39], cf. Remark 5, leads to an approximation of the field quantity f based on summation of contributions from all particles j in the domain $\Omega$

$$\begin{aligned} f(\mathbf {r}) \approx \int _{\Omega } f(\mathbf {r}') W(| \mathbf {r} - \mathbf {r}' |, h) \mathrm{d}{\mathbf {r}'} \approx \sum _{j} V_{j} f(\mathbf {r}_{j}) W(| \mathbf {r} - \mathbf {r}_{j} |, h) \end{aligned}$$

(11)

which includes a smoothing error and a discretization error [40].

Remark 5

The smoothing kernel $W(r,h)$ is a monotonically decreasing, smooth function that depends on a distance r and a smoothing length h. The smoothing length h together with a scaling factor $\kappa$ define the support radius $r_{c} = \kappa h$ of the smoothing kernel. Compact support, i.e., $W(r, h) = 0$ for $r > r_{c}$, as well as positivity, i.e., $W(r, h) \ge 0$ for $r \le r_{c}$, are typical properties of standard smoothing kernels $W(r,h)$. In addition, the normalization property requires that $\int _{\Omega } W(| \mathbf {r} - \mathbf {r}' |, h) \mathrm{d}{\mathbf {r}'} = 1$. The Dirac delta function property $\lim _{h \rightarrow 0}{ W(r, h) } = \delta (r)$ ensures an exact representation of a field quantity f in the limit $h \rightarrow 0$.

A straightforward approach in SPH to determine the gradient of a quantity f follows directly by differentiation of (11) resulting in

$$\begin{aligned} {\varvec{\nabla }}{f(\mathbf {r})} \approx \int _{\Omega } f(\mathbf {r}') {\varvec{\nabla }}{W(| \mathbf {r} - \mathbf {r}_{j} |, h)}\,\mathrm{d}{\mathbf {r}'} \approx \sum _{j} V_{j} f(\mathbf {r}_{j}) {\varvec{\nabla }}{W(| \mathbf {r} - \mathbf {r}_{j} |, h)} \, . \end{aligned}$$

(12)

Note that this (simple) gradient approximation shows some particular disadvantages. Hence, more advanced approximations for gradients are given in the literature [36] and will also be applied in the following. In sum, the concept of SPH allows to reduce partial differential equations to a system of coupled ordinary differential equations (with as many equations as particles) that is solved in the domain $\Omega$. Thereby, all field quantities are evaluated at and associated with particle positions, meaning each particle carries its corresponding field quantities. Finally, in a post-processing step a continuous field quantity f is recovered from the discrete quantities $f(\mathbf {r}_{j})$ of particles j in the domain $\Omega$ using the approximation (11) and the commonly known Shepard filter

$$\begin{aligned} \hat{f}(\mathbf {r}) \approx \frac{ \sum _{j} V_{j} f(\mathbf {r}_{j}) W(| \mathbf {r} - \mathbf {r}_{j} |, h) }{ \sum _{j} V_{j} W(| \mathbf {r} - \mathbf {r}_{j} |, h) } \, . \end{aligned}$$

(13)

Note that the denominator typically takes on values close to one inside the fluid domain and is mainly relevant for boundary regions with reduced support due to a lack of neighboring particles.

Remark 6

In the following, a quantity f evaluated for particle i at position $\mathbf {r}_{i}$ is written as $f_{i} = f(\mathbf {r}_{i})$. The short notation $W_{ij} = W(r_{ij}, h)$ denotes the smoothing kernel W evaluated for particle i at position $\mathbf {r}_{i}$ with neighboring particle j at position $\mathbf {r}_{j}$, where $r_{ij} = |\mathbf {r}_{ij}| = |\mathbf {r}_{i} - \mathbf {r}_{j}|$ is the absolute distance between particles i and j. The derivative of the smoothing kernel W with respect to the absolute distance $r_{ij}$ is denoted by $\partial {W}/\partial {r_{ij}} = \partial {W(r_{ij}, h)}/\partial {r_{ij}}$.

Remark 7

Herein, a quintic spline smoothing kernel $W(r, h)$ as defined in [21] with smoothing length h and compact support of the smoothing kernel with support radius $r_{c} = \kappa h$ and scaling factor $\kappa = 3$ is used.

Initial particle spacing

Within this contribution, the domain $\Omega$ is initially filled with particles located on a regular grid with particle spacing $\Delta {}x$, thus in the d-dimensional space each particle initially occupies an effective volume $V_{eff}=(\Delta {}x)^{d}$. A particle in the fluid domain $\Omega ^{f}$ is called a fluid particle i, whereas a particle in the solid domain $\Omega ^{s}_{k}$ of a rigid body k is called a rigid particle r. It follows, that each rigid body is fully resolved being spatially discretized as clusters of particles. Naturally, the choice of the particle spacing $\Delta {}x$ influences the accuracy of the interface representation between fluid and solid domain. The mass of a particle is initially assigned using the reference density of the respective phase, i.e., $\rho ^{f}_{0}$ for the fluid phase and $\rho ^{s}_{0}$ for the solid phase, and the effective volume $V_{eff}$.

Remark 8

Within this work, the smoothing length h of the smoothing kernel $W(r, h)$, cf. Remark 7, is set equal to the initial particle spacing $\Delta {}x$. Consequently, in a convergence analysis with decreasing particle spacing $\Delta {}x$ the ratio ${\Delta {}x}/{h}$ remains constant [40].

Parallelization via spatial decomposition of the domain

For the problems studied herein, an efficient parallel computational framework capable of handling systems constituted of a large number of particles is required. This requires addressing in particular two aspects, namely, an efficient particle neighbor pair detection, and a parallel load distribution strategy while keeping the communication overhead at an acceptable level. In the literature, several approaches for parallel computational frameworks utilizing particle-based methods have been proposed, e.g., [41,42,43,44,45,46,47]. In the present work, a spatial decomposition approach with neighbor pair detection utilizing a combination of cell-linked lists and Verlet-lists based on [42] is applied. The general idea of the spatial decomposition approach is briefly explained in the following, however, for detailed information, the interested reader is referred to the original publication [42]. In addition, the concept is extended herein to consider the motion of rigid bodies spatially discretized as clusters of particles.

Fig. 2

A fluid–solid interaction problem consisting of a rigid body k with affiliated rigid particles r and surrounding fluid particles i distributed over several processors p according to a spatial decomposition approach

The evaluation of particle interactions in SPH requires knowledge of neighboring particles within a geometrically limited interaction distance, i.e., within the support radius $r_{c}$ of the smoothing kernel. Thus, the computational domain is divided into several cubic cells forming a uniform lattice, while each particle is uniquely assigned to one of those cells according to its current spatial position, cf. Fig. 2. The size of the cells is chosen such that neighboring particles are either located in the same cell or in adjacent cells, i.e., the size of the cells is at least equal to the support radius $r_{c}$ of the smoothing kernel.

Following a spatial decomposition approach, the cells together with assigned particles are distributed over all involved processors, i.e., forming so-called processor domains. To keep the computational load balanced between all processors and to minimize the communication overhead, cubic processor domains are defined such that each contains (nearly) the same number of particles. The cells occupied by each processor are called owned cells. On each processor the position of particles located in its processor domain, i.e., the position of so-called owned particles, is evolved. This requires the evaluation of interactions of owned particles with their neighboring particles. However, the correct evaluation of particle interactions close to processor domain boundaries requires that each processor has information not only about its owned particles but also about particles in cells adjacent to its processor domain. To this end, each processor is provided full information not only about its own domain but additionally about a layer of ghosted cells (with ghosted particles) around its own domain. Keeping the information about ghosted cells and particles continuously updated requires communication between processors.

Remark 9

To exemplify the cost of communication overhead, consider a perfectly cubic processor domain occupying $n_{o}$ owned cells. Consequently, assuming one layer of ghosted cells surrounding the processor domain, a total of $n_{g} = ( \root 3 \of {n_{o}} + 2 )^{3} - n_{o}$ cells are ghosted. That is, the communication overhead scales with the ratio ${n_{g}}/{n_{o}}$ of ghosted cells $n_{g}$ to owned cells $n_{o}$. Furthermore, the (average) number of particles per cell, and, consequently, also the communication overhead, scale with the ratio ${r_{c}}/{\Delta {}x}$ of the support radius $r_{c}$ and the initial particle spacing $\Delta {}x$.

As a consequence of the spatial decomposition approach, the affiliated rigid particles r of a rigid body k might be distributed over several processors, cf. Fig. 2. However, note that the balance of linear and angular momentum, cf. Eqs. (6)–(7), describing the motion of a rigid body k are given with respect to the center of mass position $\mathbf {r}^{s}_{k}$. Thus, the evaluation of mass quantities, i.e., mass $m_k$, center of mass position $\mathbf {r}^{s}_{k}$, and mass moment of inertia $\mathbf {I}_{k}$, as well as the evaluation of resultant force $\mathbf {f}_{k}$ and torque $\mathbf {m}_{k}$ acting on a rigid body k, requires special communication between all processors hosting rigid particles r belonging to rigid body k and the single processor owning rigid body k.

Modeling fluid flow using weakly compressible SPH

For modeling fluid flow using SPH, several different formulations each with its own characteristics and benefits can be derived. Here, the instationary Navier–Stokes equations (1) and (2) are discretized by a weakly compressible approach [36, 39, 48]. This section gives a brief overview of this formulation applied already in [18, 49]. For ease of notation, in the following the index $(\cdot )^{f}$ denoting fluid quantities is dropped.

Density summation

The density of a particle i is determined via summation of the respective smoothing kernel contributions of all neighboring particles j within the support radius $r_{c}$

$$\begin{aligned} \rho _{i} = m_{i} \sum _{j} W_{ij} \end{aligned}$$

(14)

with mass $m_{i}$ of particle i. This approach is typically denoted as density summation and results in an exact conservation of mass in the fluid domain, which can be shown in a straightforward manner considering the commonly applied normalization of the smoothing kernel to unity. It shall be noted that the density field may alternatively be obtained by discretization and integration of the mass continuity equation (1) [39].

Momentum equation

The momentum equation (2) is discretized following [23, 50] including a transport velocity formulation to suppress the problem of tensile instability. It will be briefly recapitulated in the following. The transport velocity formulation relies on a constant background pressure $p_{b}$ that is applied to all particles and results in a contribution to the particle accelerations for in general disordered particle distributions. However, these additional acceleration contributions vanish for particle distributions fulfilling the partition of unity of the smoothing kernel, thus fostering these desirable configurations. For the sake of brevity, the definition of the modified advection velocity and the additional terms in the momentum equation from the aforementioned transport velocity formulation are not discussed in the following and the reader is referred to the original publication [50]. Altogether, the acceleration $\mathbf {a}_{i} = \mathrm{d}{\mathbf {u}_{i}}/\mathrm{d}{t}$ of a particle i results from summation of all acceleration contributions due to interaction with neighboring particles j and a body force as

$$\begin{aligned} \mathbf {a}_{i} = \frac{1}{m_{i}} \sum _{j} (V_{i}^{2}+V_{j}^{2}) \left[ - \tilde{p}_{ij} \frac{{\partial }W}{{\partial }r_{ij}} \mathbf {e}_{ij} + \tilde{\eta }_{ij} \frac{\mathbf {u}_{ij}}{r_{ij}} \frac{{\partial }W}{{\partial }r_{ij}} \right] + \mathbf {b}_{i} \, , \end{aligned}$$

(15)

with volume $V_{i} = m_{i}/\rho _{i}$ of particle i, unit vector $\mathbf {e}_{ij} = {\mathbf {r}_{i} - \mathbf {r}_{j}}/{|\mathbf {r}_{i} - \mathbf {r}_{j}|} = {\mathbf {r}_{ij}}/{r_{ij}}$, relative velocity $\mathbf {u}_{ij} = \mathbf {u}_{i}-\mathbf {u}_{j}$, and density-weighted inter-particle averaged pressure and inter-particle averaged viscosity

$$\begin{aligned} \tilde{p}_{ij} = \frac{\rho _{j}p_{i}+\rho _{i}p_{j}}{\rho _{i} + \rho _{j}} {{\quad }\text {and}\quad }\tilde{\eta }_{ij} = \frac{2\eta _{i}\eta _{j}}{\eta _{i}+\eta _{j}} \, . \end{aligned}$$

(16)

In the following, the acceleration contribution of a neighboring particle j to particle i is, for ease of notation, denoted as $\mathbf {a}_{ij}$, where $\mathbf {a}_{i} = \sum _{j} \mathbf {a}_{ij} + \mathbf {b}_{i}$. The above given momentum equation (15) exactly conserves linear momentum due to pairwise anti-symmetric particle forces

$$\begin{aligned} m_{i} \mathbf {a}_{ij} = - m_{j} \mathbf {a}_{ji} \, , \end{aligned}$$

(17)

which follows from the property $\partial {W}/\partial {r_{ij}} = \partial {W}/\partial {r_{ji}}$ of the smoothing kernel.

Equation of state

Following a weakly compressible approach, the density $\rho _{i}$ and pressure $p_{i}$ of a particle i are linked via the equation of state

$$\begin{aligned} p_{i}(\rho _{i}) = c^{2} (\rho _{i} - \rho _{0}) = p_{0} \left( \frac{\rho _{i}}{\rho _{0}} - 1\right) \end{aligned}$$

(18)

with reference density $\rho _{0}$, reference pressure $p_{0} = \rho _{0} c^{2}$ and artificial speed of sound c. Note that this commonly applied approach can only capture deviations from the reference pressure, i.e., $p_{i}(\rho _{0}) = 0$, and not the total pressure. To limit density fluctuations to an acceptable level, while still avoiding too severe time step restrictions, cf. Eq. (41), strategies are discussed in [21] how to choose a reasonable value for the artificial speed of sound c. Accordingly, in this work the artificial speed of sound c is set allowing an average density variation of approximately $1\%$.

Boundary and coupling conditions

Herein, both rigid wall boundary conditions as well as rigid body coupling conditions, are modeled following [23]. In the former case, at least $q = \mathrm {floor}({r_{c}}/{\Delta {}x})$ layers of boundary particles b are placed parallel to the fluid boundary $\Gamma ^{f}_{D}$ with a distance of ${\Delta {}x}/{2}$ outside of the fluid domain $\Omega ^{f}$ in order to maintain full support of the smoothing kernel. In the latter case, rigid particles r of rigid bodies k are considered, while naturally describing the fluid–solid interface $\Gamma ^{fs}$. In both cases, a boundary particle b or a rigid particle r contribute to the density summation (14) and to the momentum equation (15) evaluated for a fluid particle i considered as neighboring particle j. The respective quantities of boundary particles b respectively rigid particles r are extrapolated from the fluid field based on a local force balance as described in [23]. Consequently, striving for conservation of linear momentum, cf. Eq. (17), the force acting on a rigid particle r stemming from interaction with fluid particle i, cf. Eq. (15), is given as

$$\begin{aligned} \mathbf {f}_{ri} = - m_{i} \mathbf {a}_{ir} \, . \end{aligned}$$

(19)

Remark 10

The $\mathrm {floor}$ operator is defined by $\mathrm {floor}(x) := \max \, \{ k \in \mathbb {Z} \mid k \le x\}$ and returns the largest integer that is less than or equal to its argument x.

Modeling the motion of rigid bodies discretized by particles

Within this formulation, each rigid body k is composed of several rigid particles r that are fixed relative to a rigid body frame, i.e., there is no relative motion among rigid particles of a rigid body. Thus, the rigid particles of a rigid body are not evolved in time individually, but follow the motion of the rigid body described by the balance of linear and angular momentum, cf. Eqs. (6)–(7). As a consequence of the spatial decomposition approach, special communication between all processors hosting rigid particles r of a rigid body k in terms of the evaluation of mass quantities, or resultant forces and torques, is required. For ease of notation, in the following the index $(\cdot )^{s}$ denoting solid quantities is dropped.

Orientation of rigid bodies

The orientation $\varvec{\psi }_{k}$ of a rigid body k, described by one respectively three degrees of freedom in two- and three-dimensional space, is previously introduced without explicitly defining a specific parameterization of the underlying rotation, e.g., via Euler angles or Rodriguez parameters. Moreover, as stated in Remark 2, explicit evolution of the orientation $\varvec{\psi }_{k}$ requires special Lie group time integrators. A straightforward approach to overcome aforementioned issues is to describe the orientation of a rigid body via quaternion algebra, cf. Remark 12. Consequently, in the following it is assumed that at all times the orientation $\varvec{\psi }_{k}$ of a rigid body k can be uniquely described by a unit quaternion $\mathbf {q}_{k}$, cf. Fig. 3. Once a local rigid body frame is defined, this allows to transform the relative position of rigid particles $\mathbf {r}_{rk}$, cf. Remark 11, from that rigid body frame to the reference frame, e.g., a global cartesian system.

Fig. 3

Orientation of a rigid body k with rigid particles r and their relative position $\mathbf {r}_{rk}$ described via a unit quaternion $\mathbf {q}_{k}$ at different times

Remark 11

Note that the relative position of rigid particles $\mathbf {r}_{rk}$ expressed in the rigid body frame is, in general, a known (and constant) quantity, that only needs to be updated in case the center of mass position $\mathbf {r}_{k}$ changes, i.e., due to phase transitions.

Remark 12

For the sake of brevity, the principals of quaternion algebra are not delineated herein. It remains the definition of operator $\circ$ denoting quaternion multiplication as used in the following.

Parallel evaluation of mass-related quantities

In a first step, on each processor p the processor-wise mass ${}_{p}{m}_{k}$ and center of mass position ${}_{p}{\mathbf {r}}_{k}$ of a rigid body k are computed as

$$\begin{aligned} {}_{p}{m}_{k} = \sum _{r} m_{r} {{\quad }\text {and}\quad }{}_{p}{\mathbf {r}}_{k} = \frac{\sum _{r} m_{r} \mathbf {r}_{r}}{\sum _{r} m_{r}} \end{aligned}$$

(20)

considering the mass $m_{r}$ and position $\mathbf {r}_{r}$ of all affiliated rigid particles r being located in the computational domain of processor p, cf. Fig. 4 for an illustration. Accordingly, the processor-wise mass moment of inertia ${}_{p}{\mathbf {I}}_{k}$ of a rigid body k follows componentwise (in index notation) as

$$\begin{aligned} {}_{p}{I}_{k,ij} = \sum _{r} \left[ I_{r} \, \delta _{ij} + \left[ \sum _{q} ( {}_{p}{r}_{k,q} - r_{r,q} )^{2} \delta _{ij} - ( {}_{p}{r}_{k,i} - r_{r,i} ) ( {}_{p}{r}_{k,j} - r_{r,j} ) \right] m_{r} \right] \end{aligned}$$

(21)

with mass $m_{r}$ and mass moment of inertia $I_{r}$ of a rigid particle r, cf. Remark 13, and Kronecker delta $\delta _{ij}$, cf. Remark 14. The computed processor-wise quantities, i.e., mass ${}_{p}{m}_{k}$, center of mass position ${}_{p}{\mathbf {r}}_{k}$, and mass moment of inertia ${}_{p}{\mathbf {I}}_{k}$, are communicated to the owning processor of rigid body k. In a second step, on the owning processor the total mass $m_{k}$ and center of mass position $\mathbf {r}_{k}$ of rigid body k are computed over all processors p as

$$\begin{aligned} m_{k} = \sum _{p} {}_{p}{m}_{k} {{\quad }\text {and}\quad }\mathbf {r}_{k} = \frac{\sum _{p} {}_{p}{m}_{k} {}_{p}{\mathbf {r}}_{k}}{\sum _{p} {}_{p}{m}_{k}} \end{aligned}$$

(22)

making use of the received processor-wise quantities. Similar to (21) the mass moment of inertia $\mathbf {I}_{k}$ of rigid body k follows componentwise (in index notation) as

$$\begin{aligned} I_{k,ij} = \sum _{p} \left[ {}_{p}{I}_{k,ij} + \left[ \sum _{q} ( r_{k,q} - {}_{p}{r}_{k,q} )^{2} \delta _{ij} - ( r_{k,i} - {}_{p}{r}_{k,i} ) ( r_{k,j} - {}_{p}{r}_{k,j} ) \right] {}_{p}{m}_{k} \right] \end{aligned}$$

(23)

again considering the received processor-wise quantities. Finally, the determined global quantities, i.e., mass $m_{k}$, center of mass position $\mathbf {r}_{k}$, and mass moment of inertia $\mathbf {I}_{k}$, are communicated from the owning processor to all hosting processors of rigid body k.

Fig. 4

Parallel distribution of a rigid body k with rigid particles r over several processors p illustrating the evaluation of mass $m_{k}$, center of mass position $\mathbf {r}_{k}$, and mass moment of inertia $\mathbf {I}_{k}$ via processor-wise mass ${}_{p}{m}_{k}$, center of mass position ${}_{p}{\mathbf {r}}_{k}$, and mass moment of inertia ${}_{p}{\mathbf {I}}_{k}$

Remark 13

The mass moment of inertia $I_{r}$ of a rigid particle r with mass $m_{r}$ is computed based on the effective volume $V_{eff}=(\Delta {}x)^{d}$ with initial particle spacing $\Delta {}x$. In two-dimensional space ($d=2$) assuming circular disk-shaped particles results in $I_{r} = 0.5 m_{r} r_{eff}$ with effective radius $r_{eff} = {\Delta {}x}/{\sqrt{\pi }}$. Accordingly, in three-dimensional space ($d=3$) assuming spherical-shaped particles results in $I_{r} = 0.4 m_{r} r_{eff}$ with effective radius $r_{eff} = \root 3 \of {{0.75}/{\pi }} \Delta {}x$.

Remark 14

The Kronecker delta $\delta _{ij}$ used in Eqs. (21) and (23) to compute the mass moments of inertia is defined by $\delta _{ij} = {{\left\{ \begin{array}{ll} 1 \quad \text {if}\quad i=j \, , \\ 0 \quad \text {otherwise.}\quad \end{array}\right. }}$

Remark 15

The computation of the mass moment of inertia, cf. Eqs. (21) respectively (23), is based on the Huygens-Steiner theorem, also called parallel axis theorem.

Parallel evaluation of resultant force and torque

To begin with, the resultant coupling and contact force acting on a rigid particle r of rigid body k is given as

$$\begin{aligned} \mathbf {f}_{r} = \sum _{i} \mathbf {f}_{ri} + \sum _{\hat{k}} \sum _{\hat{r}} \mathbf {f}_{r\hat{r}} \end{aligned}$$

(24)

with coupling forces $\mathbf {f}_{ri}$ stemming from interaction with neighboring fluid particles i, cf. Eq. (19), and contact forces $\mathbf {f}_{r\hat{r}}$ stemming from interaction with rigid particles $\hat{r}$ of contacting rigid bodies $\hat{k}$, cf. Eq. (27). Similar to the computation of mass-related quantities as described previously, the resultant force $\mathbf {f}_{k}$ and torque $\mathbf {m}_{k}$ acting on a rigid body k are determined considering the parallel distribution of the affiliated rigid particles r on hosting processors p, cf. Fig. 4. Thus, in a first step the processor-wise resultant force ${}_{p}{\mathbf {f}}_{k}$ and torque ${}_{p}{\mathbf {m}}_{k}$ acting on rigid body k are computed as

$$\begin{aligned} {}_{p}{\mathbf {f}}_{k} = \sum _{r} \mathbf {f}_{r} {{\quad }\text {and}\quad }{}_{p}{\mathbf {m}}_{k} = \sum _{r} \mathbf {r}_{rk} \times \mathbf {f}_{r} \end{aligned}$$

(25)

with $\mathbf {r}_{rk} = \mathbf {r}_{r} - \mathbf {r}_{k}$ while considering the resultant forces $\mathbf {f}_{r}$ acting on all rigid particles r being located in the computational domain of processor p. For correct computation of the processor-wise resultant torque ${}_{p}{\mathbf {m}}_{k}$ the knowledge of the global center of mass position $\mathbf {r}_{k}$ is required on all processors. Finally, the computed processor-wise forces ${}_{p}{\mathbf {f}}_{k}$ and torques ${}_{p}{\mathbf {m}}_{k}$ are communicated to the owning processor of rigid body k and summed up to the global resultant force and torque acting on rigid body k

$$\begin{aligned} \mathbf {f}_{k} = \sum _{p} {}_{p}{\mathbf {f}}_{k} {{\quad }\text {and}\quad }\mathbf {m}_{k} = \sum _{p} {}_{p}{\mathbf {m}}_{k} \, . \end{aligned}$$

(26)

Remark 16

In the case of a computation on a single processor, the evaluation of mass $m_{k}$, center of mass position $\mathbf {r}_{k}$, and mass moment of inertia $\mathbf {I}_{k}$ of a rigid body k follow directly from Eqs. (20) and (21), while the resultant force $\mathbf {f}_{k}$ and torque $\mathbf {m}_{k}$ directly follow from Eq. (25), in each case without the need for special communication.

Contact evaluation between neighboring rigid bodies

For the modeling of frictionless contact between neighboring rigid bodies k and $\hat{k}$, a contact normal force law based on a spring-dashpot model, similar to [29], is employed. The contact force is acting between pairs of neighboring rigid particles r and $\hat{r}$ of contacting rigid bodies, i.e., for distances $r_{r\hat{r}} < \Delta {}x$ with $r_{r\hat{r}} = |\mathbf {r}_{r\hat{r}}| = |\mathbf {r}_{r} - \mathbf {r}_{\hat{r}}|$. Accordingly, the contact force acting on a particle r of rigid body k due to contact with a particle $\hat{r}$ of neighboring rigid body $\hat{k}$ is given as

$$\begin{aligned} \mathbf {f}_{r\hat{r}} = {\left\{ \begin{array}{ll} - \Big [ \min \Big (0, k_{c} ( r_{r\hat{r}} - \Delta {}x ) + d_{c} ( \mathbf {e}_{r\hat{r}} \cdot \mathbf {u}_{r\hat{r}} ) \Big ) \Big ] \, \mathbf {e}_{r\hat{r}} &{} \quad \text {if}\quad r_{r\hat{r}} < \Delta {}x \, , \\ 0 &{} \quad \text {otherwise,}\quad \end{array}\right. } \end{aligned}$$

(27)

with unit vector $\mathbf {e}_{r\hat{r}} = {\mathbf {r}_{r\hat{r}}}/{r_{r\hat{r}}}$, stiffness constant $k_{c}$, and damping constant $d_{c}$. The $\min$ operator in Eq. (27) ensures that only repulsive forces between the rigid particles are considered, also known as tension cut-off.

Remark 17

Contact between a rigid body k and a rigid wall is modeled similar to Eq. (27) considering rigid particles r of a rigid body k and boundary particles b of a discretized rigid wall.

Remark 18

The applied contact evaluation between rigid bodies is for simplicity based on a contact normal force law evaluated between rigid particles while neglecting frictional effects. Generally, following a macroscopic approach of contact mechanics with non-penetration constraint, the normal distance between the contacting bodies, typically determined via closest point projections, is the contact-relevant kinematic quantity. Accordingly, the concept applied in this work, can be interpreted as a microscale approach based on a repulsive/steric interaction potential [51] defined between pairs of rigid particles of contacting rigid bodies. In the current work, this approach has been chosen for reasons of simplicity and numerical robustness. An extension to a macroscale approach, i.e., a normal distance-based contact interaction [12, 52, 53], is possible in a straightforward manner. In addition, also a momentum-based energy tracking method [54] was recently applied for collision modeling of fully resolved rigid bodies [55, 56].

Discretization of the heat equation using SPH

Thermal conduction in the combined fluid and solid domain governed by the heat equation (8) is discretized using smoothed particle hydrodynamics following a formulation proposed by Cleary and Monaghan [57]

$$\begin{aligned} c_{p,a} \frac{\mathrm{d}T_{a}}{\mathrm{d}t} = \frac{1}{\rho _{a}} \sum _{b} V_{b} \frac{4\kappa _{a} \kappa _{b}}{\kappa _{a} + \kappa _{b}} \frac{T_{ab}}{r_{ab}} \frac{{\partial }W}{{\partial }r_{ab}} \end{aligned}$$

(28)

with volume $V_{b} = {m_{b}}/{\rho _{b}}$ of particle b and temperature difference $T_{ab} = T_{a} - T_{b}$ between particle a and particle b. The discretization of the conductive term is especially suited for problems involving a different thermal conductivity among the fields [57]. In the equation above, the index $(\cdot )^{\phi }$ with $\phi \in \{f,s\}$ for fluid and solid field is dropped for ease of notation. Accordingly, the particles a and b may denote fluid particles i as well as rigid particles r, respectively.

Modeling thermally driven reversible phase transitions

Due to the Lagrangian nature of SPH, each (material) particle carries its phase information. This allows for direct evaluation of the discretized heat equation (28) for fluid and rigid particles with corresponding phase-specific parameters of the particle a itself and of neighboring particles b. Phase transitions in the form of melting of a rigid body occurs, in case the temperature $T_{r}$ of a rigid particle r exceeds the transition temperature $T_{t}$. The former rigid particle r changes phase to become a fluid particle i. Conversely, phase transitions in form of solidification occurs, in case the temperature $T_{i}$ of a fluid particle i falls below the transition temperature $T_{t}$ and the former fluid particle i becomes a rigid particle r.

Consequently, each time a rigid body k is subject to phase transition, its mass $m_{k}$, center of mass position $\mathbf {r}_{k}$, and mass moment of inertia $\mathbf {I}_{k}$ are updated. In addition, the velocity $\mathbf {u}_{k}$ after phase transition is determined based on quantities prior to phase transition indicated by index $(\cdot )^{\prime }$ as

$$\begin{aligned} \mathbf {u}_{k} = \mathbf {u}_{k}^{\prime } + \varvec{\omega }_{k} \times ( \mathbf {r}_{rk} - \mathbf {r}_{rk}^{\prime }) \end{aligned}$$

(29)

following rigid body motion with (unchanged) angular velocity $\varvec{\omega }_{k}$.

Time integration following a velocity-Verlet scheme

The discretized fluid and solid field are both integrated in time applying an explicit velocity-Verlet time integration scheme in kick-drift-kick form, also denoted as leapfrog scheme, that is of second order accuracy and reversible in time when dissipative effects are absent [36]. Again, for ease of notation, in the following the indices $(\cdot )^{f}$ and $(\cdot )^{s}$ denoting fluid respectively solid quantities are dropped. Altogether, for the fluid field the positions $\mathbf {r}_{i}$ of fluid particles i are evolved in time, while for the solid field the center of mass positions $\mathbf {r}_{k}$ and the orientations $\varvec{\psi }_{k}$ of all rigid bodies k are evolved in time. However, the positions $\mathbf {r}_{r}$ of rigid particles r are not evolved in time but directly follow the motion of corresponding affiliated rigid bodies k.

In a first kick-step, the accelerations $\mathbf {a}_{i}^{n} = (\mathrm{d}{\mathbf {u}_{i}}/\mathrm{d}{t})^{n}$, as determined in the previous time step n, are used to compute the intermediate velocities

$$\begin{aligned} \mathbf {u}_{i}^{n+1/2} = \mathbf {u}_{i}^{n} + \frac{\Delta {}t}{2} \, \mathbf {a}_{i}^{n} \end{aligned}$$

(30)

of fluid particles i, where $\Delta {}t$ is the time step size. Similar, for rigid bodies k the linear and angular accelerations $\mathbf {a}_{k}^{n} = (\mathrm{d}^{2}{\mathbf {r}_{k}}/\mathrm{d}{t}^{2})^{n}$ respectively $\varvec{\alpha }_{k}^{n} = (\mathrm{d}{\varvec{\omega }_{k}}/\mathrm{d}{t})^{n}$ are used to compute the intermediate linear and angular velocities

$$\begin{aligned} \mathbf {u}_{k}^{n+1/2} = \mathbf {u}_{k}^{n} + \frac{\Delta {}t}{2} \, \mathbf {a}_{k}^{n} {{\quad }\text {and}\quad }\varvec{\omega }_{k}^{n+1/2} = \varvec{\omega }_{k}^{n} + \frac{\Delta {}t}{2} \, \varvec{\alpha }_{k}^{n} \, . \end{aligned}$$

(31)

In a drift-step, the positions (and orientations) of fluid particles i and rigid bodies k are updated to time step $n+1$ using the intermediate velocities. Accordingly, the positions of fluid particles i follow as

$$\begin{aligned} \mathbf {r}_{i}^{n+1} = \mathbf {r}_{i}^{n} + \Delta {}t \, \mathbf {u}_{i}^{n+1/2} \end{aligned}$$

(32)

and the center of mass positions of rigid bodies k as

$$\begin{aligned} \mathbf {r}_{k}^{n+1} = \mathbf {r}_{k}^{n} + \Delta {}t \, \mathbf {u}_{k}^{n+1/2} \, . \end{aligned}$$

(33)

The orientations of rigid bodies k are updated making use of quaternion algebra. First, the angular orientation increments from time step n to time step $n+1$ are determined using the intermediate angular velocities of rigid bodies k following

$$\begin{aligned} \varvec{\phi }_{k}^{n,n+1} = \Delta {}t \, \varvec{\omega }_{k}^{n+1/2} \, . \end{aligned}$$

(34)

Next, the angular orientation increments are described by so-called transition quaternions $\mathbf {q}_{k}^{n,n+1}$. Finally, quaternion multiplication, cf. Remark 12, gives the updated orientations of rigid bodies k at time step $n+1$

$$\begin{aligned} \mathbf {q}_{k}^{n+1} = \mathbf {q}_{k}^{n,n+1} \circ \mathbf {q}_{k}^{n} \, . \end{aligned}$$

(35)

Once the updated orientations (and thus also the updated rigid body frames) are known, the relative positions of rigid particles $\mathbf {r}_{rk}^{n+1}$ can be transformed from the rigid body frame to the reference frame. The velocities and the positions of rigid particles r are updated, considering the underlying rigid body motion of the corresponding rigid bodies k, in consistency with the applied time integration scheme following

$$\begin{aligned} \mathbf {u}_{r}^{n+1/2} = \mathbf {u}_{k}^{n+1/2} + \varvec{\omega }_{k}^{n+1/2} \times \mathbf {r}_{rk}^{n+1} {{\quad }\text {and}\quad }\mathbf {r}_{r}^{n+1} = \mathbf {r}_{k}^{n+1} + \mathbf {r}_{rk}^{n+1} \, . \end{aligned}$$

(36)

Using the positions $\mathbf {r}_{a}^{n+1}$ and the intermediate velocities $\mathbf {u}_{a}^{n+1/2}$ of fluid and rigid particles $a \in \{i,r\}$, the densities $\rho _{i}^{n+1}$ of fluid particles i are computed via Eq. (14). The densities $\rho _{r}$ of rigid particles r are not evolved and remain constant. The temperature rates $(\mathrm{d}{T_{a}}/\mathrm{d}{t})^{n+1}$ of fluid and rigid particles a are then updated on the basis of Eq. (28) with the temperatures $T_{a}^{n}$ as well as the positions $\mathbf {r}_{a}^{n+1}$ and densities $\rho _{a}^{n+1}$. Finally, the temperatures of fluid and rigid particles a are computed as

$$\begin{aligned} T_{a}^{n+1} = T_{a}^{n} + \Delta {}t \left( \frac{\mathrm{d}T_{a}}{\mathrm{d}t}\right) ^{n+1} \, . \end{aligned}$$

(37)

The accelerations $\mathbf {a}_{i}^{n+1}$ of fluid particles i, cf. Eq. (15), and the forces $\mathbf {f}_{r}^{n+1}$ acting on rigid particles r, cf. Eq. (24), are concurrently computed using the positions $\mathbf {r}_{a}^{n+1}$, the intermediate velocities $\mathbf {u}_{a}^{n+1/2}$, and the densities $\rho _{a}^{n+1}$ of fluid and rigid particles a. Consequently, the resultant forces $\mathbf {f}_{k}^{n+1}$ and torques $\mathbf {m}_{k}^{n+1}$ acting on rigid bodies k together with mass-related quantities give the linear and angular accelerations $\mathbf {a}_{k}^{n+1}$ respectively ${\varvec{\alpha }}_{k}^{n+1}$, cf. Eqs. (6) and (7). In a final kick-step, the velocities of fluid particles i at time step $n+1$ are computed as

$$\begin{aligned} \mathbf {u}_{i}^{n+1} = \mathbf {u}_{i}^{n+1/2} + \frac{\Delta {}t}{2} \, \mathbf {a}_{i}^{n+1} \, , \end{aligned}$$

(38)

while the linear and angular velocities of rigid bodies k are

$$\begin{aligned} \mathbf {u}_{k}^{n+1} = \mathbf {u}_{k}^{n+1/2} + \frac{\Delta {}t}{2} \, \mathbf {a}_{k}^{n+1} {{\quad }\text {and}\quad }\varvec{\omega }_{k}^{n+1} = \varvec{\omega }_{k}^{n+1/2} + \frac{\Delta {}t}{2} \, \varvec{\alpha }_{k}^{n+1} \, . \end{aligned}$$

(39)

Accordingly, the velocities of rigid particles r are determined following the motion of the corresponding rigid bodies k to

$$\begin{aligned} \mathbf {u}_{r}^{n+1} = \mathbf {u}_{k}^{n+1} + \varvec{\omega }_{k}^{n+1} \times \mathbf {r}_{rk}^{n+1} \, . \end{aligned}$$

(40)

To maintain stability of the time integration scheme, the time step size $\Delta {}t$ is restricted by the Courant–Friedrichs–Lewy (CFL) condition, the viscous condition, the body force condition, the contact condition, and the conductivity condition, refer to [21, 50, 58, 59] for more details,

$$\begin{aligned} \Delta {}t \le \min \left\{ 0.25\frac{h}{c+|\mathbf {u}_{max}|}, 0.125\frac{h^{2}}{\nu }, 0.25\sqrt{\frac{h}{|\mathbf {b}_{max}|}}, 0.22 \sqrt{\frac{m_{r}}{k_{c}}}, 0.1\frac{\rho c_{p} h^{2}}{\kappa } \right\} \, , \end{aligned}$$

(41)

with maximum fluid velocity $\mathbf {u}_{max}$ and maximum body force $\mathbf {b}_{max}$.

Numerical examples

The purpose of this section is to investigate the proposed numerical formulation for solving fluid–solid and contact interaction problems examining several numerical examples in two and three dimensions involving multiple mobile rigid bodies, two-phase flow, and reversible phase transitions. To begin with, several numerical examples of a single rigid body in a fluid flow, considering different spatial discretizations, are studied and compared to reference solutions. In a next step, two examples close to potential application scenarios of the proposed formulation in the fields of engineering and biomechanics are investigated. Finally, the capabilities of the proposed parallel computational framework are demonstrated performing a strong scaling analysis. The parameter values in the numerical examples, unless indicated otherwise, are given in a consistent set of units and presented in non-dimensional form.

Spatial discretization of a rigid circular disk

In the following, a rigid circular disk of diameter $D = 2.5 \times 10^{-3}$ with density $\rho ^{s} = 1.0 \times 10^{3}$, motivated by two subsequent examples, is discretized with different values of the initial particle spacing $\Delta {}x$. The mass $m^{s}$ and the mass moment of inertia $I^{s}$ (with respect to the axis of symmetry) of the circular disk are computed with the proposed formulation and shown in Fig. 5. With decreasing initial particle spacing $\Delta {}x$ the values for mass $m^{s}$ and mass moment of inertia $I^{s}$ converge to the analytical solution confirming the proposed formulation. To illustrate, the resulting spatial discretizations of the circular disk with rigid particles are shown in Fig. 6. Clearly, the approximation of the circular shape of the disk is of better accuracy for decreasing initial particle spacing $\Delta {}x$. To keep the computational effort at a feasible level, in the two subsequent examples the domain is discretized with an initial particle spacing of $\Delta {}x = 2.0 \times 10^{-4}$ and $\Delta {}x = 1.0 \times 10^{-4}$.

Fig. 5

Spatial discretization of a rigid circular disk: mass $m^{s}$ and mass moment of inertia $I^{s}$ of a rigid circular disk for different values of the initial particle spacing $\Delta {}x$ (solid line) compared to the analytical solution (dashed line)

Fig. 6

Spatial discretization of a rigid circular disk: resulting spatial discretization of a rigid circular disk with rigid particles for selected values of the initial particle spacing $\Delta {}x$

A rigid circular disk floating in a shear flow

The following numerical examples are concerned with the motion of a rigid circular disk floating in a shear flow. First, the principal setup of the problem along with numerical parameters is described, thereafter, two distinct cases are considered in detail. For validation, the results obtained with the proposed formulation are compared to [8] also applying SPH to discretize the fluid and the solid field.

A rigid circular disk of diameter $D = 2.5 \times 10^{-3}$ with density $\rho ^{s} = 1.0 \times 10^{3}$ is allowed to move freely in a rectangular channel of length $L = 5.0 \times 10^{-2}$ and height $H = 1.0 \times 10^{-2}$, cf. Fig. 7. The remainder of the channel is occupied by a Newtonian fluid with density $\rho ^{f} = 1.0 \times 10^{3}$ and kinematic viscosity $\nu ^{f} = 5.0 \times 10^{-6}$. The bottom and top channel walls move with velocity ${u_{w}}/{2}$ in opposite direction inducing a shear flow in the channel. The Reynolds number of the problem is given as $Re = {u_{w} D^{2}}/{4 \nu ^{f} H}$ [8, 60] taking into account the diameter of the circular disk D and the channel height H. At the left and right end of the channel, periodic boundary conditions are applied, cf. Remark 19.

For the fluid phase, an artificial speed of sound $c = 0.25$ is chosen, resulting in a reference pressure $p_{0} = 62.5$ of the weakly compressible model. The background pressure $p_{b}$ of the transport velocity formulation is set equal to the reference pressure $p_{0}$. The motion of the bottom and top channel walls is modeled using moving boundary particles. The problem is solved for different values of the initial particle spacing $\Delta {}x$ for times $t \in [0, 60.0]$ with time step size $\Delta {}t$ obeying respective conditions (41).

Remark 19

Imposing a periodic boundary condition in a specific spatial direction allows for particle interaction evaluation across opposite domain borders. Moreover, particles leaving the domain on one side are re-entering on the opposite side.

Fig. 7

A rigid circular disk floating in a shear flow: geometry and boundary conditions of two different cases

Case 1: Migration of a floating rigid circular disk to the center line of a channel

This case is based on studies [60, 61] stating that a rigid circular disk floating in a shear flow in a channel migrates to the center line of the channel independent of its initial position and initial velocity. Herein, the rigid circular disk is initially at rest placed at vertical position $r_{y} = 2.5 \times 10^{-3}$ in the channel, cf. Fig. 7. The channel walls move in opposite direction with a velocity magnitude of ${u_{w}}/{2} = 0.01$ resulting in the Reynolds number $Re = 0.625$ of the problem.

The obtained vertical position $r_{y}$ and the horizontal velocity $u_{x}$ of the center of the circular disk in the channel over time t are displayed in Fig. 8 for two different values of the initial particle spacing $\Delta {}x$. The circular disk migrates to the center line of the channel as expected, showing no significant difference between the results obtained with different initial particle spacings $\Delta {}x$. In addition, a comparison to the results of [8] shows very good agreement for the dynamics of the solution.

Fig. 8

Migration of a floating rigid circular disk to the center line of a channel: vertical position $r_{y}$ and horizontal velocity $u_{x}$ of the center of the floating circular disk in the channel computed with the proposed formulation and an initial particle spacing of $\Delta {}x = 2.0 \times 10^{-4}$ (red dashed line) and $\Delta {}x = 1.0 \times 10^{-4}$ (black solid line) compared to the reference solution [8] (crosses)

Case 2: Interaction of a floating rigid circular disk with a fixed rigid circular disk

In the presence of a rigid circular disk that is fixed at the center line of the channel, a rigid circular disk floating in a shear flow migrates to a specific position of equilibrium independent of its initial position and velocity as stated in [62]. Herein, the fixed and the floating rigid circular disks are initially placed on the center line of the channel at horizontal position $r_{x} = \pm 3.75 \times 10^{-3}$, cf. Fig. 7. The channel walls move in opposite direction with a velocity magnitude of ${u_{w}}/{2} = 0.012$ resulting in the Reynolds number $Re = 0.75$ of the problem.

Figure 9 shows the obtained trajectory, i.e., vertical position $r_{y}$ over horizontal position $r_{x}$, and horizontal velocity $u_{x}$ of the center of the floating circular disk in the channel for two different values of the initial particle spacing $\Delta {}x$. The results obtained with initial particle spacing $\Delta {}x = 1.0 \times 10^{-4}$ are in good agreement to the reference solution [8]. However, the results obtained with initial particle spacing $\Delta {}x = 2.0 \times 10^{-4}$ show fluctuations of the horizontal velocity $u_{x}$, which is why also the trajectory deviates from the reference solution [8]. This can be explained with disturbances of the density field due to relative particle movement [21], that are more pronounced with a coarser spatial discretization, i.e., with larger initial particle spacing $\Delta {}x$.

Fig. 9

Interaction of a floating rigid circular disk with a fixed rigid circular disk: trajectory and horizontal velocity $u_{x}$ of the center of the floating circular disk in the channel computed with the proposed formulation and an initial particle spacing of $\Delta {}x = 2.0 \times 10^{-4}$ (red dashed line) and $\Delta {}x = 1.0 \times 10^{-4}$ (black solid line) compared to the reference solution [8] (crosses)

A rigid circular disk falling in a fluid column

A rigid circular disk of diameter $D = 2.5 \times 10^{-3}$ with density $\rho ^{s} = 1.25 \times 10^{3}$ is initially at rest placed on the y-axis at vertical position $r_{y} = 1.0 \times 10^{-2}$ in a closed rectangular box of height $H = 6.0 \times 10^{-2}$ and width $W = 2.0 \times 10^{-2}$, cf. Fig. 10. The remainder of the box is occupied by a Newtonian fluid with density $\rho ^{f} = 1.0 \times 10^{3}$ and kinematic viscosity $\nu ^{f} = 1.0 \times 10^{-5}$. A gravitational acceleration of magnitude $|\mathbf {g}| = 9.81$ shall act on both the fluid and solid field in negative y-direction. Following [8] this is modeled considering the buoyancy effect, i.e., body force $\mathbf {b}^{s} = { ( \rho ^{s} - \rho ^{f} ) }/{\rho ^{s}} \mathbf {g}$ is acting on the solid field while no body force $\mathbf {b}^{f} = 0.0$ is applied on the fluid field (each per unit mass). It is worth noting that, naturally, it would also be possible to directly set the gravitational acceleration for the fluid and solid field. For validation, the results obtained with the proposed formulation are compared to [8] also applying SPH to discretize the fluid and the solid field.

Fig. 10

A rigid circular disk falling in a fluid column: geometry and boundary conditions of the problem

For the fluid phase, an artificial speed of sound $c = 0.5$ is chosen, resulting in a reference pressure $p_{0} = 250.0$ of the weakly compressible model. The background pressure $p_{b}$ of the transport velocity formulation is set equal to the reference pressure $p_{0}$. The stiffness and damping constant applied for contact evaluation are set to $k_{c} = 1.0 \times 10^{8}$ and $d_{c} = 1.0 \times 10^{2}$. The walls of the box are modeled using boundary particles. The problem is solved for different values of the initial particle spacing $\Delta {}x$ for times $t \in [0, 0.8]$ with time step size $\Delta {}t$ obeying respective conditions (41).

The obtained vertical velocity and horizontal position of the center of the circular disk in the box over time t are displayed in Fig. 11 for two different values of the initial particle spacing $\Delta {}x$ compared to the reference solution [8]. The results obtained with different initial particle spacing $\Delta {}x$ show only minor differences. The terminal velocity of the rigid circular disk is slightly smaller than given in the reference solution [8]. It shall be noted that, in contrast to [8], contact of the rigid circular disk and the wall of the box is explicitly considered. Consequently, the rigid circular disk comes at rest when approaching the bottom wall of the box.

Fig. 11

A rigid circular disk falling in a fluid column: vertical position $r_{y}$ and vertical velocity $u_{y}$ of the center of the circular disk in the fluid column computed with the proposed formulation and an initial particle spacing of $\Delta {}x = 2.0 \times 10^{-4}$ (red dashed line) and $\Delta {}x = 1.0 \times 10^{-4}$ (black solid line) compared to the reference solution [8] (crosses)

Melting and solidification of powder grains in a melt pool

In metal powder bed fusion additive manufacturing (PBFAM), structural components are created utilizing a laser or electron beam to melt and fuse metal powder, layer per layer, to form the final part. PBFAM has the potential to enable new paradigms of product design, manufacturing and supply chains. However, due to the complexity of PBFAM processes, the interplay of process parameters is not completely understood, creating the need for further research, amongst others in the field of computational melt pool modeling [15, 16]. For this purpose, an SPH formulation for thermo-capillary phase transition problems with a focus on metal PBFAM melt pool modeling has recently been proposed [18]. For simplicity, this and other state-of-the-art approaches [18,19,20] in the field consider powder grains that are spatially fixed. In the real physical process, however, it is observed that, depending on the processing conditions, melt evaporation and thereby induced vapor and gas flows in the build chamber may result in powder grains entrainment and ejection, i.e., a considerable degree of material re-distribution during the melting process. On the one hand, this effect considerably affects process stability and mechanisms of defect creation, on the other hand, it can not be represented by state-of-the-art approaches restricted to immobile powder grains [18,19,20]. In the following, a two- and a three-dimensional example each with different geometry, boundary conditions, and material parameters are considered. The purpose of these examples is not to study PBFAM in detail but to showcase the general applicability of the proposed formulation to capture the dynamics of mobile powder grains undergoing temperature-induced phase transitions, i.e., melting and solidification, while being exposed to gravitational acceleration and gas flow. To this end, surface tension and wetting effects as well as the influence of evaporation-induced recoil pressure, as discussed in [18], are neglected. The focus is set on the investigation of highly dynamic motion and interaction of powder grains with each other, the liquid melt phase and a surrounding gas phase, undergoing reversible phase transitions, i.e., melting and solidification.

Powder grains exposed to a highly dynamic gas flow

This two-dimensional example is solely intended to demonstrate the model capabilities for this type of application, while keeping the overall example simple in this method-focused contribution. For this reason, non-physical parameter values and boundary conditions are chosen in the following.

A rectangular box is composed of two chambers, each with width 20.0 and height 12.0, that are connected by an opening spanning the upper half of the box. An inlet and outlet of width 3.0 are located at the top left and top right end of the box. Powder grains of a solid metal phase (density $\rho ^{s} = 1.0$, heat capacity $c_{p}^{s} = 1.0$, thermal conductivity $\kappa ^{s} = 10.0$) with diameters between 2.5 and 4.4 are placed initially at rest inside the left chamber of the box. The initial positions of the powder grains can, e.g., be obtained in a pre-processing step based on the discrete element method (DEM) and a cohesive powder model [29, 63]. The remainder of the box is initially filled with a gas phase (Newtonian fluid, density $\rho ^{g} = 0.1$, kinematic viscosity $\nu ^{g} = 100.0$, heat capacity $c_{p}^{g} = 0.01$, thermal conductivity $\kappa ^{g} = 0.1$). The temperature is initialized to $T^{s}_{0} = 25.0$ within the solid metal phase and to $T^{g}_{0} = 50.0$ within the gas phase. In the upper half of the box walls the temperature is fixed to $\hat{T} = 50.0$ at all times. In the lower half of the box walls the temperature is set to $\hat{T} = 100.0$ until time $t \le 0.5$, and to $\hat{T} = 0.0$ for time $t > 0.5$. Refer to Figs. 12 and 13 containing an illustration of the initial configuration. Reversible phase transitions between solid metal phase and liquid metal phase (Newtonian fluid, density $\rho ^{l} = 1.0$, kinematic viscosity $\nu ^{l} = 100.0$, heat capacity $c_{p}^{l} = 1.0$, thermal conductivity $\kappa ^{l} = 10.0$) is assumed to occur at a transition temperature of $T_{t} = 50.0$. With the goal to evoke drag forces acting on the powder grains, for times $t > 0.25$ a parabolic inflow respectively outflow of the gas phase with mean velocity 420.0 is prescribed at the inlet and outlet of the box. A gravitational acceleration of magnitude $|\mathbf {g}| = 1.0 \times 10^{4}$ is acting downwards, set as body force (per unit mass) of all involved phases.

Fig. 12

Powder grains exposed to a highly dynamic gas flow: time series of the obtained results with temperature field ranging from 0.0 (blue) to 100.0 (red)

Fig. 13

Powder grains exposed to a highly dynamic gas flow: time series of the obtained results with magnitude of the velocity field ranging from 0.0 (blue) to 600.0 (red)

For both fluid phases (liquid metal and gas), the reference pressure of the weakly compressible model is set to $p_{0} = 16.0 \times 10^{6}$, and the background pressure $p_{b}$ of the transport velocity formulation is set equal to the reference pressure $p_{0}$. The stiffness and damping constant applied for contact evaluation are set to $k_{c} = 1.0 \times 10^{8}$ and $d_{c} = 1.0 \times 10^{2}$. The wall of the box is modeled using boundary particles. The inflow and outflow conditions are modeled similar as described in [49]. The problem is solved with initial particle spacing $\Delta {}x = 0.1$ for times $t \in [0, 0.75]$ with a time step size of $\Delta {}t = 0.625 \times 10^{-5}$.

A time series of illustrations of the obtained results is given in Figs. 12 and 13. The solid metal phase is visualized in gray color. The particles discretizing the liquid metal phase are displayed in black color. In the background, the temperature respectively velocity field of the combined liquid metal and gas phase are displayed. Thereto, both fields were post-processed applying SPH approximation (13) and visualized by a color code. In Fig. 12 additionally the temperature of the walls is shown. First, the powder grains are heated and gradually start melting into liquid metal where in close contact to the hot wall. Eventually, after time $t=0.25$, powder grains are subjected to the gas flow through the box. Some (partially melted) powder grains are swept into the right chamber of the box, where melting after contact with the hot wall continues. A non-smooth and strongly distorted interface topology between liquid metal and gas phase develops, especially in the right chamber of the box, because surface tension and wetting effects are neglected. Finally, with the temperature in the lower half of the box set to $\hat{T} = 0.0$ after time $t = 0.5$, the liquid metal phase is cooled down drastically and eventually resolidifies.

Three-dimensional powder melting setup with representative material parameters

In this example, the focus is set on a three-dimensional setup while considering representative material parameters for stainless steel and the surrounding gas phase taken from [18] with the purpose to demonstrate the general applicability of the proposed formulation for metal PBFAM melt pool modeling. However, note that some characteristic phenomena relevant for the latter, e.g., surface tension and wetting effects, are still neglected.

Powder grains of a solid metal phase (density $\rho ^{s} = {7430}\, \mathrm{kg}\, \mathrm{m}^{-3}$, heat capacity $c_{p}^{s} = {965}\,\mathrm{J}\, \mathrm{kg}^{-1}\, \mathrm{K}^{-1}$, thermal conductivity $\kappa ^{s} = {35.95}\,\mathrm{W \, m}^{-1}\, \mathrm{K}^{-1}$) with diameters between ${16}\,{\upmu }\mathrm{m}$ and ${20}\,{\upmu }\mathrm{m}$ are placed initially at rest inside a cuboid box with dimensions of ${50}\,{\upmu }\mathrm{m} \times {50}\,{\upmu }\mathrm{m} \times {30}\,{\upmu }\mathrm{m}$. The initial positions of the powder grains can, e.g., be obtained in a pre-processing step based on the discrete element method (DEM) and a cohesive powder model [29, 63]. The remainder of the box is initially filled with a gas phase (Newtonian fluid, density $\rho ^{g} = {74.3}\, \mathrm{kg\, m}^{-3}$, dynamic viscosity $\eta ^{g} = {6.0e-4} \,\mathrm{kg \, m}^{-1} \mathrm{s}^{-1}$, heat capacity $c_{p}^{g} = {10}\, \mathrm{J \, kg}^{-1} \mathrm{K}^{-1}$, thermal conductivity $\kappa ^{g} = {2.6}\times 10^{-2}\,\mathrm{W \, m}^{-1}\mathrm{K}^{-1}$). The temperature is initialized to $T^{s}_{0} = {500}\,\mathrm{K}$ and $T^{g}_{0} = {500}\,\mathrm{K}$ within the solid metal phase and the gas phase. In all box walls the temperature is set to $\hat{T} = {1800}\,\mathrm{K}$ until time $t \le {0.15}\,\mathrm{ms}$, and to $\hat{T} = {500}\,\mathrm{K}$ for time $t > {0.15}\,\mathrm{ms}$. Refer to Fig. 14 containing an illustration of the initial configuration. Reversible phase transitions between solid metal phase and liquid metal phase (Newtonian fluid, density $\rho ^{l} = {7430}\,\mathrm{kg}\, \mathrm{m}^{-3}$, dynamic viscosity $\eta ^{l} = {6.0} \times 10^{-3}\, \mathrm{kg \, m}^{-1} \mathrm{s}^{-1}$, heat capacity $c_{p}^{l} = {965}\,\mathrm{J \, kg}^{-1} \mathrm{K}^{-1}$, thermal conductivity $\kappa ^{l} = {35.95}\,\mathrm{W \, m}^{-1} \mathrm{K}^{-1}$) is assumed to occur at a transition temperature of $T_{t} = {1700}\,\mathrm{K}$. A gravitational acceleration of magnitude $|\mathbf {g}| = {9.81}\,\mathrm{m s}^{-2}$ is acting downwards, set as body force (per unit mass) of all involved phases.

For both fluid phases (liquid metal and gas), the reference pressure of the weakly compressible model is set to $p_{0} = {1.0} \times 10^{7}\, \mathrm{N \, m}^{-2}$ and the background pressure of the transport velocity formulation to $p_{b} = 5 p_{0}$. The stiffness and damping constant applied for contact evaluation are set to $k_{c} = {1.0}\,\mathrm{kg}\, \mathrm{s}^{-2}$ and $d_{c} = {1.0} \times \mathrm{kg \,s}^{-1}$. The wall of the box is modeled using boundary particles. The complete domain is discretized by particles with initial particle spacing $\Delta {}x = {1.0}\,{\upmu }\mathrm{m}$ resulting in a total of approximately $1.13 \times 10^{5}$ particles. The problem is solved for times $t \in [0, 0.175] \, \mathrm{ms}$ with a time step size of $\Delta {}t = {5.0} \times 10^{-7}\mathrm{ms}$.

Fig. 14

Three-dimensional powder melting setup with representative material parameters: time series of the obtained results with temperature field ranging from ${1700}\,\mathrm{K}$ (blue) to ${1800}\,\mathrm{K}$ (red)

A time series of illustrations of the obtained results is given in Fig. 14. The particles discretizing the solid metal phase are displayed in gray color. The particles discretizing the liquid metal phase are colored based on the temperature field with transition temperature $T_{t} = {1700}\,\mathrm{K}$ as lower value and temperature of the hot box walls $\hat{T} = {1800}\,\mathrm{K}$ as upper value. In the background, the temperature field of the combined liquid metal and gas phase is post-processed applying SPH approximation (13) and visualized by a color code until time $t \le {0.15}\,\mathrm{ms}$. The powder grains in close contact to the hot box walls are heated and gradually start melting into liquid metal. Under the influence of gravity the (partially melted) powder grains are gradually displacing the liquid metal. A non-smooth and strongly distorted interface topology between the liquid metal and gas phase develops. This can be explained by the fact that surface tension and wetting effects are neglected in this specific example. After time $t = {0.15}\,\mathrm{ms}$ the temperature of the box walls is suddenly set to $\hat{T} = {500}\, \mathrm{K}$. As a consequence, the liquid metal close to the cool box walls rapidly resolidifies.

In sum, both examples demonstrates that highly dynamic motion of arbitrarily-shaped powder grains as relevant, e.g., for metal PBFAM melt pool modeling, can be captured along with melting and solidification by the proposed formulation in a robust manner. Consequently, the proposed formulation can be recommended as a useful extension of the SPH formulation for mesoscale melt pool modeling [18], or other current state-of-the-art approaches, e.g., [19, 20], allowing for more detailed studies of PBFAM processes.

Gastric disintegration of food boluses

Examination of gastric fluid mechanics plays an important role for modeling digestion of food in the human stomach. The digesta are characterized by a multiphasic nature consisting of fluid (gastric juice and chyme) and solid (food boluses) phases [17]. Intragastric fluid motion is driven by the propagation of so-called antral contraction waves (ACWs), i.e., circular constrictions of the gastric wall due to smooth muscle contractions [64]. The ACWs are initiated at the pacemaker region of the stomach and travel along the greater curvature towards the pylorus both mixing and grinding the digesta. Concurrently, absorption of gastric juice fosters chemical and mechanical breakdown of food boluses into chyme [65]. At low viscosity, i.e., following intragastric dilution of the digesta with gastric juice, retropulsive jet-like fluid motion between the ACWs can be observed [66,67,68].

This example aims to demonstrate the capability of the proposed formulation to replicate typical gastric flow patterns including phase transitions. As compared to this complex application scenario, the configuration of the example is kept simple to focus on the principal effects. Consequently, non-physiological parameter values and boundary conditions are applied. Consider a rectangular box of width 40.0 and height 16.0 (coordinate system in the center) with a mobile constriction. A total of 60 food boluses (density $\rho ^{s} = 1.0$, diffusivity $D^{s} = 0.25$), represented by mobile rigid bodies with diameters between 1.6 and 2.8, are placed at random positions inside the box. The remainder of the box is initially filled with gastric juice (Newtonian fluid, density $\rho ^{g} = 1.0$, kinematic viscosity $\nu ^{g} = 100.0$, diffusivity $D^{g} = 1.0$). Both the food boluses and the gastric juice are initially at rest. The initial configuration of the example is contained in Fig. 15. Over time, food boluses disintegrate into chyme (Newtonian fluid, density $\rho ^{c} = 1.0$, kinematic viscosity $\nu ^{c} = 200.0$, diffusivity $D^{c} = 0.25$). Herein, this is modeled considering the transport of a concentration C within the food boluses and chyme, resembling some kind of moisture penetration, by solving a diffusion equation, cf. Remark 4. Accordingly, the concentration within the food boluses is initialized with $C_{0} = 0.0$, while the concentration within the gastric juice is fixed to $\hat{C} = 1.0$ at all times. Phase transitions from food boluses to chyme is assumed to occur at a transition concentration of $C_{t} = 0.8$. The propagation of an ACW is modeled by the movement of the mobile constriction in the box with a time dependent horizontal velocity of $-14.5 \pi \sin (\pi t)$ from horizontal position 14.5 to $-14.5$, cf. Fig. 15.

For both fluid phases (gastric juice and chyme), an artificial speed of sound $c = 1.0 \times 10^{3}$ is chosen, resulting in a reference pressure $p_{0} = 1.0 \times 10^{6}$ of the weakly compressible model, with background pressure of the transport velocity formulation set to $p_{b} = 5 p_{0}$. The stiffness and damping constant applied for contact evaluation are set to $k_{c} = 1.0 \times 10^{8}$ and $d_{c} = 1.0 \times 10^{2}$. The wall of the box and the mobile constriction are modeled using (moving) boundary particles. The problem is solved with initial particle spacing $\Delta {}x = 0.1$ for times $t \in [0, 1.0]$ with a time step size of $\Delta {}t = 1.25 \times 10^{-5}$.

Figure 15 shows a time series of illustrations of the obtained results. The food boluses are visualized in gray color. The particles discretizing the chyme are displayed in black color. In the background, the velocity field of both gastric juice and chyme is post-processed applying SPH approximation (13) and visualized by a color code. Clearly, the typical retropulsive jet-like fluid motion induced by the moving constriction can be observed. As a consequence, the food boluses are entrained with the fluid flow through the opening while coming into contact with each other. At the same time, disintegration of food boluses into chyme gradually takes place. After time $t = 0.875$ some food boluses are completely dissolved. A detailed view of the region at the mobile constriction is given in Fig. 16 for selected points in time. Here, the particles discretizing the food boluses and the chyme are colored based on the concentration field, for distinction, utilizing two different color maps with transition concentration $C_{t}$ as upper respectively lower value. A progressive mixing of gastric juice and chyme can be observed primarily driven by the fluid motion.

Fig. 15

Gastric disintegration of food boluses: time series of the obtained results with magnitude of the velocity field ranging from 0.0 (blue) to 120.0 (red)

Fig. 16

Gastric disintegration of food boluses: detailed view of region at the mobile constriction for selected points in time. The concentration C is ranging from 0.0 (black) to $C_{t} = 0.8$ (white) within the food boluses and from $C_{t} = 0.8$ (blue) to 1.0 (red) within the chyme. In the background, the magnitude of the velocity field is ranging from 0.0 (blue) to 120.0 (red)

Note that the main purpose of this example is to show the robustness of the proposed formulation in the context of highly dynamic fluid flow and phase transitions, e.g., as occurring in the form of retropulsive jet-like fluid motion during digestion of food in the human stomach. For the sake of simplicity, non-physiological parameter values are applied. Amongst others, the time scales of ACW propagation and disintegration of food boluses are in a mismatch. In addition, the employed phenomenological digestion model does not explicitly resolve the influence of chemical and mechanical breakdown taking place in reality. In conclusion, this example demonstrates that typical gastric flow patterns including phase transitions are fully captured in a stable and robust manner.

Strong scaling analysis of parallel computational framework

The purpose of this example is to demonstrate the capability and efficiency of the proposed parallel computational framework in handling systems constituted of a large number of particles. To this end, a three-dimensional example consisting of a total of approximately $3.79 \times 10^{6}$ particles is examined on two different parallel systems. Conclusions are drawn concerning the parallel behavior of the parallel computational framework in three dimensions.

A total of 216 spherical-shaped mobile rigid bodies with diameter $D = 2.5$ and density $\rho ^{s} = 10.0$ are placed on a regular grid in a cubic box of edge length $L = 30.0$. The rigid bodies are initially at rest and not in contact with each other or the walls of the box. The remainder of the box is occupied by a Newtonian fluid initially at rest with density $\rho ^{f} = 1.0$ and kinematic viscosity $\nu ^{f} = 1.0$. A gravitational acceleration of magnitude $|\mathbf {g}| = 1.0$ is acting in downward direction, i.e., the body forces (per unit mass) of fluid and solid field are given to $\mathbf {b}^{f} = \mathbf {g}$ and $\mathbf {b}^{s} = \mathbf {g}$.

For the fluid phase, an artificial speed of sound $c = 50.0$ is chosen, resulting in a reference pressure $p_{0} = 2.5 \times 10^{3}$ of the weakly compressible model, with background pressure $p_{b}$ of the transport velocity formulation set equal to the reference pressure $p_{0}$. The stiffness and damping constant applied for contact evaluation are set to $k_{c} = 1.0 \times 10^{3}$ and $d_{c} = 1.0 \times 10^{2}$. The wall of the box is modeled using boundary particles. The complete domain is discretized by particles with initial particle spacing $\Delta {}x = 0.2$ resulting in a total of approximately $3.79 \times 10^{6}$ particles, thereof $3.15 \times 10^{6}$ fluid particles, $2.20 \times 10^{5}$ rigid particles, and $4.21 \times 10^{5}$ boundary particles. Following a spatial decomposition approach, the computational domain is divided into $48 \times 48 \times 48$ cubic cells of edge length 0.65 resulting in approximately 34 particles per cell. The problem is solved for times $t \in [0, 30.0]$ with a time step size of $\Delta {}t = 1.0 \times 10^{-3}$.

For the purposes of illustration, the spherical-shaped rigid bodies within the box are shown in Fig. 17 for the initial setup at $t = 0.0$ and later points in time. In addition, the velocity field of the surrounding fluid is post-processed applying SPH approximation (13) and visualized by a color code with opacity. The rigid bodies are falling freely in the viscous fluid under gravity due to density ratio ${\rho ^{s}}/{\rho ^{f}} = 10.0$ until contact with the bottom wall of the box occurs, as first observed after $t \approx 2.5$, or with neighboring rigid bodies, as first observed after $t \approx 5.0$. The rigid bodies begin piling up at the bottom wall of the box and are nearly at rest at $t = 30.0$.

Fig. 17

Strong scaling analysis of parallel computational framework: position of spherical-shaped rigid bodies within the box for different points in time with magnitude of the velocity field ranging from 0.0 (blue) to 1.25 (red)

To showcase the capability and efficiency of the parallel computational framework, a strong scaling analysis is performed utilizing two different parallel systems: The first one consisting of 32 nodes with $2 \times 12$ cores (Intel Xeon E5-2680 v3 Haswell, 2.5 GHz) and the second one consisting of 8 nodes with $2 \times 8$ cores (Intel Xeon E5-2630 v3 Haswell, 2.4 GHz). The parallel behavior of the proposed computational framework is given in Fig. 18, illustrating the obtained solver time per time step and the parallel efficiency given in percent of linear scaling. The parallel efficiency is computed as ${t_{1}}/{ ( n \cdot t_{n} ) } \cdot 100\%$, where $t_{1}$ and $t_{n}$ are the times to solve the problem on one node respectively n nodes.

Fig. 18

Strong scaling analysis of parallel computational framework: solver time per time step (left) and parallel efficiency given in percent of linear scaling (right) for a three-dimensional problem consisting of approximately $3.79 \times 10^{6}$ particles on up to 768 cores

The parallel computational framework scales almost linearly on both parallel systems for up to 128 cores respectively 192 cores. In this regime a parallel efficiency of more than $60\%$ can be observed. For larger numbers of cores, the scalability deteriorates and the parallel efficiency drops to under $50\%$. This can be explained with an increasing communication overhead, cf. Remark 9. Comparable results of a strong scaling analysis for an SPH implementation are given, e.g., in [45] (${r_{c}}/{\Delta {}x} = 2.5$) and [69] (${r_{c}}/{\Delta {}x} = 2.4$), however, in contrast to this example (${r_{c}}/{\Delta {}x} = 3.0$) with a smaller ratio of the support radius $r_{c}$ and the initial particle spacing $\Delta {}x$, resulting in a lower influence on the communication overhead, cf. Remark 9. As a conclusion one can state that the parallel computational framework is capable of efficiently solving systems constituted of a large number of particles in three dimensions on multiple cores. Looking at the parallel behavior, the obtained results confirm that the proposed framework meets all requirements necessary for detailed and accordingly computationally expensive studies.

Conclusion

In this work, an approach for fluid–solid and contact interaction problems including thermo-mechanical coupling and reversible phase transitions is presented. All fields are spatially discretized using smoothed particle hydrodynamics (SPH). Being a mesh-free discretization scheme, SPH is, compared to mesh-based methods, especially suitable in the context of continually changing interface topologies and dynamic phase transitions by avoiding additional methodological and computational effort to capture such phenomena. A detailed concept for the parallelization of the computational framework, especially for an efficient evaluation of rigid body motion, is an essential part of this work.

The accuracy and robustness of the proposed formulation are demonstrated by several numerical examples studying a single rigid body in fluid flow. The obtained numerical results are in very good agreement with the literature. Also two complex examples close to potential applications scenarios in the fields of engineering and biomechanics were studied. First, motivated by metal PBFAM melt pool modeling, melting and solidification of powder grains subject to highly dynamic fluid motion was simulated. Second, inspired by multiphysics modeling of the human stomach, gastric disintegration of food boluses is considered. Both examples confirm that highly dynamic motion of arbitrarily-shaped rigid bodies embedded in a complex fluid flow and including reversible phase transitions can be captured by the proposed framework in a stable and robust manner. Finally, the parallel computing abilities of the proposed computational framework were demonstrated by a strong scaling analysis of a three-dimensional example with $3.79 \times 10^{6}$ particles revealing a parallel efficiency of more than $60\%$ on up to 192 cores.

To the best of the authors’ knowledge, the proposed parallel computational framework is the first of its kind modeling rigid body motion while simultaneously considering thermal conduction, reversible phase transitions, and multiple (liquid and gas) phases. For the sake of simplicity, some characteristic phenomena for thermo-capillary flow and, especially, for metal PBFAM melt pool modeling are not yet addressed in this work, however, straightforward to include considering the authors’ previous work. Besides, the applied contact evaluation between rigid bodies is based on a contact normal force law evaluated between rigid particles while neglecting frictional and adhesive effects. An extension to a more sophisticated contact evaluation is part of ongoing research. In summary, the proposed formulation has the ability to accurately model a host of complex multiphysics problems, and it can thus be expected to become a valuable tool for detailed studies in engineering, e.g., metal additive manufacturing, and biomechanics, e.g., digestion of food in the human stomach.