1 Introduction

Solid tablets are commonly used in various applications such as medicine, agriculture, or water purification. In the pharmaceutical industry, solid tablets are one of the most popular dosage forms. The reason behind this is that tablets are relatively easy to manufacture, store and transport, and provide a convenient way to administer drugs [1]. The majority of the pharmaceutical tablets are coated as a part of the manufacturing process. The coating serves several functions such as odour and taste masking, increasing swallowability, physical and chemical protection during storage and transportation, increase aesthetic appeal and control drug release profile [2, 3].

Pan coating is one of the first and most used methods for applying a film coating onto a tablet. The process consists of a batch of tablets deposited onto a rotating bed. The spray covers the tablets with a coating solution. Heated air evaporates the solvent leaving a solid film on the surface of the tablets [4, 5]. Residual moisture in the porous tablet can cause swelling, which affects coating morphology and may accelerate drug degradation and affect the dissolution behaviour [6]. Hence, understanding and predicting the spreading, absorption and evaporation mechanisms are crucial for reducing moisture content in the tablet.

The three different mechanisms of water transport and removal occur at different time scales. A liquid droplet impacting on a porous surface induces a spreading on the surface and absorption inside the porous material due to capillary forces. In particular, droplet spreading is a fast process of the order of milliseconds while evaporation takes place in the range between a few seconds and minutes [7]. We will focus our attention on the dynamics of the first phases: spreading (or kinematic) phase and absorption (or capillary) phase. The dynamics of spreading and absorption depend on different physical parameters. The spreading rate depends on the properties of the liquid (density, viscosity and surface tension), initial droplet geometry and impact velocity, wettability (contact angle) and surface roughness [8]. Absorption is controlled by the properties of the liquid and the porous substrate (porosity, pore size, wettability) [9]. During the process of coating, a droplet can rebound, splash, or simply deposit when impacting the tablet, depending on the velocity, size and viscosity of the atomized droplet [10]. In the present analysis, we will focus on the simple deposition of droplets, when the impact velocity is low and the spreading and absorption dynamics are strongly related to the tablet surface roughness, and, more generally, to its microstructure [11, 12].

The liquid droplet impact on porous surfaces has been studied previously for several different applications. For instance, inkjet printing and spray coating of paper [13,14,15], spray cooling [16], environmental applications [17] and spray coating of tablets [7, 10]. These and other studies have tried to characterise the behaviour of liquid droplets upon the impact with porous media. Park et al. [18] formulated a mathematical theory to estimate the maximum spreading rate at low impact velocity. Attané et al. [19] developed a simplified 1-D model based on the energy equation. A collection of experiments and numerical simulations to evaluate the impact and spreading of a pure liquid water droplet on pharmaceutical tablets can be found in the work by Shaari [20]. A Volume-Of-Fluid (VOF) approach was employed to predict the spreading, the model can however not be applied for the absorption phase.

The absorption of deposited droplets in porous media has been rarely investigated given the extreme experimental or numerical resolutions required to detect water down to the pore scale. Given the previous limitation, experimental absorption studies have measured the volume of the depleting droplet [9, 13, 21]. In the previous studies, the dynamics of the water front inside the porous media were assumed to follow a classic Washburn equation. However, experiments and numerical simulations have shown deviation from Washburn’s law under different conditions [22,23,24] and direct observations of the liquid content redistribution in the porous medium are required. Various experimental techniques have been recently applied to measure the liquid mass inside porous substrate such as X-ray radiography [25], magnetic resonance imaging (MRI) [26] and neutron radiography [27]. The latter study characterises the droplet dynamics inside natural stones from spreading to evaporation in a time-resolved manner. Their results show a classic Washburn behaviour for the water content during the absorption phase.

Given the difficulties monitoring droplet absorption in porous media experimentally in terms of cost and accuracy, numerical simulations can play an important role in the understanding of such a complex phenomenon. Several numerical methodologies have been used for simulating flows into porous substrates, such as volume averaged techniques [28], the immersed boundary and Lattice-Boltzmann methods [29, 30]. The former are based on a continuum approach where the pore-scale quantities are represented via effective macroscopic parameters. Therefore, the prediction of the fluid dynamics is limited by the formulation and accuracy of these effective parameters. On the other hand, immersed boundary and Lattice-Boltzmann methodologies aim at solving the pore-scale dynamics and provide a more comprehensive and accurate representation of the flow behaviour inside the porous microstructure. In particular, the Lattice-Boltzmann method offers several advantages, such as intrinsic parallelism and easy implementation of boundary conditions [31]. It also allows to efficiently simulate multiphase flows by taking into account surface tension effects. In this work, we use a Lattice-Boltzmann single-component multiphase model based on the Shan-and-Chen formulation [32] in order to solve the two-phase flow that characterises droplet spreading and absorption into a porous substrate.

