2.1 Geometry of the flow

The problem of an impacting drop consists of a fluid drop launched with a velocity $U_{0}$ at a distance of $H_{0}$ from a rigid wall. In general, the normal impact velocity $U_{i}$ and the diameter $D_{0}$ of the drop are used in this problem to construct the dimensionless parameters, as will be presented in the next section. A scheme of the flow problem for the normal impact situation is illustrated in figure 1(a), while the drop impact on an inclined surface at a generic angle $\unicode[STIX]{x1D703}$ is sketched in figure 1(b). In the impacting drop problem, two important aspects analysed as time functions in the simulations/experiments are the spreading diameter $D(t)/D_{0}$ and the minimum bounce height $H(t)/D_{0}$ . These quantities are illustrated in figure 2.

Figure 1.

Description of the flow problem: (a) normal drop impact and (b) impact on an inclined surface.

Figure 2.

Illustration of the quantities: (a) the spreading diameter and (b) the minimum height of the drop.

In the present paper, we have fixed the drop diameter, for instance, $D_{0}=0.02$  m, based on the previous experimental work of Luu & Forterre (Reference Luu and Forterre2009) and Blackwell et al. (Reference Blackwell, Deetjen, Gaudio and Ewoldt2015).

The height of the drop from the plate is also fixed at $H_{0}=0.03$  m while the normal impact velocity is varied within the range $0.77~\text{m}~\text{s}^{-1}\leqslant U_{i}\leqslant 5.06~\text{m}~\text{s}^{-1}$ .

2.2 Physical guidelines

The model used in the present work can be seen as a generalization of the model for elasto-viscoplastic fluids of Saramito (Reference Saramito2007, Reference Saramito2009), but now accounting for thixotropic effects. There are two ways in which the model employed in the present work can be reduced to an elasto-viscoplastic model. One is to consider that quantities like viscosity and elastic modulus are still dependent on the structure parameter, but there is no delay between the current structure level and the equilibrium one. In other words, changes in stress induce instantaneous changes in structure impacting immediately the value of the mechanical properties. This first form is obtained by taking the limit where the thixotropic equilibrium time vanishes, i.e. when the time scale for structure changes is much shorter than the time scale of the considered processes. Another form of elasto-viscoplastic model is obtained in the other limit, when the material has constant mechanical properties that do not depend on the structure level. This second case is obtained by the model employed here by considering the limit where the thixotropic equilibrium time is very large, i.e. changes in structure are so slow that they do not happen in the time scales of the considered processes. Therefore, with respect to models for elasto-viscoplastic fluids, there is one phenomenon that is added, namely thixotropy, i.e. changes in the stress level have a finite time to impact the structure level. It is worth noticing that this additional phenomenon has influence on all the others since the mechanical properties of the material associated with elasticity and viscoplasticity are dependent on the structure parameter, whose changes are controlled by its degree of thixotropy. Therefore, a whole theory on how to introduce thixotropic effects is needed.

The general basis for the construction of a theory of how thixotropy affects the elasto-viscoplastic nature of the material is described in de Souza Mendes & Thompson (Reference de Souza Mendes and Thompson2012). Here, we highlight the main aspects. As commonly adopted in the literature, (e.g. Barnes Reference Barnes1997), a structure parameter $\unicode[STIX]{x1D706}\in [0,1]$ is the normalized quantity employed to describe the structural state of the material, where $\unicode[STIX]{x1D706}=0$ corresponds to the minimum structural level, and $\unicode[STIX]{x1D706}=1$ corresponds to the maximum structuring level. Making the mechanical properties functions of the structural state and applying an evolution equation for $\unicode[STIX]{x1D706}$ completes the backbone to taking into account thixotropic effects. The structural dependency of the viscosity and the elastic modulus of a viscoelastic model are defined so that the viscous and elastic effects are smaller in less structured states. The evolution equation depends on a driving potential between the current and the equilibrium structure states, where this last quantity is defined by the current stress level and by considering the flow curve as the locus of the attractor points corresponding to the equilibrium states. A breakdown process takes place only when the stress intensity surpasses the yield stress, providing the model with the main feature of a yield stress fluid.

2.3 The constitutive model

2.3.1 The total stress

The total stress is decomposed into three parts as

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D748}=-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D749}^{N}+\unicode[STIX]{x1D749}^{M},\end{eqnarray}$$

where $p$ is the pressure, $\unicode[STIX]{x1D749}^{N}=2\unicode[STIX]{x1D702}_{\infty }\unicode[STIX]{x1D63F}$ , where $\unicode[STIX]{x1D63F}$ is the rate of deformation tensor and $\unicode[STIX]{x1D702}_{\infty }$ is the viscosity related to the fully unstructured material. In addition, $\unicode[STIX]{x1D749}^{M}$ is a non-Newtonian part of the stress tensor that obeys an upper convected Maxwell equation of the following form

(2.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D749}^{M}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\unicode[STIX]{x1D749}^{M})-[(\unicode[STIX]{x1D735}\boldsymbol{u})\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{M}+\unicode[STIX]{x1D749}^{M}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}]=-\frac{G(\unicode[STIX]{x1D706})}{\unicode[STIX]{x1D702}(\unicode[STIX]{x1D706})}\unicode[STIX]{x1D749}^{M}+2G(\unicode[STIX]{x1D706})\unicode[STIX]{x1D63F},\end{eqnarray}$$

where $G(\unicode[STIX]{x1D706})$ and $\unicode[STIX]{x1D702}(\unicode[STIX]{x1D706})$ are the structure elastic modulus and structure viscosity. In the present work, the dependence of these material properties with respect to the structure parameter are given by

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}(\unicode[STIX]{x1D706})=\left[\left(\frac{\unicode[STIX]{x1D702}_{0}}{\unicode[STIX]{x1D702}_{\infty }}\right)^{\unicode[STIX]{x1D706}}-1\right]\unicode[STIX]{x1D702}_{\infty }, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle G(\unicode[STIX]{x1D706})=G_{0}\,f(\unicode[STIX]{x1D706})=G_{0}\exp \left[m\left(\frac{1}{\unicode[STIX]{x1D706}}-1\right)\right], & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D702}_{0}$ and $G_{0}$ are the values of viscosity and elastic modulus in the fully structured state, respectively, and $m$ is a positive constant.

2.3.2 Evolution of the structure parameter

For the evolution for the structure parameter, $\unicode[STIX]{x1D706}$ , an additional equation is generally formulated specifying the rate of formation and rate of destruction of micro-structure. The following equation is obtained from simple assumptions in these terms with the stress being identified as the destructing agent, i.e.

