A Discrete Multi-Physics Model to Simulate Fluid Structure Interaction and Breakage of Capsules Filled with Liquid under Coaxial Load

: This paper investigated the mechanical response (including breakage and release of the internal liquid) of single core–shell capsules under compression by means of discrete multi-physics. The model combined Smoothed Particle Hydrodynamics for modelling the ﬂuid and the Lattice Spring Model for the elastic membrane. Thanks to the meshless nature of discrete multi-physics, the model can easily account for the fracture of the capsule’s shell and the interactions between the internal liquid and the solid shell. The simulations replicated a parallel plate compression test of a single core–shell capsule. The inputs of the model were the size of the capsule, the thickness of the shell, the geometry of the internal structure, the Young’s modulus of the shell material, and the ﬂuid’s density and viscosity. The outputs of the model were the fracture type, the maximum force needed for the fracture, and the force–displacement curve. The data were validated by reproducing equivalent experimental tests in the laboratory. The simulations accurately reproduced the breakage of capsules with different mechanical properties. The proposed model can be used as a tool for designing capsules that, under stress, break and release their internal liquid at a speciﬁc time.


Introduction
A capsule is a solid material enclosing an active substance, often in liquid form. The internal structure of the capsules can vary significantly depending on the encapsulation method and the materials used [1][2][3]. It can have a simple single-core shell [4], a multicore shell, or a porous structure where the liquid is trapped [5]. Often, the purpose of the capsule is to ensure a controlled release of the liquid inside. With this aim, they are widely used in a variety of fields [6], including cosmetics [7], food [8], textiles [9], pharmacy [10], agriculture [11], and, more recently, for self-healing materials such as self-repairing concrete [12] or asphalt [13,14].
Mechanical characterisation is essential to the design of capsules for a timely release of an active substance [15]. It provides information about material parameters such as the Young's modulus and the stress-strain relationship that are key factors in designing capsules fit for their intended application. Currently, capsules are designed for end-use applications by trial and error. Computer simulations have the potential to significantly accelerate their design by providing a method for predicting the behaviour of the capsule in advance.
Conventional modelling techniques like the Finite Element Method (FEM) can be used to calculate stress in shells [16], but they are less effective when the simulation includes the fracture of the shell. Modelling fluid-solid interactions, occurring after the fracture   Here, only a very general introduction to SPH is given; for additional details, see [27]. The general idea of SPH is to approximate a partial differential equation over a group of movable computational particles that are not collocated over a grid or a mesh [28]. Each particle has an associated mass, momentum, and energy and the motion of the particle is calculated by integrating Newton's second law. Particle properties such as density, velocity, energy per unit of mass, pressure, or heat flux are obtained by interpolating neighbouring particles by means of: where f (r) is a generic property defined over the volume, r is the position vector where the property is measured, W is the so-called Kernel function, and h is a smoothing length [27]. We used the Lucy kernel [29]. The Navier-Stokes equation, for instance, can be discretised over a series of computational particles, obtaining: where P is the pressure, t represents time, v is the velocity, m is the mass, ρ is the density associated with particles i and j, and F E refers to the external forces acting on the fluid (e.g., gravity). The viscosity term ( ) used in this study is defined in [30]. Equation (2) requires an equation of state that relates density and pressure. In this work, we used Tait's equation [31]: where c 0 is a reference speed of sound and ρ 0 is a reference density. To ensure weak compressibility, c 0 was chosen to be at least 10 times larger than the highest fluid velocity. Π ij in Equation (2) introduces the viscosity tensor; here, we used the so-called artificial viscosity [30], defined as: where v ij = v j − v i and ρ ij = ρ j + ρ i . The dimensionless parameters α and b ensure the stability of the simulation. In our case, b = 0.01, while: where ν is the kinematic viscosity. At each time step, the local density is updated according to the partition of unity equation [30]:

Lattice Spring Model (LSM)
The LSM was used to model the shell of the capsule. Here, only a very general introduction to LSM is given; for additional details see [32]. The computational particles that constitute the solid shell are linked together by spring-like bonds and their trajectories calculated with the Newtonian equations of motion: Processes 2021, 9, 354 4 of 13 where U tot is the potential used to model the bonds between interconnected particles. In this study, linear spring was used, which corresponds to the harmonic potential: where r 0 is the equilibrium distance and k b is the Hookean constant. The Young's modulus of the materials relates to k b , as explained in Section 2.3. Equation (8) represents the forces that maintain two particles at a specific distance (r 0 ). The distribution of the shell's computational particles is determined by defining a quadratic mesh (see Figure 1) over an empty sphere with external diameter 3.1 × 10 −5 m and an internal diameter of 2.6 × 10 −5 m. Nearest neighbour particles and next-nearest neighbour particles are linked with computational springs. Breakage is simulated assuming that, if the distance between two particles exceeds a maximal value (r max ), the bond is broken and the two particles separated (Figure 2).  How kb and rmax are calibrated to reproduce the capsule mechanical properties is discussed later on.