Our numerical approach is briefly discussed in Sect. 2. In Sect. 3, we will describe the geometry and properties of the porous medium and the three different cases under investigation. Sect. 4 contains a description of the results extracted from the numerical simulations. Finally, Sect. 5 presents a discussion and conclusions.

2 Methodology

We perform numerical simulations by means of the Lattice-Boltzmann method (LBM). The LBM solves the Boltzmann transport equation in discretised phase space, that is, in a regular cartesian lattice grid. The adoption of a regular grid for solving the flow provides some computational advantages, such as the ease of handling complex geometries and of parallelisation [31]. In particular, we use a 3-dimensional lattice grid defined by 18 nodes on the convex shell and a centroid (D3Q19). The Lattice-Boltzmann multiphase algorithm has been previously used and validated in complex microstructures [33, 34]. The Lattice-Boltzmann equation reads as:

$$\begin{aligned} f_r (\varvec{x} + \varvec{c}_r \delta t, t+\delta t) - f_r(\varvec{x},t) = - \tau ^{-1} ( f_r(\varvec{x},t)-f^{eq}_r (\varvec{x},t) ) + F_r , \end{aligned}$$
(1)

where \(f_r(\varvec{x},t)\) is the distribution function at position \(\varvec{x}\) and time t along the r-th direction of the lattice grid; \(\varvec{c}_r\) is the so-called discrete velocity vector along the r-th direction, \(F_r\) is the body force and \(f_r^{eq}(\varvec{x},t)\) is the equilibrium distribution function. The relaxation time \(\tau\) is a function of the kinematic viscosity \(\nu\) and the speed of sound \(c_s\). The body force is proportional to the gravitational acceleration g and is defined following the scheme proposed by Guo et al. [35]:

$$\begin{aligned} F_r = \Bigg ( 1-\frac{1}{2\tau } \Bigg ) w_r \Bigg ( \frac{\varvec{c}_r-\varvec{u}}{c_s^2} +\frac{\varvec{c}_r\cdot \varvec{u}}{c_s^4}\varvec{c}_r \Bigg ) (\rho g) \end{aligned}$$
(2)

with \(w_r\) representing the weighting coefficients [31] and \(\rho\), \(\varvec{u}\) the density and velocity vector. Within the Lattice-Boltzmann framework, interfacial effects in multiphase flow, such as surface tension, are formulated on the mesoscopic molecular level. In particular, we adopt the Schan-and-Chen scheme that introduce a density-dependent Van-der-Waals-like interaction force [32] mimicking gas-liquid surface tension \(\gamma\). This force is formulated as:

$$\begin{aligned} F_{sc} = - G\varPsi (\mathbf{x },t) \sum _r w_r \varPsi (\varvec{x} + \varvec{c}_r \delta t,t)\varvec{c}_r \end{aligned}$$
(3)

with G the interaction strength between the phases (\(G=-5.5\) in the present simulations) and \(\varPsi (\rho )= 1-e^{-1/\rho }\) is the density-dependent pseudo-potential function. The interaction force between the phases yields a non-ideal equation of state \(P(\rho )\) of the form:

$$\begin{aligned} P(\rho )=\rho c_s^2 + \frac{G}{2}c_s^2\ {\varPsi (\rho )}^2 \ . \end{aligned}$$
(4)

The desired equilibrium contact angle is imposed through spatial averaging of the density-dependent potential function, at the gas-liquid-solid contact point [36]. The pseudo-potential at the solid wall \(\varPsi _{wall}\) is then given by:

$$\begin{aligned} \varPsi _{wall}(\rho ) = N^{-1} \sum _N \varPsi (\rho ) + \varDelta _w, \end{aligned}$$
(5)