(2.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\unicode[STIX]{x1D706})=\frac{1}{t_{eq}}\left(1-\frac{\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}_{eq}(\unicode[STIX]{x1D70E})}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{eq}(\unicode[STIX]{x1D70E})$ is a function of the current stress intensity value. A fundamental concept of the determination of $\unicode[STIX]{x1D706}_{eq}(\unicode[STIX]{x1D70E})$ is the election of the flow curve as a set of steady-state points which is the locus of the attractor equilibrium points of the current non-equilibrium state. The thixotropic characteristic time, $t_{eq}$ , is the inverse of the rate of structure change in the state of maximum destruction. When $t_{eq}$ is high, thixotropic effects are more pronounced, i.e. a higher delay in structure change would occur as a consequence of a change in stress. Using the intensity of the deviatoric part of the stress tensor computed as

(2.6) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{dev}=\sqrt{{\textstyle \frac{1}{2}}\,\text{tr}\left[\left(\unicode[STIX]{x1D748}-{\textstyle \frac{1}{3}}\,\text{tr}(\unicode[STIX]{x1D748})\unicode[STIX]{x1D644}\right)^{2}\right]},\end{eqnarray}$$

and inputting this value in the flow curve, we can find an equilibrium deformation rate value, $\dot{\unicode[STIX]{x1D6FE}}_{eq}$ , corresponding to the current stress tensor. The equilibrium viscosity is computed by

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{eq}=\frac{\unicode[STIX]{x1D70E}_{dev}}{\dot{\unicode[STIX]{x1D6FE}}_{eq}}.\end{eqnarray}$$

Finally, we determine $\unicode[STIX]{x1D706}_{eq}(\unicode[STIX]{x1D70E}_{dev})$ corresponding to $\unicode[STIX]{x1D702}_{eq}$ . This structure parameter is obtained through the inverse of the viscosity function expressed by (2.3)

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{eq}=\frac{\ln \unicode[STIX]{x1D702}_{eq}-\ln \unicode[STIX]{x1D702}_{\infty }}{\ln \unicode[STIX]{x1D702}_{0}-\ln \unicode[STIX]{x1D702}_{\infty }}.\end{eqnarray}$$

The value of the current structure parameter that feeds (2.5) is computed by

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\frac{\ln [\unicode[STIX]{x1D702}(\unicode[STIX]{x1D706})+\unicode[STIX]{x1D702}_{\infty }]-\ln \unicode[STIX]{x1D702}_{\infty }}{\ln \unicode[STIX]{x1D702}_{0}-\ln \unicode[STIX]{x1D702}_{\infty }},\end{eqnarray}$$

where

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D702}(\unicode[STIX]{x1D706})+\unicode[STIX]{x1D702}_{\infty }=\frac{\unicode[STIX]{x1D70E}_{dev}}{\dot{\unicode[STIX]{x1D6FE}}}.\end{eqnarray}$$

2.3.3 Flow curve

The flow curve is generally obtained from a viscometric rheometry and constitutes a relation between the shear stress, $\unicode[STIX]{x1D70F}$ , and the deformation rate at steady-state conditions, i.e. equilibrium conditions. In the present work we will adopt the following equation that has shown to be useful in fitting rheometric data

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D701}(\dot{\unicode[STIX]{x1D6FE}})(\unicode[STIX]{x1D70F}_{y}+K\dot{\unicode[STIX]{x1D6FE}}^{n})+\unicode[STIX]{x1D702}_{\infty }\dot{\unicode[STIX]{x1D6FE}},\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{y}$ is the apparent yield stress of the material, $K$ and $n$ are the consistency and exponent of the power-law term in the flow curve, and $\unicode[STIX]{x1D702}_{\infty }$ is the viscosity at high shear rates. The function $\unicode[STIX]{x1D701}(\dot{\unicode[STIX]{x1D6FE}})$ is a regularization factor given by

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D701}(\dot{\unicode[STIX]{x1D6FE}})=\left\{\begin{array}{@{}ll@{}}1, & \text{if}~\unicode[STIX]{x1D70F}\geqslant \unicode[STIX]{x1D70F}_{y},\\ {\displaystyle \frac{(\unicode[STIX]{x1D702}_{0}-\unicode[STIX]{x1D702}_{\infty })\dot{\unicode[STIX]{x1D6FE}}}{\unicode[STIX]{x1D70F}_{y}+K\dot{\unicode[STIX]{x1D6FE}}^{n}}}, & \text{if}~\unicode[STIX]{x1D70F}<\unicode[STIX]{x1D70F}_{y},\end{array}\right.\end{eqnarray}$$

which is equivalent to considering a bi-viscosity model.

2.4 Governing equations and boundary conditions

For the isothermal motion of an incompressible material, the governing equations for mass and momentum are given by

(2.13) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left[\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{u})\right]=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}+\unicode[STIX]{x1D70C}\boldsymbol{g}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D748}$ is the total stress and $\boldsymbol{g}$ is the acceleration due to gravity.

The adopted boundary conditions for solving (2.13)–(2.14) are

(2.15) $$\begin{eqnarray}\boldsymbol{u}=\mathbf{0},\quad \text{on}~\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{d},\end{eqnarray}$$

and

(2.16) $$\begin{eqnarray}\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{t},\quad \text{on}~\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{t},\end{eqnarray}$$

where $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{d}$ is the boundary where the velocity $\boldsymbol{u}$ is prescribed and $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{t}$ is the domain boundary for imposing the traction condition. In (2.16), $\boldsymbol{n}$ denotes the outward unit normal vector to the boundary $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{t}$ while $\boldsymbol{t}$ is the traction vector. From (2.16) we can obtain the normal stress and tangential stress boundary conditions, given by

(2.17) $$\begin{eqnarray}\displaystyle & \boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}^{\text{T}}=\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{t}=F_{n},\quad \text{on}~\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{t}, & \displaystyle\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & \boldsymbol{m}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}^{\text{T}}=\boldsymbol{m}\boldsymbol{\cdot }\boldsymbol{t}=F_{m},\quad \text{on}~\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{t}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{m}$ is the unit tangent vector while $F_{n}$ and $F_{m}$ are known functions according to the physical problem which will be simulated.

For the structure parameter we solved (2.5) in $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{t}$ whilst homogeneous Neumann boundary conditions are employed at the solid wall.

2.5 Dimensionless formulation

Due to the complexity of studying the combination of the rheological properties of an EVPT material and the flow properties of the impacting drop problem, it is convenient to describe the governing equations in a non-dimensional form. To begin with, we define the following dimensionless variables