Coupling the Two Models
Besides the interaction between two liquid particles (model with SPH) and two solid particles (model with LSM), the model also requires specifying the interaction between liquid and solid particles that occurs at the solid-liquid interface. Here, only a general introduction to coupling SPH and LSM is given; for more details, see [33]. In continuum mechanics, three conditions must be satisfied: σs•n = σf•(−n) (continuity of stresses condition) (11) where n is the vector normal to the boundary, u is the displacement of the solid, v is the velocity of the liquid, and σs is the stress tensor in the solid and σf in the fluid. The advantage of using a particle method is that we only solve the Newton equation of motion for every particle. The only difference is in the force acting on the particle. Equations (9)-(11) are no exceptions. They do not need a special treatment, but only forces that replicate the effect of the boundary conditions. For Equation (9), this is achieved by using the repulsive potential of the type: How k b and r max are calibrated to reproduce the capsule mechanical properties is discussed later on.

Coupling the Two Models
Besides the interaction between two liquid particles (model with SPH) and two solid particles (model with LSM), the model also requires specifying the interaction between liquid and solid particles that occurs at the solid-liquid interface. Here, only a general introduction to coupling SPH and LSM is given; for more details, see [33]. In continuum mechanics, three conditions must be satisfied: σ s ·n = σ f ·(−n) (continuity of stresses condition) (11) where n is the vector normal to the boundary, u is the displacement of the solid, v is the velocity of the liquid, and σ s is the stress tensor in the solid and σ f in the fluid. The advantage of using a particle method is that we only solve the Newton equation of motion for every particle. The only difference is in the force acting on the particle. Equations (9)-(11) are no exceptions. They do not need a special treatment, but only forces that replicate the effect of the boundary conditions. For Equation (9), this is achieved by using the repulsive potential of the type: where the constants , σ, a, and b are given in Table 1. This repulsion force does not influence the simulation itself. It is only used to avoid particle compenetrating and its parameters are determined to achieve this goal without requiring a smaller time step. The cut-off r < σ in Equation (12) ensures that, when the particles are not in contact, the force is zero. The parameters in Equation (9) were chosen by trial and error to guarantee no penetration without decreasing the timestep. If the repulsion potential is too weak, fluid will penetrate inside the solid; if it is too strong, the simulation will require small timesteps to remain stable.
No-slip conditions are modelled by superimposing virtual fluid particles above the solid particles at the interface, providing zero relative velocity at the solid-liquid boundary for the liquid particles that interact with the interface. Finally, in particle methods, the continuity-of-stress is automatically satisfied (it is a consequence of the law of actionreaction) and does not require to be specificity implemented.
The calculations were carried out with LAMMPS, a molecular dynamics software package that also accounts for other particle methods, such as SPH. The reader is referred to [34,35] for additional details on the numerical formulations and schemes used in the simulation of, respectively, the structure and the fluid.
The model does not account for temperature changes. However, this can be easily implemented in DMP [18].

Geometry
The starting geometry in the simulations was a three-dimensional single liquid core capsule with a diameter of 0.031 mm. As the purpose of this simulation was to simulate compression of single particles between two parallel plates, the drivers were defined as rectangular geometries that apply force to the capsule. The parallel plates were modelled with solid particles (Figure 3). The lower plate was stationary, while the upper plate moved as a rigid body under the effect of an applied force that compressed the capsules during the compression test. Processes 2021, 9, x FOR PEER REVIEW 6 of 13 Figure 3. Initial simulation geometry. Figure 3 shows the geometry of the shell and the parallel plates at the beginning of the simulation. A transversal cut of the capsule shows the fluid particles inside the capsule. In order to prevent penetration, repulsive forces were used between the fluid and the shell, the fluid and the plates, and the shell and the plates by applying no-penetration and no-slip boundary conditions, as explained in the previous section. The shell was built with three particle layers. Simulations with additional layers (i.e., more particles) were carried out and show similar results. In this study, we only considered elastic capsules. Viscoelastic can be modelled by implementing the method proposed in [36]. From Figure  3, it can be appreciated that the shell was not completely filled with fluid. The empty space was created by the small deformation of the capsule due to gravity and the resulting settling of the fluid particles. In theory, it is possible to add more particles to replenish the capsule. However, we decided not to do this, since the same phenomenon can also occur in real particles. This space is normally within 1%-2% of the total volume. We verified that within this range, and the results were not affected.