where N indicates the nearest fluid computational nodes at the three-phase contact line and \(\varDelta _w\) is a numerical parameter that determines the contact angle. For sake of clarity, \(\varDelta _w\) is a measure of the affinity between liquid and solid; a positive value indicates a hydrophilic solid phase while a negative value a more hydrophobic one.

The macroscopic flow quantities density and velocity \((\rho ,\varvec{u})\) are extracted from the formulation of the hydrodynamic moments [35]:

$$\begin{aligned} \begin{aligned} \rho&= \sum _r f_r \\ \rho \varvec{u}&=\sum _r \varvec{c}_r f_r + 1/2\ \rho \varvec{g} + 1/2\ F_{sc} \ , \end{aligned} \end{aligned}$$
(6)

and the equilibrium distribution functions are given by:

$$\begin{aligned} \begin{aligned} f^{eq}_r&= w_r \rho \ \Bigg ( 1 - \frac{\varvec{u}_{eq}\cdot \varvec{u}_{eq}}{2c_s^2} \Bigg ) \ ,\ r=1\\ f^{eq}_r&= w_r \rho \ \Bigg (1+\frac{\varvec{c}_r \cdot \varvec{u}_{eq}}{c_s^2} + \frac{(\varvec{c}_r \cdot \varvec{u}_{eq})^2}{2c_s^4} - \frac{\varvec{u}_{eq}\cdot \varvec{u}_{eq}}{2c_s^2} \Bigg ) \ , \ r=2-19 \end{aligned} \end{aligned}$$
(7)

where the equilibrium velocity is defined as \(\varvec{u}_{eq}=\varvec{u} + (\tau - 1/2 ) F_{sc}\) [32].

3 System description

We investigate the dynamic spreading and absorption of a liquid droplet of the size \(d_0\), density \(\rho\) and dynamic viscosity \(\mu\), into a porous medium composed of mono-disperse randomly overlapping spherical particles of size \(d_s\). At time zero, the droplet is placed on the top of the porous medium in a box of height \(4\ d_s\) filled with gas. The gravitational force acts along the vertical direction z and periodic and free-slip boundary conditions are imposed along x, y and z, respectively. A typical configuration of the system is shown in Fig. 1 with the liquid droplet (blue) being absorbed inside the porous tablet.

Fig. 1
figure 1

Typical configuration of the system where the liquid droplet (blue) is being absorbed on the top of the porous medium. (Color figure online)

Full size image

The porous medium height along the direction of droplet absorption z is \(L_z=5.6\ d_s\) and the horizontal dimensions are \(L_x=L_y=7.5 \ d_s\). The porosity is \(\epsilon =0.25\) and the mean pore size \(d_p=0.28\ d_s\). The pore size distribution is shown in Fig. 2. Values of porosity in real pharmaceutical tablets are lower (0.05-0.1) but simulations of such systems would require high computational costs. We found the value of 0.25 to be a good compromise between computational speed and detailed description of the problem. We assume that the droplet impact velocity is low and that the corresponding kinetic energy can be considered negligible. In such cases, the dynamic behaviour of the droplet deposition and absorption is primarily governed by the balance between gravitational and capillary forces, as well as the viscous forces acting inside the pores, and it is strongly affected by the porous microstructure.

Fig. 2
figure 2

Pore diameter distribution of the porous domain

Full size image

The system under investigation is governed by several dimensionless numbers. We define the Bond number and the ratio between the viscosities of the two-phase flow as:

$$\begin{aligned} \begin{aligned} {\textit{Bo}}&=\frac{\rho g d_0^2}{\gamma } \ ,\\ M&=\frac{\mu }{\mu _g} \ , \end{aligned} \end{aligned}$$
(8)