(2.19a-g ) $$\begin{eqnarray}t^{\ast }=t\frac{U_{c}}{L_{c}},\quad \unicode[STIX]{x1D735}^{\ast }=L_{c}\unicode[STIX]{x1D735},\quad \boldsymbol{u}^{\ast }=\frac{\boldsymbol{u}}{U_{c}},\quad {\unicode[STIX]{x1D749}^{M}}^{\ast }=\frac{\unicode[STIX]{x1D749}^{M}}{\unicode[STIX]{x1D70F}_{c}},\quad p^{\ast }=\frac{p}{\unicode[STIX]{x1D70C}U_{c}^{2}},\quad \boldsymbol{g}^{\ast }=\frac{\boldsymbol{g}}{g},\quad \unicode[STIX]{x1D63F}^{\ast }=\unicode[STIX]{x1D63F}\frac{L_{c}}{U_{c}},\end{eqnarray}$$

where $t$ is the time, $\unicode[STIX]{x1D735}$ is the gradient operator, $\boldsymbol{u}$ is the velocity vector field. The characteristic stress $\unicode[STIX]{x1D70F}_{c}$ is defined (see Thompson & Soares Reference Thompson and Soares2016) by evaluating the stress in the non-regularized version of the extended Hershel–Bulkley model of (2.11). Therefore,

(2.20) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{c}=\unicode[STIX]{x1D70F}_{y}+K\left(\frac{U_{c}}{L_{c}}\right)^{n}+\unicode[STIX]{x1D702}_{\infty }\frac{U_{c}}{L_{c}},\end{eqnarray}$$

where the length and velocity scales $L_{c}$ , respectively, $U_{c}$ , are chosen to be the normal impact velocity $U_{i}$ , respectively, the drop diameter  $D_{0}$ . Next, the dimensionless numbers related to the inertia, gravity, elasticity, plasticity and thixotropy of the fluid are defined by

(2.21) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle Re=\frac{\unicode[STIX]{x1D70C}U_{i}^{2}}{\unicode[STIX]{x1D70F}_{c}}~(\mathit{Reynold}\,\mathit{number}),\quad Fr=\frac{U_{i}}{\sqrt{gD_{0}}}~(\mathit{Froude}\,\mathit{number}),\\ \displaystyle \unicode[STIX]{x1D70F}_{y}^{\ast }=\frac{\unicode[STIX]{x1D70F}_{y}}{\unicode[STIX]{x1D70F}_{c}}~(\mathit{plastic}\,\mathit{number}),\\ \displaystyle Wi=\frac{\unicode[STIX]{x1D70F}_{c}}{G_{0}}~(\mathit{Weissenberg}\,\mathit{number}),\quad t_{eq}^{\ast }=t_{eq}U_{i}/D_{0}~(\mathit{Thixotropic}\,\mathit{number}),\end{array}\right\}\end{eqnarray}$$

Remark.

Note that (2.20) can be reformulated as

(2.22) $$\begin{eqnarray}1=\unicode[STIX]{x1D70F}_{y}^{\ast }+K^{\ast }+\unicode[STIX]{x1D702}_{\infty }^{\ast },\end{eqnarray}$$

where

(2.23a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{y}^{\ast }=\frac{\unicode[STIX]{x1D70F}_{y}}{\unicode[STIX]{x1D70F}_{c}},\quad K^{\ast }=\frac{K}{\unicode[STIX]{x1D70F}_{c}}\left(\frac{U_{c}}{L_{c}}\right)^{n},\quad \unicode[STIX]{x1D702}_{\infty }^{\ast }=\frac{\unicode[STIX]{x1D702}_{\infty }}{\unicode[STIX]{x1D702}_{c}},\end{eqnarray}$$

with the characteristic viscosity defined as $\unicode[STIX]{x1D702}_{c}\equiv \unicode[STIX]{x1D70F}_{c}L_{c}/U_{c}$ . Therefore, in the present model, $\unicode[STIX]{x1D702}_{\infty }^{\ast }$ denotes the unstructured dimensionless value of the viscosity, while the consistency index is denoted by  $K^{\ast }$ .

At this point, it is important to notice that the arguments for choosing $\unicode[STIX]{x1D70F}_{c}$ as a characteristic viscous stress to compose the Reynolds number can be used to support the choice adopted for the Weissenberg number. This happens because the elastic Mach number, $M$ , that in viscoelastic fluids can be obtained from $M=\sqrt{Re\,Wi}$ can be computed by the same expression in elasto-viscoplastic fluids. The elastic number $El=Wi/Re$ is also embedded with the same physical meaning as conceived, if we adopt as a relaxation time $\unicode[STIX]{x1D702}_{c}/G_{0}$ .

Finally, omitting the asterisk, the system of equations for the incompressible and isothermal free surface flow used in the present paper is

(2.24) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\qquad & \displaystyle\end{eqnarray}$$

(2.25) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{u})=-\unicode[STIX]{x1D735}p+\frac{1}{Re}(\unicode[STIX]{x1D702}_{\infty }\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{M})+\frac{1}{Fr^{2}}\boldsymbol{g},\qquad & \displaystyle\end{eqnarray}$$

(2.26) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D749}^{M}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\unicode[STIX]{x1D749}^{M})-[(\unicode[STIX]{x1D735}\boldsymbol{u})\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{M}+\unicode[STIX]{x1D749}^{M}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}]=\frac{f(\unicode[STIX]{x1D706})}{Wi}\left[-\frac{1}{\unicode[STIX]{x1D702}(\unicode[STIX]{x1D706})}\unicode[STIX]{x1D749}^{M}+2\unicode[STIX]{x1D63F}\right]\!.\qquad & \displaystyle\end{eqnarray}$$

As expected for a Weissenberg number, $Wi\rightarrow 0\Rightarrow \unicode[STIX]{x1D749}^{M}\rightarrow 2\unicode[STIX]{x1D702}(\unicode[STIX]{x1D706})\unicode[STIX]{x1D63F}$ , i.e. the fluid is purely viscous.

The form of the dimensionless equations associated with the evolution of the structure parameter, the mapping between viscosity and structure parameter (current and equilibrium ones), and the stress intensity, have exactly the same form as their dimensional counterparts.

The boundary conditions in dimensionless form are given by:

1. (i) No-slip boundary conditions

   (2.27a,b ) $$\begin{eqnarray}u=0\quad \text{and}\quad v=0.\end{eqnarray}$$