Main Simulation Parameters
The parameters used in the simulations can be divided into three groups (i.e., SPH, MSM, and boundaries) and are shown in Table 1.
The parameter that determines the elastic behaviour of the model is kb, the Hookean coefficient of the computational springs. The actual Young's modulus of the capsules can be calculated from compression experiments with the technique proposed by Yap et al. [37]. In our case, it corresponds to E = 4.85 MPa. Once the Young's modulus is known, the value of kb can be calculated with the following equation [38]: where E is the Young's modulus, t0 is the thickness of the shell of the capsule, and n is the number of layers of the shell. The timestep used in the simulation was Δt = 1 × 10 −10 s. We carried out a sensitivity analysis of the results with the timestep to make sure the results were independent of the timestep adopted. The parallel plates applied a progressive pressure to the capsule until the capsule broke, replicating the experimental conditions. From the simulations, therefore, we can plot the force-displacement data and, for validation purposes, compare the plot with that obtained from the experiments.

Parametric Study and Influence of kb and rmax
The capsule's deformability and strength are governed by two parameters: the Hookean coefficient of the computational springs (kb) and the maximum distance between solid particles after which the computational spring breaks (rmax). In order to analyse the influence of these parameters on the deformability, strength, and type of fracture of the capsule, the main simulation, fitted to the experimental force-displacement curve, was  Figure 3 shows the geometry of the shell and the parallel plates at the beginning of the simulation. A transversal cut of the capsule shows the fluid particles inside the capsule. In order to prevent penetration, repulsive forces were used between the fluid and the shell, the fluid and the plates, and the shell and the plates by applying no-penetration and no-slip boundary conditions, as explained in the previous section. The shell was built with three particle layers. Simulations with additional layers (i.e., more particles) were carried out and show similar results. In this study, we only considered elastic capsules. Viscoelastic can be modelled by implementing the method proposed in [36]. From Figure 3, it can be appreciated that the shell was not completely filled with fluid. The empty space was created by the small deformation of the capsule due to gravity and the resulting settling of the fluid particles. In theory, it is possible to add more particles to replenish the capsule. However, we decided not to do this, since the same phenomenon can also occur in real particles. This space is normally within 1%-2% of the total volume. We verified that within this range, and the results were not affected.

Main Simulation Parameters
The parameters used in the simulations can be divided into three groups (i.e., SPH, MSM, and boundaries) and are shown in Table 1.
The parameter that determines the elastic behaviour of the model is k b , the Hookean coefficient of the computational springs. The actual Young's modulus of the capsules can be calculated from compression experiments with the technique proposed by Yap et al. [37]. In our case, it corresponds to E = 4.85 MPa. Once the Young's modulus is known, the value of k b can be calculated with the following equation [38]: where E is the Young's modulus, t 0 is the thickness of the shell of the capsule, and n is the number of layers of the shell. The timestep used in the simulation was ∆t = 1 × 10 −10 s. We carried out a sensitivity analysis of the results with the timestep to make sure the results were independent of the timestep adopted. The parallel plates applied a progressive pressure to the capsule until the capsule broke, replicating the experimental conditions. From the simulations, therefore, we can plot the force-displacement data and, for validation purposes, compare the plot with that obtained from the experiments.

Parametric Study and Influence of k b and r max
The capsule's deformability and strength are governed by two parameters: the Hookean coefficient of the computational springs (k b ) and the maximum distance between solid particles after which the computational spring breaks (r max ). In order to analyse the influence of these parameters on the deformability, strength, and type of fracture of the capsule, the main simulation, fitted to the experimental force-displacement curve, was