where \(\mu _g\) is the gas dynamic viscosity. The sizes of particles that compose typical pharmaceutical tablets can vary significantly on the basis of the specific technological manufacturing process. The particle size can vary as \(\sim O(10^{-2}-1)\)\({\mathrm {mm}}\) [37, 38]. The droplet-to-particle diameter ratio can be estimated as \(d_0/d_s\sim O(0.1-1000)\) [10]. The coating of pharmaceutical tablets is therefore characterised by the deposition and absorption of droplets whose diameters are comparable to the diameter of the particles composing the porous substrate. If we consider the rheological characteristics of the droplet similar to the ones of water (\(\rho \sim 10^3\ {\mathrm {kg}}/{\mathrm {m}}^3\), \(\mu \sim 10^{-3}\ {\mathrm {kg}}/{\mathrm {m}}/{\mathrm {s}}\), and \(\gamma \sim 0.07\)\({\mathrm {N}}/{\mathrm {m}}\)), the Bond number results in the order of \({\textit{Bo}}\sim 10^{-3}\) for small droplets and \({\textit{Bo}}\sim 10\) for large droplets. Furthermore, the dynamic viscosity ratio (water-air) results in \(M\sim 50\). The droplet dynamic depends on the equilibrium contact angle \(\theta\), a value that depends on the material properties of the particles. A non-dimensional number characterising a flow in porous media is the Darcy number \(Da=k/d_0^2\), where k is the permeability of the medium. We assume a Newtonian behaviour for the droplet while in reality its rheology is more complex since dissolved polymers and undissolved solid particles are present in the coating solution [2].

We performed three numerical simulations with varying the Bond number as \({\textit{Bo}} = 1.5\) and 3.5, the contact angle as \(\theta =65^\circ\) and \(\theta =80^\circ\), and the droplet to particle diameter ratio as \(d_0/d_s=2\) and 3. These values in the three considered cases are summarised in Table 1. It should be noted that such values of the Bond number are of particular interest since they represent situations where gravitational and capillary forces are of the same order of magnitude and both contribute to the droplet dynamics. The viscosity ratio is set as \(M\sim 20\), which is a good compromise between accuracy in representing a water-air system and stability of the numerical algorithm. The domain discretisation consists of \(300\times 300\) computational cells in the horizontal directions and 200 in the vertical direction.

Table 1 Non-dimensional parameters of the three simulated cases. The viscosity ratio is set as \(M\sim 20\)
Full size table

4 Results

A qualitative scenario of the spreading and absorption processes for the three different cases is shown in Fig. 3. The plots depict the water content in one middle plane section in the vertical plane. The porous domain is located on the lower parts of each panel (negative ordinate), and the solid part has been removed to identify the water inside the tablet better. Cases A and B share similar dynamics from a qualitative point of view, whereas Case C shows different dynamics. Spreading and absorption phases are characterised by two different time scales. The spreading phase is fast, whereas significantly slower dynamics govern the absorption. The two phases are classically separated by a non-dimensional time \(t^*=t/t_c\) defined as the ratio between time and the so-called capillary spreading time scale [39]:

$$\begin{aligned} t_c=\sqrt{\frac{\rho d_0^3}{\gamma }}. \end{aligned}$$
(9)

In our simulations, the capillary time scale ranges from 22 (Case C) to 42 ms (Cases A-B). All plots in a horizontal line of Fig. 1 are taken at the same physical time and, consequently, their corresponding capillary spreading time \(t^*\) may differ. Figure 3 shows that after a non-dimensional time \(t^*\approx 1\), the shape of the water front changes slowly in time. The most hydrophilic configuration (Case B) exhibits a larger water penetration in the porous medium. Case C differs from the others since water absorption mainly proceeds in the vertical direction compared to the horizontal plane. We will show that these dynamics will affect the absorption statistics and are linked to the fact that the water content of a small droplet can be contained in fewer pores compared to larger droplets.

Fig. 3
figure 3

Instantaneous snapshots of the shape of the water content after the impact on the porous surface: case A (left panels), case B (middle panels), case C (right panels)

Full size image

To quantitatively describe the spreading and absorption dynamics, we will introduce and measure the time behaviour of several characteristic geometrical lengths shown in Fig. 4. The left panel shows the initial droplet shape before the impact with the initial diameter of \(d_0\). The right panel describes the relevant geometrical scales of spreading and absorption. Spreading can be described by measuring the equivalent wetted area diameter \(d_{CL}\) and the maximum droplet height from the surface h. The wetted area diameter \(d_{CL}\) is obtained as the diameter of the circle of the equivalent wetted area at the tablet surface. Absorption can be related to the wetting front position \(h_p\) measured as the maximum water penetration depth in the porous medium.

Fig. 4
figure 4

Visual representation of the characteristic geometrical lengths of spreading and absorption. Left panel: droplet before the impact. Right panel: deformed droplet after the impact

Full size image