2. (ii) Free surface stress conditions

   (2.28) $$\begin{eqnarray}\displaystyle & & \displaystyle -p+\frac{2\unicode[STIX]{x1D702}_{\infty }}{Re}\left[\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}n_{x}^{2}+\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}y}n_{y}^{2}+\left(\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}x}\right)n_{x}n_{y}\right]\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{1}{Re}[(\unicode[STIX]{x1D70F}_{xx}^{M})n_{x}^{2}+2(\unicode[STIX]{x1D70F}_{xy}^{M})n_{x}n_{y}+(\unicode[STIX]{x1D70F}_{yy}^{M})n_{y}^{2}]=0;\end{eqnarray}$$
   and
   (2.29) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x1D702}_{\infty }}{Re}\left[2\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}n_{x}m_{x}+\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}y}n_{y}m_{y}+\left[\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}x}\right](n_{y}m_{x}+n_{x}m_{y})\right]\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{1}{Re}[((\unicode[STIX]{x1D70F}_{xx}^{M})-(\unicode[STIX]{x1D70F}_{yy}^{M}))n_{x}n_{y}+(\unicode[STIX]{x1D70F}_{xy}^{M})(n_{y}^{2}-n_{x}^{2})]=0,\end{eqnarray}$$
   where $\boldsymbol{n}=(n_{x},n_{y})$ denotes the outward unit normal vector to the boundary while $\boldsymbol{m}=(m_{x},m_{y})$ is the unit tangent vector.

Figure 3.

The elasto-viscoplastic thixotropic model of the present work is compared with the experimental and model results provided by Luu & Forterre (Reference Luu and Forterre2009) for the kaolin K51.

2.6 Overview of the computational technique

In this paper, we use a numerical scheme which combines the finite difference discretization and the marker-and-cell method implemented with in-house code. In summary, the two relevant stages of the numerical approach are as follows:

1. (i) Solving the governing equations. The conservation equations (2.24)–(2.25) are solved via a projection method which allows the solution in a segregated form. The constitutive equation (2.26) and the dimensionless counterpart of the evolution equation given by (2.5) are time discretized using the second-order Runge–Kutta scheme.

2. (ii) Moving free surface interface. In the context of the one-fluid model, the front-tracking method is employed here for dealing with the movement of the drop interface. The marker particles that describe the free surface are convected solving the equation

   (2.30) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{x}}{\text{d}t}=\boldsymbol{u}.\end{eqnarray}$$

More details about our code can be found in Oishi et al. (Reference Oishi, Thompson and Martins2016) while a mesh verification of the numerical scheme is presented in the appendix A.

4.1 Analysis of evolution of the height and diameter of the drop

Initially, we have concentrated our parametric studies on considering the influence of the elasticity, plasticity and inertia, fixing the thixotropic equilibrium time at $t_{eq}=7.5$ . The Weissenberg numbers have been selected to be $Wi=0.05$ , 0.5, 2, while the yield stress numbers have been selected to be $\unicode[STIX]{x1D70F}_{y}=0.1$ , 0.3, 0.5, 0.7. In the first round of simulations, we have considered lower normal impact velocities, which leads to the Reynolds numbers $Re=2.05$ , 3.36, 5.71. The temporal evolution of the minimum height and the spread factor are presented in figures 7 and 8, respectively. The results are displayed in a Cartesian manner with the columns corresponding to different Reynolds numbers, while the rows correspond to the different Weissenberg numbers. The different curves inside each figure of $H(t)/D_{0}$ and $D(t)/D_{0}$ correspond to changes in the plastic number. This manner of displaying the results enables the reader to go from one subfigure to another by different paths, providing insights into the combined behaviour. The analysis conducted below does not cover the universe of all possibilities.

Figure 7.

Temporal evolution of the normalized minimum bounce height $H(t)/D_{0}$ for $t_{eq}=7.5$ varying the following dimensionless parameters: $Wi=0.05$ , 0.5, 2 (rows), $Re=2.05$ , 3.36, 5.71 (columns) and $\unicode[STIX]{x1D70F}_{y}=0.1$ , 0.3, 0.5, 0.7 (inside of the figures).

Figure 8.

Temporal evolution of the normalized drop diameter $D(t)/D_{0}$ for $t_{eq}=7.5$ varying the following dimensionless parameters: $Wi=0.05$ , 0.5, 2 (rows), $Re=2.05$ , 3.36, 5.71 (columns) and $\unicode[STIX]{x1D70F}_{y}=0.1$ , 0.3, 0.5, 0.7 (inside of the figures).

From figure 7 one can quickly get some rough information as to whether the drop bounces or not after the first impact. In this connection, we can see that increasing the Weissenberg number favours the scenario where the drop attaches to the floor. In the lowest row, corresponding to the highest Weissenberg numbers ( $Wi=2$ ), there are no situations where the drop bounces for the present value of $t_{eq}$ irrespective of the Reynolds or plastic numbers. From the results in the first row, we can see that the first rebound can achieve a significant height, especially for the lowest plastic number, $\unicode[STIX]{x1D70F}_{y}=0.1$ (red line). High values of $Wi$ imply high relative flexibility, hence, high elastic deformations. It is not a coincidence that the lowest row of the dual figure, figure 8, subfigures 8(g), 8(h) and 8(i), have higher maximum spreading factors when compared to their counterparts with lower Weissenberg numbers. The Reynolds number has a positive effect on the maximum spreading, due to the initial kinetic energy of the system, but the effect of the plastic number is low. Another major tendency that is shown in figure 7 is that increasing the plastic number also contributes to the absence of bouncing. It is tempting to attribute this to an increase in dissipation alone, which could rise due to an increase of the relative yield stress. However, because the plastic number only changes the balance with respect to the weights of yield stress, the power-law term, and the constant viscous term, perhaps another effect contributes to produce this result. In fact, from (2.26), we can see that the effective Weissenberg number, $Wi_{eff}(\unicode[STIX]{x1D706})\equiv Wi/f(\unicode[STIX]{x1D706})$ , is local and takes into account the structural level of the fluid at each point. Since lower values of the plastic number represent situations in which it is easier to break down the fluid, this induces higher values of the elastic modulus, which in turn lowers the effective Weissenberg number. As discussed previously, this favours bouncing at this level of  $t_{eq}$ .