Validation of the Simulation of Capsule Compression between Two Parallel Plates
A parallel compression test of a single core-shell capsule under compressive load was reproduced in the simulations. The force-displacement curve obtained from these simulations was compared to the corresponding data obtained experimentally [25] for a core-shell with the same size and internal structure. For validation purposes, we compared the results with only one type of capsule available in the literature [25]. However, it must be noted that the parameters k b describing the mechanical properties of the capsule are calculated by first-principles Equation (13), and are is not the result of fitting a free parameter to the experimental data.
In order to determine the optimal resolution of the core-shell model, various simulations with an increasing number of computational particles were compared. Each simulation had a different layer used to model the shell. The best compromise between accuracy and computational times was found with three layers. Figure 4 compares the force-displacement curve obtained with the model and with the experiments. Furthermore, an additional simulation with the capsule without fluid inside was carried out in order to understand the influence of the fluid on the capsule s strength.
The results show that the compression can be divided into three regions: (i) an elastic region where the force is proportional to δ 3/2 (where δ is the displacement), as predicted by Hertzian theory, and all of the load is supported by the external shell; (ii) an inelastic region where the stress reaches the fluid that adds additional resistance to the compression; (iii) a final region where the capsule breaks and collapses.
These three different stages were observed in both the experimental and the DMP force-displacement curves. The DMP results are very similar to the experimental results and reached the same maximal force required to break the capsule. The force-displacement curve from the simulation of the shell without fluid shows a completely different shape, which indicates that the internal fluid affects significantly the deformation of the capsule, especially at high deformations. It is important to highlight that the model does not account for plastic deformation. Locally, the rupture of every single bond was "fragile." However, this still produced a gradual rupture at the global level. In Figure 4, for instance, we defined the region before the rupture as "inelastic" because we did not observe capsule break-up. However, some of the bonds in the solid structure began to rupture, weakening the structure and producing the observed inelastic behaviour. Figure 5 shows that the DMP model is consistent with the fracture and the capsule deformation for different coaxial load values applied during the compression test of the real capsule tested in the laboratory [25]. In order to save computational time, the resulting speed of the driver was higher in the simulations shown in Figure 5. However, quasi-static conditions were verified by carrying out simulations at different driver speeds, enabling a comparison between the force-displacement curves.

Influence of the Fluid on the Capsule's Strength
Taking into consideration the previous results, it can be said that the internal liquid has an important effect on the force-displacement relationship. Figure 6 shows an image sequence of the simulation of the parallel compression test of a single liquid core-shell capsule. Figure 6a shows the pressure of the internal fluid at the beginning of the test. As the capsule was compressed by the upper plate load and deformed, the pressure increased due to the load applied by the upper plate (Figure 6b,c). Finally, when the capsule started to break, the pressure of the fluid inside the capsule reduced because of the fluid released through the broken shell (Figure 6d).

Influence of the Fluid on the Capsule's Strength
Taking into consideration the previous results, it can be said that the internal liquid has an important effect on the force-displacement relationship. Figure 6 shows an image sequence of the simulation of the parallel compression test of a single liquid core-shell capsule. Figure 6a shows the pressure of the internal fluid at the beginning of the test. As the capsule was compressed by the upper plate load and deformed, the pressure increased due to the load applied by the upper plate (Figure 6b,c). Finally, when the capsule started to break, the pressure of the fluid inside the capsule reduced because of the fluid released through the broken shell (Figure 6d).