Droplet spreading is characterised by the time evolution of the wetted area diameter and the maximum height from the surface which can be normalised with the initial droplet diameter. The two non-dimensional quantities are called spreading ratio (\(\beta =d_{CL}/d_0\)) and normalised maximum height (\(H=h/d_0\)) and are plotted in Fig. 5 in logarithmic scales. The plots highlight the role of the capillary spreading time \(t_c\) in separating spreading from absorption. In particular, the spreading phase terminates when \(\beta\) and H reach a quasi-steady value in correspondence of the capillary time scale (\(t\approx t_c\) or equivalently \(t^* \approx 1\)).

The spreading ratio can be approximated by a power law:

$$\begin{aligned} \beta =C t^{*\alpha }, \end{aligned}$$
(10)

where C and \(\alpha\) are two positive coefficients. Various experimental and numerical studies have shown that the coefficient C and exponent \(\alpha\) depend on the contact angle, fluid properties, and porosity [40,41,42]. We find a power-law exponent of \(\alpha =0.5\) for all the cases in the initial phase. The absorption phase (\(t^* >1\)) shows a constant value of \(\beta\) for all the simulations.

A lower contact angle (Case B-most hydrophilic case) results in a slightly higher maximum spreading. To test experimentally whether this deduction is correct is problematic because it is not possible to change surface wettability without changing the liquid or the porous sample, which would change other properties besides the contact angle. In contrast with our results, simulations by Reis et al. [43] found that a lower contact angle reduced the spreading ratio. Nevertheless, other Lattice-Boltzmann studies done by [41, 44], and [42] performed on simplified porous media (regular grid of cylindrical or square pores) have also shown that a lower contact angle increases the spreading rate. The disagreement between the former numerical model [43] and the Lattice-Boltzmann studies can be explained by a different model used for describing the contact angle dynamics and the surface wettability. In the former model, the contact angle dynamics are defined by means of two different values, one for the advancing contact line and one for the receding. The underlying assumption is that the dynamic contact angle does not depend on the velocity of the contact line, an assumption that is only valid for high velocities of the order of meters per second. Meanwhile, in the Lattice-Boltzmann framework, the solid-phase wettability is implemented by imposing a fictitious wall density [45] or–as we do here–a specific value of the pseudopotential function that defines the desired contact angle at equilibrium [46]. Furthermore, the Lattice-Boltzmann method makes use of a diffusive interface approach where the density between different phases changes with continuous variations. The dynamics of the contact angle and its variation from the equilibrium value are then affected by the presence of external forces, such as gravity, a condition that makes the contact angle depend on the velocity [47].

Larger droplets have a higher maximum spreading since a decrease in the droplet diameter reduces the effect of inertia and thus the relative spreading as confirmed in experiments [27] and numerical simulations [15].

The right panel of Fig. 5 shows the normalised maximum droplet height. Again, the spreading phase terminates at time scales comparable with the capillary spreading time.

Fig. 5
figure 5

Time evolution of the spreading ratio (left panel) and normalised maximum height (right panel)

Full size image

The absorption process is characterised by the amount of liquid penetrating inside the porous medium in time. The normalised maximum penetration depth \(H_{p}\) is shown in the left panel of Fig. 6. Many experiments available in the current literature, such as [27], report an absorption phase that starts after the droplet spreading terminates. On the contrary, our data show that the absorption phase starts together with the spreading at times smaller than the capillary spreading time. This contradiction arises because current experimental instruments do not have enough spatial and temporal resolution to capture small volumes of liquid inside the porous media at times of the order of few milliseconds. Moreover, the experiments have been conducted in porous materials (i.e. rocks [27]) that have a lower porosity compared to our setup. Our results suggest that the capillary time cannot be considered as the starting time of the absorption phase, and a different relevant time function of porosity should be defined to characterise this process. An extensive parametric study at different porosities is needed to determine such a characteristic time.

Figure 6 shows that the quantity \(H_{p}\) follows a classic Washburn law behaviour [48]:

$$\begin{aligned} H_{p}=k \sqrt{t}, \end{aligned}$$
(11)

a power law with the exponent of 0.5 and a constant \(k=\sqrt{\gamma r_c {\hbox {cos}} \theta /2\mu }\), called capillary coefficient, where \(r_c\) is the capillary radius. The most hydrophilic porous medium (Case B) enhances water absorption. As mentioned in the Introduction, some experiments and numerical simulations show deviations from Washburn’s law [22,23,24]. Here, we recover the classic result comparable with recent high-resolution experimental data [27]. A recent paper [49] provides criteria for applying the Washburn law based on experimental imbibition tests.