It is interesting to note that the influence of the Reynolds number is non-monotonic for $H(t)/D_{0}$ . Comparing subfigures 7(a), 7(b) and 7(c), we can see at least two effects that are worth mentioning. The first is that from $Re=2.05$ to 3.36, this low $Wi$ number produces opposite tendencies, depending on the plastic number. For example, this increase in $Re$ leads to stronger bouncing for low $\unicode[STIX]{x1D70F}_{y}$ (0.1 or 0.3), in the sense that the subsequent rebounds achieve a higher $H/D_{0}$ . On the other hand, for $\unicode[STIX]{x1D70F}_{y}=0.5$ , the same change in the Reynolds number leads to weaker bouncing. A second point worth noticing is the fact the top and the middle rows of figure 7 exhibit stronger bouncing at intermediate Reynolds numbers, $Re=3.36$ , when compared to $Re=2.05$ and 5.71 for the lower plastic number values $\unicode[STIX]{x1D70F}_{y}=0.1$ and 0.3. In other words, a non-monotonic effect is found for the increase of  $Re$ . The increase in the kinetic energy that is gained from $Re=2.05$ to 3.36 is stored more as an elastic energy for lower values of  $\unicode[STIX]{x1D70F}_{y}$ . However, a subsequent increase in the Reynolds number to $Re=5.71$ inhibits the bouncing of the fluids with different plastic number, so that the $\unicode[STIX]{x1D70F}_{y}=0.1$ case has a lower maximum height, while in the other cases there is no bouncing. These results can be explained by the following rationale. The energy per unit volume that is stored due to elastic effects is, roughly, given by $E_{el}=0.5\,G\unicode[STIX]{x1D6FE}_{e}^{2}$ , where $\unicode[STIX]{x1D6FE}_{e}$ is the elastic deformation. Lower values of the plastic number induces a higher elastic modulus and lower elastic deformation. With a small increase in the Reynolds number, the effect of increasing $G$ overcomes the decrease in the elastic deformation. With a higher increase in the Reynolds number, the breakdown process increases the dissipation and a smaller portion of the deformation is elastic.

From a general analysis of figure 8, we can observe that for most of the cases, the first spreading stage is not affected by the plastic number. This can be inferred by the spreading factor that, from the moment the drop reaches the floor until the moment where the maximum spreading factor is achieved, there is a coincidence among the different plastic number cases. The exception comes with higher values of the Reynolds and plastic numbers, where the drop spreads more intensively and the high value of the plastic number can prevent the receding stage (see the case where $\unicode[STIX]{x1D70F}_{y}=0.7$ in figure 8 c). In general, after achieving the maximum spreading factor, the receding stage starts with the different plastic number values still over the same master curve. Gradually, from the highest to the lowest plastic number value, the spreading factor curves abandon this master curve. This happens because higher values of the plastic number imply a higher resistance to receding, which inhibits the drop from returning to its original configuration. From figures 7(a) and 8(a) we can see that the different cases where bouncing occurs have similar contact times during the first contact period and that the spreading factors are also similar until the drop abandons the floor for the first time. Since the height achieved by the different fluids is not the same, this tendency does not happen for subsequent rebounds. In fact, we can see a decrease in the number of rebounds in this case of low $Re$ and low  $Wi$ : three for $\unicode[STIX]{x1D70F}_{y}=0.1$ , two for $\unicode[STIX]{x1D70F}_{y}=0.3$ , one for $\unicode[STIX]{x1D70F}_{y}=0.5$ and zero for $\unicode[STIX]{x1D70F}_{y}=0.7$ . The effect of the plastic number on the spreading factor is monotonic and well defined in the $Wi=0.5$ – $Re=5.71$ case, where a clear separation between the curves occurs after the maximum spreading. The higher resistance of fluids with higher plastic numbers leads to oscillations at the floor around larger diameters. When the Weissenberg number is increased, the plastic number has a weak effect, showing that elasticity damps the effect of the yield stress. For the lowest $Wi$ , keeping $Re=5.71$ , we can see a decrease in the tendency of an oscillatory motion. The non-monotonic result of the case of the lowest plastic number can be attributed to the bouncing of this drop during the captured period. These general results are in accordance with some experimental data provided by Luu & Forterre (Reference Luu and Forterre2009), where the elastic character of the yield stress fluid mimics some of the results obtained by inelastic fluids where interfacial effects have an elastic nature.

The effects of thixotropy, combined with elasticity, and inertia are analysed considering a high value of the plastic number, $\unicode[STIX]{x1D70F}_{y}=0.7$ , shown in figures 9 and 10 for minimum height and spreading factor, respectively. For these investigations, we have taken various Weissenberg numbers, $Wi=0.05$ , 0.5, 2, various Reynolds numbers, $Re=2.05$ , 5.71, 13.49, and various thixotropic equilibrium times, $t_{eq}=5$ , 7.5, 10, 20.

The first thing that attracts attention when one examines figure 9 is that there are few cases with rebounds. However, this high a value of the plastic number inhibits rejuvenation, which in turn keeps the level of elasticity high. In the material analysed in the present paper, this implies higher elastic deformations, which can be verified from the levels of the spreading factor depicted in figure 10, which exceed significantly those obtained in when $\unicode[STIX]{x1D70F}_{y}=0.1$ (not shown here). The few bouncing results are associated to high values of  $t_{eq}$ , as discussed previously. When there is not enough time for the structure to change, the drop kicks just after it impacts the solid surface, as an elastic ball.

The lower $Wi$ row shows almost non-oscillating results. As the Reynolds number is increased, gradually the final configuration is changed from the maximum to the minimum spreading factor, since lower thixotropic equilibrium times induce the drop to remain at the maximum relative diameter, while higher values of $t_{eq}$ enable the fluid to retract; but once the minimum spreading factor is achieved, there is not enough energy to overcome the yield threshold. The intermediate Weissenberg number row is just a wavy version of the top row, where now elastic effects are able to induce oscillations in the drop. The case of high $Wi$ , high Reynolds number, and high $t_{eq}$ shows the beginning of nonlinear combinations of elastic, inertial and thixotropic effects.

Figure 9.

Temporal evolution of the normalized minimum bounce height $H(t)/D_{0}$ for $\unicode[STIX]{x1D70F}_{y}=0.7$ , varying the following dimensionless parameters: $Wi=0.05$ , 0.5, 2 (rows), $Re=2.05$ , 5.71, 13.49 (columns) and $t_{eq}=5$ , 7.5, 10, 20 (inside of the figures).

Figure 10.

Temporal evolution of the normalized drop diameter $D(t)/D_{0}$ for $\unicode[STIX]{x1D70F}_{y}=0.7$ , varying the following dimensionless parameters: $Wi=0.05$ , 0.5, 2 (rows), $Re=2.05$ , 5.71, 13.49 (columns) and $t_{eq}=5$ , 7.5, 10, 20 (inside of the figures).

4.2 Bouncing–non-bouncing outputs of the normal drop impact

A qualitative aspect of the problem that is of particular practical importance in the drop impact problem is if the drop will bounce or stick in the substrate. In order to reduce the four-dimensional parameter space of the dimensionless numbers analysed in this work, we followed the rationale described next.