Deformability and Strength of the Capsule (a Parametric Study)
In this section, a parametric study was carried out to understand how the behaviour of the capsule is affected by its mechanical properties: in particular, by kb (which is related to the Young's modulus by Equation (13) and rmax (which provides the maximal local displacement before rupture). We studied three values of kb (i.e., 14, 4.7, and 1.6 N m −1 ) and rmax (i.e., 1.03•r0, 1.2•r0, and 1.27•r0), and all nine possible pairings of these parameters.
Lower rmax values imply that the distances required to break the computational springs are shorter. This can be seen, for instance, in Figure 7b. Given the same kb, the capsule broke when the displacement reached 45% for rmax = 1.03•r0 and 54% for rmax = 1.2•r0. The load that initiated the breakage of the capsule was also lower (0.54 versus 1.46 mN). This effect was observed for all three values of kb, but it was more evident at high values of kb.

Deformability and Strength of the Capsule (a Parametric Study)
In this section, a parametric study was carried out to understand how the behaviour of the capsule is affected by its mechanical properties: in particular, by k b (which is related to the Young's modulus by Equation (13) and r max (which provides the maximal local displacement before rupture). We studied three values of k b (i.e., 14, 4.7, and 1.6 N m −1 ) and r max (i.e., 1.03·r 0 , 1.2·r 0 , and 1.27·r 0 ), and all nine possible pairings of these parameters.
Lower r max values imply that the distances required to break the computational springs are shorter. This can be seen, for instance, in Figure 7b. Given the same k b , the capsule broke when the displacement reached 45% for r max = 1.03·r 0 and 54% for r max = 1.2·r 0 . The load that initiated the breakage of the capsule was also lower (0.54 versus 1.46 mN). This effect was observed for all three values of k b , but it was more evident at high values of k b . As Figures 7 and 8 show, lower values of rmax mean that the local deformations required to break a computational bond are lower. As a result of this, breakage occurred at lower displacements in the simulations where rmax was lower. High values of kb mean high Young's moduli (see Equation (13)). The shell was more rigid and, when it broke, the fracture started at a single location and propagated from there, forming a single wedgeshaped gash. When kb was low, the capsule was flexible and showed little resistance to compression. In this case, the capsule tended to break at various locations at the same time during compression. As Figures 7 and 8 show, lower values of r max mean that the local deformations required to break a computational bond are lower. As a result of this, breakage occurred at lower displacements in the simulations where r max was lower. High values of k b mean high Young's moduli (see Equation (13)). The shell was more rigid and, when it broke, the fracture started at a single location and propagated from there, forming a single wedgeshaped gash. When k b was low, the capsule was flexible and showed little resistance to compression. In this case, the capsule tended to break at various locations at the same time during compression. Processes 2021, 9, x FOR PEER REVIEW 11 of 13

Conclusions
This study proposed a Discrete Multi-Physics model to simulate capsules formed by a shell containing a fluid. The liquid core was implemented with Smoothed Particle Hydrodynamics and the solid shell with the Lattice Spring Model.
The model was validated against experimental compression of a single core-shell capsule between two parallel plates. The results showed that the model replicates with good accuracy the behaviour of the real capsule. Three different stages were highlighted in the force-displacement curve occurring during the compression. In the first stage, the load was supported by the external shell and the behaviour was Hertzian (i.e., the force was proportional to δ 3/2 ). During the second stage, the stress reached the fluid, which added additional resistance to the compression. Finally, in the third stage, the capsule broke and collapsed.
After the model had been validated, we carried out a parametric study that showed how the breakage of the capsule depended on its material properties. By changing the material properties, we varied (i) the maximal stress the capsule can resist before breaking, and (ii) the mode of breakage. In some cases, the shell cracked at one point that opened a wide gash in the capsule. In other cases, the shell fractured at multiple locations. This affected the way the liquid was released from the capsule. In the first scenario, the release was not uniform and was concentrated at one side of the capsule. In the second, the liquid was released from multiple points distributed around the capsule.
The proposed model can help the design of capsules with specific features. For instance, we can adjust the mechanical properties of the capsule to make sure that it breaks only when a precise stress level is reached. The model also simulates the liquid release occurring after breakage. In this way, we can design the capsule to make sure that the liquid is released effectively. In fact, the capsule breaks in different ways according to its mechanical properties, and only specific combinations of these properties allow for a uniform release of the liquid after breakage.

Conclusions
This study proposed a Discrete Multi-Physics model to simulate capsules formed by a shell containing a fluid. The liquid core was implemented with Smoothed Particle Hydrodynamics and the solid shell with the Lattice Spring Model.
The model was validated against experimental compression of a single core-shell capsule between two parallel plates. The results showed that the model replicates with good accuracy the behaviour of the real capsule. Three different stages were highlighted in the force-displacement curve occurring during the compression. In the first stage, the load was supported by the external shell and the behaviour was Hertzian (i.e., the force was proportional to δ 3/2 ). During the second stage, the stress reached the fluid, which added additional resistance to the compression. Finally, in the third stage, the capsule broke and collapsed.
After the model had been validated, we carried out a parametric study that showed how the breakage of the capsule depended on its material properties. By changing the material properties, we varied (i) the maximal stress the capsule can resist before breaking, and (ii) the mode of breakage. In some cases, the shell cracked at one point that opened a wide gash in the capsule. In other cases, the shell fractured at multiple locations. This affected the way the liquid was released from the capsule. In the first scenario, the release was not uniform and was concentrated at one side of the capsule. In the second, the liquid was released from multiple points distributed around the capsule.
The proposed model can help the design of capsules with specific features. For instance, we can adjust the mechanical properties of the capsule to make sure that it breaks only when a precise stress level is reached. The model also simulates the liquid release occurring after breakage. In this way, we can design the capsule to make sure that the liquid is released effectively. In fact, the capsule breaks in different ways according to its mechanical properties, and only specific combinations of these properties allow for a uniform release of the liquid after breakage.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.