The right panel of Fig. 6 describes the time behaviour of the normalised absorbed volume defined as the ratio between the water volume inside the porous medium and the initial droplet volume. This quantity is essential since it is directly proportional to the amount of liquid mass inside the porous domain. Surprisingly, our results suggest two different power laws: one for times shorter (\(\alpha =1.5\)) and another for times larger (\(\alpha =0.5\)) than the capillary time scale. The difference in the exponents is due to the different behaviour of the absorbed water. During the first phase, the droplet spreads on the external surface of the porous medium and at the same time penetrates the porous structures in the horizontal and vertical directions. The horizontal absorption, or radial capillarity, has been quantified by Marmur [50], who found that the radially expanding liquid front grows with the square root of time. Combining the Marmur and Washburn results, we can easily derive that the absorbed volume should scale as \(V_{abs}\propto t^{1.5}\). When the droplet spreading terminates (\(t^*>1\)), and most of the vertical momentum has been converted in the horizontal plane, radial capillary stops and the water front only penetrates in the vertical direction inducing a classic Washburn scaling for the absorbed volume (or mass). The presence of the scaling law for \(t^*<1\) depends probably on the porosity of the medium and is linked with an earlier beginning of the absorption phase for times smaller than the capillary time. Case C reaches a steady configuration at times comparable with the capillary time and does not exhibit the scaling for \(t^*>1\).

Fig. 6
figure 6

Time evolution of the normalised maximum penetration depth (left panel) and normalised absorbed volume (right panel)

Full size image

The previous conclusions are confirmed by the time behaviour of the absorption area defined as \(A_{abs}=V_{abs}/h_p\) and plotted in Fig. 7. The quantity has been normalised with the frontal area of the initial droplet \(\pi d_0^2/4\) (left panel) and the wetted surface \(\pi d_{CL}^2/4\) (right panel). The plot in the left panel shows a power law with the exponent \(\alpha =1\) during the initial phase confirming a radial capillary phenomenology when the absorption rate is the same in the lateral and vertical direction. For longer times, in the absorption regime, the plot confirms that the water mainly penetrates in the vertical direction since the absorbed area remains constant in time. The small differences between Cases A and Case B are due to the difference in contact angle. The main difference is between Case C and the other cases since the absorption area decreases after the capillary time. This observation indicates that the absorption is narrower and is linked to the particular configuration of the water front, as shown in the right panel of Fig. 3.

The right panel of Fig. 7 displays the absorption area normalised with the wetted surface. During the initial phase, scaling arguments predict a constant value of this observable since both the absorption area and the square of the spreading ratio increase linearly with time. During the final phase, the normalised quantity increases to reach a constant value (Cases A and B) except for the case with the initial smaller droplet (Case C) showing a decrease in time. This behaviour can again be explained by observing the shape of the absorbed droplet in the right panel of Fig. 3. The liquid volume of the small droplet can be contained in a few pores. Therefore, the droplet penetration is relatively narrower and deeper (the absorption area decreases with time). This phenomenon is related to the ratio between the pore sizes and the initial droplet diameter. As the droplet approaches the pore size, fewer pores are required for absorption. Similar observations have been made in [24]. The authors found that, as the Darcy number increases, the spreading ratio decreases while the penetration depth increases. Also, if the Darcy number is large enough, the droplet will not stop penetrating as the resistance imposed by the substrate will not be sufficient to stop it. The Darcy number is inversely proportional to the droplet size, which means that as the droplet becomes smaller the spreading ratio decreases and the penetration depth increases. Using the Darcy number values for all cases found in Table 1, the same trend described in [24] was observed.

Fig. 7
figure 7

Time evolution of the normalised absorption area. This observable can be normalised by using the frontal area of the initial droplet (left panel) or the wetted surface (right panel)

Full size image

The left plot of Fig. 8 shows the normalised vertical (aligned with the gravity direction) droplet velocity evolution in time. The droplet velocity is normalised with the quantity \(\sqrt{gd_0}\). During the spreading phase, the droplet decelerates as its axial momentum is converted to a lateral one and small bounces can be observed. The droplet velocity becomes zero at longer times when a quasi-steady configuration has been reached.