Considering that the elasticity of the material is responsible for effects that in the Newtonian case are attributed to elastic interfacial causes, we can apply what was learned from the literature and elect a counterpart of the Ohnesorge number as a parameter to be analysed. From its original definition an elastic Ohnesorge number can be defined as $Oh_{e}\equiv \sqrt{Wi/Re}$ , which preserves a lack of dependence with the characteristic velocity. A second combination of dimensionless numbers is obtained from the observation of the results described in the previous subsection. For the range of the tested parameters, we could identify a non-bouncing tendency when $Re$ , $Wi$ and $\unicode[STIX]{x1D70F}_{y}$ are increased. The opposite trend would happen by increasing  $t_{eq}$ . The importance of the elastic Mach number, $M\equiv U/\sqrt{G/\unicode[STIX]{x1D70C}}=\sqrt{Wi\,Re}$ , highlighted in Luu & Forterre (Reference Luu and Forterre2009), lead us to construct the following dimensionless number

(4.1) $$\begin{eqnarray}Y=\frac{M\unicode[STIX]{x1D70F}_{y}}{t_{eq}}.\end{eqnarray}$$

The ratio $\unicode[STIX]{x1D70F}_{y}/t_{eq}$ can be interpreted as a dissipation per unit volume due to the plastic nature of the material during a time scale equivalent to the thixotropic equilibrium time of the material. Therefore, it is reasonable to expect that increasing this ratio, the material has a tendency to avoid bouncing. In figure 11, we can see a plot of $Y$ versus $Oh_{e}$ for the normal impact case. Rebound and non-rebound outputs were separated, although not perfectly. A critical value for the elastic Ohnesorge number was identified, i.e. for $Oh_{ec}>0.5$ no bouncing occurred (results for higher $Oh_{e}$ values not shown in this figure corroborate this finding). For $Y>\approx 0.08$ no bouncing results are verified, but we can identify some non-bouncing results below this value. In the appendix B, we provide maps of the bouncing–non-bouncing outputs for the four dimensionless numbers analysed in the present work.

Figure 11.

Normal impact.

4.3 Analysis of the evolution of the shape and structure of the drop

The last set of results of the normal drop impact subsection are some selected cases that show the temporal evolution of the position and shape of the drops, filled with the structure parameter field, $\unicode[STIX]{x1D706}$ . The $\unicode[STIX]{x1D706}$ domain, $\unicode[STIX]{x1D706}\in [0,1]$ , was divided equally into seven parts and a different colour was designated for each part. The colours range from dark blue (corresponding to the lowest values of  $\unicode[STIX]{x1D706}$ ) to a yellow mustard colour (corresponding to the highest values of  $\unicode[STIX]{x1D706}$ ). The colour code is illustrated at the top of each of figures 12–14. Before impact, all drops present a uniform distribution of $\unicode[STIX]{x1D706}=1$ (fully structured state). The main purpose of showing this last set of results is to highlight the importance of thixotropy in the dynamics of the drop impact. It is interesting to note that, by keeping the dimensionless parameters $Re$ , $Wi$ , $\unicode[STIX]{x1D70F}_{y}$ fixed and changing the thixotropic time scale of the material, very different evolutions of the shape can happen as a consequence of the modified dynamics of the structure. In common, these figures contrast a high and a low value for the thixotropic equilibrium time, $t_{eq}=5$ and 20.

Figure 12 also shows an elastic effect by varying the Weissenberg number. The Reynolds number is low and $\unicode[STIX]{x1D70F}_{y}=0.7$ . At $t=3.2$ , just after impact, we can see that for $Wi=2$ (k,l), the material deforms to a higher extent than for drops with $Wi=0.05$ (a,b). The structural parameter of the cases in the left column, where $t_{eq}=5$ , are clearly lower than the cases where the thixotropic equilibrium time is higher, in the right column. This result is particularly significant in the case with low elasticity ( $Wi=0.05-t_{eq}=5$ ), since less energy is stored (compared to the case with $Wi=2-t_{eq}=5$ ) and therefore the breakdown process is more intense. When $t=4.6$ , the high- $t_{eq}$ -low- $Wi$ case (b) shows a rebound. Although the Weissenberg number was low, which induced a low degree of deformation, rejuvenation was retarded by the higher value of $t_{eq}=20$ . Therefore, the elastic effects were still present, which enabled the drop to bounce and continue breaking down without touching the solid. Due to a very low $Re$ , the difference in shape that happens as a consequence of a change in $t_{eq}$ is not significant.

Figure 12.

Numerical prediction of the structure parameter $\unicode[STIX]{x1D706}$ at different dimensionless times using $Re=2.05$ , $\unicode[STIX]{x1D70F}_{y}=0.7$ . From (column) left to right: $t_{eq}=5$ , 20 and from (row) top to bottom: $Wi=0.05$ , 2. The reader will find, in the supplementary material, available at https://doi.org/10.1017/jfm.2019.541, the file Movie 1 with the video of the evolution displayed in this figure.

Figure 13.

Numerical prediction of the structure parameter $\unicode[STIX]{x1D706}$ at different dimensionless times using $Re=13.49$ , $\unicode[STIX]{x1D70F}_{y}=0.7$ and $Wi=0.5$ . From (column) left to right: $t_{eq}=5,20$ . The reader will find in the supplementary material the video with the complete evolution of the motions shown in this figure.

A different scenario is presented in figure 13, where the evolution of the shape is dramatically dependent on the value of $t_{eq}$ . In this case, a high Reynolds number was imposed. This high kinetic energy is able to produce a stress level after impact that contributes to the breakdown process. At $t=2.0$ , when the drops are still in the spreading stage, we can observe a reasonable difference in the distribution of the structural parameter inside the left ( $t_{eq}=5$ ) and right ( $t_{eq}=20$ ) drops. As expected, the regions of the drop which are closer to the solid surface are less structured at this instant, since they were more exposed to higher stresses. When $t=2.5$ , the two drops exhibit similar spreading factors, very near the maximum ones. However, we can see a small difference in shape and significant discrepancy in structure distribution, favouring more structured states in the high $t_{eq}$ case. The thixotropic effect on the drop configuration is particularly important in the receding stage, i.e. after the drop has reached its maximum spreading factor. Examining $t=4.8$ in figure 13, we can see that before this moment, the left drop has consumed the kinetic energy by breaking the material and there is not enough stored energy to provoke retraction. In this fixed configuration, the thixotropic nature of the material induces an uniform build-up process in the left drop (which continues until $t=10$ ), since the material is exposed to an almost free stress condition. On the other hand, the right drop at $t=4.8$ has retracted due to elastic energy that was stored and that has been transformed into kinetic energy again. The edges of the drop initiate a motion away from the solid wall, which gradually weakens as time evolves (see $t=10$ ).

Figure 14.

Numerical predictions of the structural parameter $\unicode[STIX]{x1D706}$ at different dimensionless times using $Re=13.49$ , $\unicode[STIX]{x1D70F}_{y}=0.7$ and $Wi=2.0$ . Columns from left to right: $t_{eq}=5$ , 20. The reader will find in the supplementary material the video with the complete evolution of the motions shown in this figure.

When the Weissenberg number is increased from $Wi=0.5$ to 2.0, the evolutions of $t_{eq}$ take the form illustrated in figure 14. It is interesting to note that even when the speed of the rejuvenation process is high, such as when $t_{eq}=5$ (left drop), a high $Wi$ induces retraction after a maximum spreading obtained at $t\approx 6.3$ . We can see that the material remains relatively structured near the interface for this lower $t_{eq}$ value. Needless to say, this retraction is less intense than that for the $t_{eq}=20$ counterpart. Due to a rapid structure response, the analysis of $t=6.3$ reveals that during retraction, the left drop (which was experiencing low stress levels) initiates a build-up process. On the other hand, the higher $t_{eq}$ case at $t=6.3$ shows that the material continues to break down during retraction due to the stress inputs that were imposed before. The area of the drop that is in contact with the surface of the solid continues to decrease until $t=8.4$ , when the left and right drops exhibit an outward motion detaching themselves from the floor, upward in the $t_{eq}=5$ case and downward in the $t_{eq}=20$ case. From $t=8.4$ to 10, we can see that the high level of structure of the left drop has conferred a solid-like state that ‘freezes’ the material, while the right drop still has a liquid-like behaviour leading to a natural falling motion of the edges of the drop which were not in contact with the solid surface.

The reader will find, in the supplementary material, the file Movie 2 with the videos of the set of evolutions displayed in figures 13, 14.

5.1 Evolution of the shape and position

Figures 15 and 16 show two different dimensionless times, $t=2.5$ and 5.0, respectively, that correspond to the evolution of a number of cases that will be described next. In each subfigure, from (a) to (f), there are four different drops, each one corresponding to a different type of fluid, one Newtonian and three EVPT ones. The fluids are identified by the colour of their interfaces and common features of each subfigure are described in the captions of each figure. The colour code is black for the Newtonian drop, red for the drop with $t_{eq}=1$ , green for the drop with $t_{eq}=5$ and blue for the drop with $t_{eq}=20$ . In common, all subfigures are associated with cases of fixed Weissenberg number, $Wi=0.05$ . In addition to the effect of the thixotropic time, which can be analysed inside each subfigure, the comparison of the subfigures explores changes in $Re$ and plastic number.

Figure 15.

Drop interface profiles for a surface at an inclination of $\unicode[STIX]{x1D703}=30^{\circ }$ , at dimensionless time $t=2.5$ . The fixed dimensionless parameter for the EVPT material is $Wi=0.05$ ; The thixotropic times shown are $t_{eq}=1$ , 5, 20. The drop interface profiles considering a Newtonian fluid for different $Re$ are also included.

Figure 16.

Drop interface profiles for a surface at an inclination of $\unicode[STIX]{x1D703}=30^{\circ }$ , at dimensionless time $t=5$ . The fixed dimensionless parameter for the EVPT material is $Wi=0.05$ ; The times shown are $t_{eq}=1$ , 5, 20. The drop interface profiles considering a Newtonian fluid for different $Re$ are also included.

Figure 17.

Drop interface profiles for a surface at an inclination of $\unicode[STIX]{x1D703}=30^{\circ }$ at dimensionless time $t=2.5$ . The fixed dimensionless parameter for the EVPT material is $Wi=2$ ; The times shown are $t_{eq}=1$ , 5, 20. The drop interface profiles considering an Oldroyd-B fluid with $Wi=2$ for different $Re$ are also included.

Figure 18.

Drop interface profiles for a surface at an inclination of $\unicode[STIX]{x1D703}=30^{\circ }$ at dimensionless time $t=18$ . The fixed dimensionless parameter for the EVPT material is $Wi=2$ ; The times shown are $t_{eq}=1$ , 5, 20. The drop interface profiles considering an Oldroyd-B fluid with $Wi=2$ for different $Re$ are also included.

Taking as reference the instant $t=2.5$ , the left column, corresponding to subfigures 15(a), 15(c) and 15(e), have the same Reynolds number, $Re=2.05$ , while in the right column, corresponding to subfigures 15(b), 15(d) and 15(f), $Re=5.71$ . Changes in the plastic number are observed when one compares different rows. The top one, corresponding to subfigures 15(a) and 15(b), have $\unicode[STIX]{x1D70F}_{y}=0.1$ . In the middle row, figures 15(c) and 15(d), $\unicode[STIX]{x1D70F}_{y}=0.3$ ; while figures 15(e) and 15(f), have their plastic numbers fixed at $\unicode[STIX]{x1D70F}_{y}=0.7$ . This same structure is maintained in figure 16.

At $t=2.5$ , just after the impact, figure 15 shows that among the cases of low Reynolds numbers, corresponding to the left column, there is no big difference between the different drops. However, for $Re=5.71$ , we can see a more substantial difference imposed by the particular values of  $t_{eq}$ . While the lower value of  $t_{eq}$ , $t_{eq}=1$ , exhibits a similar behaviour to the Newtonian case, higher values of $t_{eq}$ are associated with less deformation and, in the case of $t_{eq}=20$ and lower values of  $\unicode[STIX]{x1D70F}_{y}$ , we can observe that a rebound occurs. It is important to note that the $t_{eq}=1$ case can be seen as a perturbation of the Newtonian case towards the EVPT material, since the dimensionless numbers that are responsible for elastic, plastic and thixotropic effects are low ( $Wi=0.05$ , $\unicode[STIX]{x1D70F}_{y}=0.1$ , $t_{eq}=1$ ). When the time is advanced to $t=5.0$ , we can see from the analysis of figure 16 that the blue drop corresponding to $t_{eq}=20$ bounces even when $Re=2.05$ , and the position and shape of the red drop ( $t_{eq}=1$ ) still gravitate around those of the Newtonian drop.