The dynamic contact angle evolution of the spreading droplet is shown in the right panel of Fig. 8. The value of the dynamic contact angle depends on the input parameters and particularly on the equilibrium contact angle. Furthermore, the dynamic contact angle depends strongly on absorption. It decreases as absorption proceeds and in the case of total absorption reduces to zero. In this study, full absorption of the droplet was never reached. Nevertheless, in the most hydrophilic case (Case B), characterised by the largest absorption, the dynamic contact angle decreases as absorption progresses (below the value of the equilibrium contact angle).

Fig. 8
figure 8

Time evolution of the normalised vertical droplet velocity (left panel) and the dynamic contact angle (right panel)

Full size image

5 Discussion and conclusions

A series of numerical simulations using the Lattice-Boltzmann method have been performed to investigate the liquid droplet spreading and absorption on a porous surface. The three-dimensional reconstruction of the porous medium has been generated as a set of random overlapping uniform spherical particles with a porosity of 0.25. Three cases differing in the following relevant parameters were run: wettability, droplet size and the Bond number. In general, hydrophilicity enhances spreading and absorption of water inside the tablets while smaller droplets reduce the absorption. These results are consistent with previous simulations and experiments.

As in previous works, we have found that the capillary time scale \(t_c\) is the appropriate time for the spreading phase where the spreading ratio follows a power-law in time with the exponent of 0.5. The pre-factor of the power-law depends on the governing parameters (porosity, contact angle, Bond number).

The most exciting results of our study are in relation to the absorption mechanisms. We have observed at short simulation times that an initial absorption co-occurs with spreading. For this reason, the capillary time is not sufficient to separate the spreading and the absorption phases. A new time scale based on porosity must be defined. Nevertheless, the capillary time can separate the two different regimes of absorption. In particular, for longer times (\(t>t_c\)), the absorbed volume (or mass) follows the classic Washburn behaviour where the water front mainly penetrates in the vertical direction, proportionally to the square root of time. In the other case, the first phase of absorption (\(t<t_c\)) is characterised by both radial and vertical capillarity. We have shown that the radial capillary is responsible for the deviation of the scaling exponent from the Washburn law. Moreover, we were able to derive the new exponent by combining the results of Washburn for the vertical capillary and Marmur for the radial capillary absorption.

Smaller droplets penetrate narrower (smaller lateral absorption) than the larger ones. This effect is induced by the amount of liquid volume of the smaller droplet that can be contained in a limited number of pores. Therefore, the drop penetration is narrower and deeper (absorption area decreases with time). The pore size and thus the particle size have a significant impact since as the droplet approaches the pore size, fewer pores are required for absorption.

The dynamic contact angle was found to decrease as the absorption proceeds as it was expected. Eventually, it would decrease to zero if full absorption occurred.

Our results can be used as a benchmark to test new mathematical models of spreading and absorption in porous media. One-dimensional models are already available in literature [19, 51, 52]. A recent work [7] focuses on one-dimensional models for pharmaceutical tablets. Regarding the capillary absorption phase, the authors [7] simplified and solved the continuity and Navier-Stokes equations using the lubrication approximation theory. Our results indicate that it is crucial to also characterise the lateral absorption inside a tablet. Therefore, additional efforts are needed to develop new three-dimensional models.

Future numerical campaigns at higher resolutions can improve the current results to generate a large amount of data with conditions much closer to a real pharmaceutical tablet coating process. The geometry of the porous substrate could be generated by reconstructing the tablet with X-ray tomography for a more realistic representation of the medium. From a numerical point of view, novel Lattice-Boltzmann methods such as Shan-Chen multi-component multiphase (MCMP) should be considered for simulation of two-phase fluid flow. This approach would allow to independently tune the dynamic viscosity ratio and density ratio of the liquid-gas phases, and, also to independently introduce inertial effects induced by the impacting velocity of the droplet. The influence of other governing parameters such as viscosity, surface tension, porosity, pore size, particle size, surface roughness, permeability should be also explored. The physical model can be improved by including evaporation and a non-Newtonian rheology model for the droplet (to mimic the polymeric water coating). We expect that the inclusion of evaporation will not alter the qualitative picture described in this work; on the contrary, non-Newtonian effects can significantly alter the dynamics of spreading and absorption. The ultimate goal is to investigate the effect of the contact angle and droplet/particle size to obtain the best size distribution that optimises the coating process and the moisture content of the tablets.