Another set of evolutions is depicted in figures 17 and 18, corresponding to $t=2.5$ and 18, respectively. Again, the left and right columns are identified with $Re=2.05$ and 5.71, respectively, and the top, middle and bottom rows are identified with $\unicode[STIX]{x1D70F}_{y}=0.1$ , $\unicode[STIX]{x1D70F}_{y}=0.3$ and $\unicode[STIX]{x1D70F}_{y}=0.7$ , respectively. The difference between this and the previous set is the Weissenberg number, which in the present set is $Wi=2$ . Another important change is that for comparison with the EVPT drops we chose as reference an Oldroyd-B viscoelastic fluid (with the same $Wi=2$ ). In figure 17, which corresponds to $t=2.5$ , we can see immediately that after this oblique impact, these drops deform more than the ones with low $Wi$ shown in figure 15. Again, we can see that in the right column the drops exhibit more spreading than in the left column, where the kinetic energy transferred to the drop is lower.

Figure 19.

Inclination angle of $15^{\circ }$ .

During the evolution until $t=18$ , rolling (with its variant, sliding) is the main result of this set of numerical experiments. The plastic and purely viscous deformations in the retraction process have consumed the elastic energy and no bouncing is to be seen. After retraction, some shapes achieve a stable configuration in that the stresses caused by gravity and the contact with the solid, together with the dynamics of aging and rejuvenation, find an equilibrium. The black drop (Oldroyd-B fluid) is not influenced by the plastic number, and therefore its shape and position do not change from one row to another. In this connection, we can see that the viscoplastic-thixotropic effects that change the viscous and elastic properties of the EVPT drops also change the dynamics of their evolution in a complex manner that makes their behaviour deviate from a purely elastic behaviour.

The reader will find, in the supplementary material, the files Movie 3 and Movie 4 with the videos of the set of evolutions displayed in this subsection.

5.2 Main outputs of the oblique drop impact

The complexity introduced by the oblique impact was not able to be captured by the quantities, $Oh_{e}$ and  $Y$ , that were used in the normal drop impact. An additional qualitative output was included: rolling. We found non-monotonic trends associated with the original dimensionless numbers, as shown in figures 19 and 20 for inclination angles of $15^{\circ }$ and $30^{\circ }$ , respectively. However, the use of $Y$ and $Oh_{e}$ could separate tendencies. We can say that, for low $Y$ values and low $Oh_{e}$ there is tendency to have a bouncing output. The rolling output is achieved for low $Y$ and high $Oh_{e}$ . The same critical elastic Ohnesorge number, $Oh_{ec}=0.5$ , confines bouncing results for both inclination angles. High $Y$ values lead to the drop to stick. For the highest thixotropic equilibrium time tested, $t_{eq}=20$ , we had bouncing and rolling results only. In addition, we found no dependence on the inclination angle for the outputs in this high $t_{eq}$ regime. Using $M/\unicode[STIX]{x1D6FE}_{y}$ , as in Luu & Forterre (Reference Luu and Forterre2009), to compose with $Oh_{e}$ the two-dimensional Cartesian space we were able to separate bouncing from rolling, as figure 21 reveals. It is not a coincidence that the same critical value for $Oh_{e}$ is found. Among the chosen $t_{eq}$ values, only in the inclined cases with the highest thixotropic equilibrium time, the results were clearly separated into two regions and a critical $Oh_{e}$ could be found; for $Oh_{ec}>0.5$ no bouncing occurred. The graphs corresponding to the inclination angles of $15^{\circ }$ and $30^{\circ }$ are shown in figure 21. In fact, for $t_{eq}=20$ there was no difference between the two tested inclination values.

Figure 20.

Inclination angle of $30^{\circ }$ .

Figure 21.

Inclination angle of $15^{\circ }$ or $30^{\circ }$ for $t_{eq}=20$ .

Figure 22.

Numerical prediction of the structural parameter $\unicode[STIX]{x1D706}$ at different dimensionless times using $Re=13.49$ , $\unicode[STIX]{x1D70F}_{y}=0.7$ , $t_{eq}=5$ and $Wi=0.5$ for a surface at an inclination of $\unicode[STIX]{x1D703}=30^{\circ }$ . The reader will find, in the supplementary material, the file Movie 5 with the video of the complete evolution of the case shown in this figure.

In appendix B we show the general outputs (stick, rolling and bouncing) according to the combinations of the values of the four dimensionless numbers explored in the present paper.

5.3 Evolution of the shape, position and structure for an oblique drop impact

Figures 22, and 23 show successive changes in time of the structural parameter field, together with the shape and position as the drop evolves after an oblique collision, with a high impact velocity ( $Re=13.49$ ). In this set of figures, the plane is exhibited with an inclination of $\unicode[STIX]{x1D703}=30^{\circ }$ and gravity acts in the vertical direction.

In figure 22, the drop exhibits a kind of ‘mushroom head’ shape just after impact, figure 22(a), where we can clearly see a much more intense rejuvenation in the regions near the wall because of the low value of  $t_{eq}$ , which induces rapid destruction. Since $Re$ is high, the drop evolves in the subsequent frame, figure 22(b), to a very stretched configuration with a very low degree of structure, quite uniformly distributed. In this configuration, the gravitational forces are not strong enough to induce more breaking and to overcome the overall viscous resistance: so, the drop almost does not change position, presenting slight changes in shape. At this stage, the build-up process begins and structuring takes place in the evolving drops.

Figure 23 shows a case where rebound occurs in the oblique impact problem. In this case, a reasonable deformation occurs during the spreading stage after impact, and the retraction process takes place with a high level of stored energy, enabling the drop to be expelled from the surface of the solid. In other words, this rebound is not a direct effect of the deformation associated with the direction normal to the inclined plane, but with a deformation, extension and retraction in the direction of the inclined plane. The instant after impact, shown in subfigure figure 23(a), has the same ‘mushroom-head’ shape, where the bottom region has less structure due to the breakdown process that occurs as a result of the impact. The next frame, figure 23(b), shows a ‘black bean’ shape, where a more severe breakdown has happened and the drop is still retracting. This retraction continues from both sides and induces an out-of-the-plane motion and in figure 23(c) the drop has just left the solid, leaving a nub near the point of its last contact. After leaving the plane, the drop does not change its shape significantly, developing a rotating trajectory and an aging process.

Figure 23.

Numerical prediction of the structural parameter $\unicode[STIX]{x1D706}$ at different dimensionless times using $Re=13.49$ , $\unicode[STIX]{x1D70F}_{y}=0.3$ , $t_{eq}=20$ and $Wi=0.5$ for a surface at an inclination of $\unicode[STIX]{x1D703}=30^{\circ }$ . The reader will find, in the supplementary material, the file Movie 6 with the video of the complete evolution of the case shown in this figure.