Numerical Model for the Cleaning of a Film-like Soil by Viscous Shifting (2025)

Model overview

The model is developed as a BCCM, which implies the standard assumptions for all BCCMs [Citation7,Citation11].

1. The height of the soil is negligible compared to the dimensions of the flow. Hence, the cleaning fluid’s flow is computed without considering the soil.

2. The cleaning progress does not influence the fluid forces acting on the soil

3. The height of the soil is negligible compared to the soil length ($h_{\mathrm{s}}\ll L_{\mathrm{s}})\text{.}$ With that, transport processes into the soil, e.g., swelling and heating, can be described using one-dimensional transport equations using the wall-normal direction.

Further assumptions will be made during the model description. For modeling, the process of viscous shifting is divided into the four subprocesses shown in . In the first step, the loads acting on the soil are calculated using a comparative stress ${\tau }_{\mathit{\text{hyd}}}\text{.}$ In the second step, the thermal behavior of the soil is modeled. This is an important aspect of viscous shifting, as temperature has a strong influence on the viscosity of most soils. In the third step, the soil rheology needs to be evaluated, which in general has the form ${\tau }_s=f({\dot{\gamma }}_s{,\vartheta }_s)\text{.}$ Here, ${\tau }_s$ is the shear stress present in the soil, ${\dot{\gamma }}_s$ the shear rate of the soil and ${\vartheta }_s$ is the soil temperature. In the final step, the computed quantities are combined to determine the movement of the soil. In the following sections, the modeling of the soil movement is discussed in detail, before modeling the thermal behavior.

Load calculation

The cleaning model requires the mechanical load acting on the soil as an input. The stress exerted by the flow is obtained from the CFD and is then imposed in the simulation of the soil. This amounts to one-way coupling between fluid and soil and is visualized in . The single-phase flow calculation uses the global coordinates $x$ and $y$ and provides the stress on the soil, without resolving its geometry. The evolution of the soil is simulated separately, accounting for the actual dimensions of the soil. It uses its own coordinates, $x_s$-$y_s\text{,}$ with $y_s$ resolving the thickness of the soil.

A single-phase flow solution, in particular the shear stress, can be obtained in different ways, by experimental measurements, by CFD simulations, or analytical solutions, if available. In a general flow situation, the load on the soil is constituted by a combination of pressure and shear force. In the present approach, this is limited to shear forces since the height of the soil is assumed to be small, so the side areas where the pressure force could attack are very small and, hence, negligible compared to the shear forces acting on the top surface of the soil. It is understood that scenarios exist in which pressure forces become relevant even in the case of film-like soils. Examples involve pressure oscillation or cavitation used for cleaning in some processes, but these situations are not considered here. In the present work a scalar reference value ${\tau }_{\mathit{\text{hyd}}}$ characterizing the magnitude of the hydrodynamic load is employed. It is determined by integrating the fluid forces along the interface between the fluid and the soil, termed $A_{\mathit{\text{int}}}$ here, to give(1) $${\tau }_{\mathrm{\text{hyd}}}=\frac{1}{A_{\mathrm{\text{int}}}}{\int }_{A_{\mathrm{\text{int}}}}^{\mathrm{}}\Vert \mathit{\tau }\cdot \mathit{n}\Vert \mathrm{d}A$$(1)

All quantities in EquationEquation (1)(1) $${\tau }_{\mathrm{\text{hyd}}}=\frac{1}{A_{\mathrm{\text{int}}}}{\int }_{A_{\mathrm{\text{int}}}}^{\mathrm{}}\Vert \mathit{\tau }\cdot \mathit{n}\Vert \mathrm{d}A$$(1) are depicted in . The magnitude is used so that shear stresses in opposing directions do not compensate each other during the integration across the interfacial area. In cases where pressure forces are relevant, these must be included in the calculation of the reference load.

Soil rheology

For the present model, it is assumed that the soil rheology can be stated in the form(2) $${\tau }_{\mathrm{s}}=f({\dot{\gamma }}_{\mathrm{s}},{\vartheta }_{\mathrm{s}})$$(2)

Here, ${\tau }_s$ is the shear stress present in the soil, ${\dot{\gamma }}_s$ the shear rate of the soil and ${\vartheta }_s$ the soil temperature.

Soil movement

To model the soil movement, the soil is discretized in flow direction into $N_{\mathit{\text{seg}}}$ segments. The discretization is depicted in . Only the flow component in the main flow direction is considered. The soil layer is assumed to be very thin. Therefore, the velocity profile in each segment $i$ is approximated with a linear function reading(3) $$u_{\mathrm{s},i}\left(y_{\mathrm{s}}\right)={\dot{\gamma }}_{\mathrm{s},i}{y}_{\mathrm{s}}$$(3) where ${\dot{\gamma }}_{s,i}$ is the shear rate in segment $i$ and $y_s$ the wall-normal coordinate. This simplification is not valid if the soil layer height is too high or if effects like wall slip become dominant. The shear stresses are assumed to be in balance at the interface of the soil layer and cleaning fluid. Assuming a soil rheology of the form stated in EquationEquation (2)(2) $${\tau }_{\mathrm{s}}=f({\dot{\gamma }}_{\mathrm{s}},{\vartheta }_{\mathrm{s}})$$(2) , the shear rate is determined by inserting the comparative stress in the rheology law (${\tau }_s={\tau }_{\mathit{\text{hyd}}}$) and solving for ${\dot{\gamma }}_s\text{.}$ Computing the soil temperature ${\vartheta }_s$ is discussed in the following section. Since the velocity profile $u_{s,i}(y_s)$ is known, it can be integrated over the cross-sectional area $A_{\mathit{\text{out}},i}$ to obtain the outgoing mass flow of a segment(4) $${\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i}={\rho }_{\mathrm{s}}{\int }_{A_{\mathrm{\text{out}},i}}^{}{u}_{\mathrm{s},i}\mathrm{d}A$$(4)

Herein, ${\rho }_s$ is the soil density, and $A_{\mathit{\text{out}},i}$ is the surface area of the outlet of cell $i\text{.}$ The incoming mass flow in a segment is assumed to be the mass flow leaving the segment directly upstream, i.e.,(5) $${\dot{m}}_{\mathrm{s},\mathrm{\text{in}},i}=\{\begin{array}{l}0i=1,\\ {\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i-1}\mathit{\text{else}}\end{array}$$(5)

The change of mass in a segment is expressed by(6) $$\frac{\mathrm{d}{m}_{\mathrm{s},i}}{\mathrm{d}t}={\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i}-{\dot{m}}_{\mathrm{s},\mathrm{\text{in}},i}$$(6)

EquationEquation (5)(5) $${\dot{m}}_{\mathrm{s},\mathrm{\text{in}},i}=\{\begin{array}{l}0i=1,\\ {\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i-1}\mathit{\text{else}}\end{array}$$(5) is discretized using the Euler forward method, yielding(7) $$m_{\mathrm{s},i}^{(n+1)}=m_{\mathrm{s},i}^{\left(n\right)}-\mathrm{\Delta }t\left({\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i}^{\left(n\right)}-{\dot{m}}_{\mathrm{s},\mathrm{\text{in}},i}^{\left(n\right)}\right)$$(7) where $\mathrm{\Delta }t$ is the time step size and the superscript $n$ refers to the approximation of a quantity at a time $t^{(n)}=n\mathrm{\Delta }t.$ This discretization imposes a stability criterion. It must be ensured that the mass within a segment is always larger than the mass being removed by EquationEquation (7)(7) $$m_{\mathrm{s},i}^{(n+1)}=m_{\mathrm{s},i}^{\left(n\right)}-\mathrm{\Delta }t\left({\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i}^{\left(n\right)}-{\dot{m}}_{\mathrm{s},\mathrm{\text{in}},i}^{\left(n\right)}\right)$$(7) . This causes the segments to never be empty, i.e., clean. One way to overcome this issue is to define a critical soil height where a segment is considered clean. The numerical value of this threshold needs to be defined problem-specific. Another workaround would be the combination with the model for adhesive detachment [Citation7,Citation12], which would define a physical criterion for the detachment of the soil from the substrate. A combination of both models is presented in [Citation42].

The way the mass of the soil segment $m_{\mathrm{s},i}$ is related to the segment height $h_{\mathrm{s},i}\text{,}$ and the outlet area $A_{\mathit{\text{out}},i}$ depends on the investigated geometry. Throughout the derivation, a unique flow direction was assumed. This is not met in all practical applications. However, the assumption does not constitute a major restriction since a more general formulation involving the flow direction can be derived using the concept of fluxes known from finite volume methods.

Thermal modeling

This section aims to derive equations that can be used to compute the evolution of the soil temperature over time. The thermal model is designed for application to the flushing processes of chocolate, where insulated pipes are used. Hence, the walls are assumed to be adiabatic. The temperature in each soil segment is assumed to be uniform. The initial phase with a non-uniform temperature profile is neglected in the present framework. A more general approach accounting for non-uniform temperature distributions and non-adiabatic walls is provided in [Citation42].

The enthalpy $H_{s,i}$ stored in each segment can be written as(8) $$H_{\mathrm{s},i}=m_{\mathrm{s},i}{c}_{p,\mathrm{s}}{\vartheta }_{\mathrm{s},i}$$(8) where $c_{p,s}$ is the heat capacity of the soil and a reference state of ${\vartheta }_{\mathit{\text{ref}}}=0°C$ is assumed. When the soil is moving, convective enthalpy transport occurs, which can be computed analogously to EquationEquation (7)(7) $$m_{\mathrm{s},i}^{(n+1)}=m_{\mathrm{s},i}^{\left(n\right)}-\mathrm{\Delta }t\left({\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i}^{\left(n\right)}-{\dot{m}}_{\mathrm{s},\mathrm{\text{in}},i}^{\left(n\right)}\right)$$(7) as(9) $${\tilde{H}}_{\mathrm{s},i}^{(n+1)}=H_{\mathrm{s},i}^n-\mathrm{\Delta }t\left({\dot{H}}_{\mathrm{s},\mathrm{\text{out}},i}^{\left(n\right)}-{\dot{H}}_{\mathrm{s},\mathrm{\text{in}},i}^{\left(n\right)}\right)$$(9) where ${\tilde{H}}_{s,i}^{n+1}$ is an intermediate value for the thermal energy since it does not account for heat convection from the main flow, yet. The enthalpy fluxes in EquationEquation (9)(9) $${\tilde{H}}_{\mathrm{s},i}^{(n+1)}=H_{\mathrm{s},i}^n-\mathrm{\Delta }t\left({\dot{H}}_{\mathrm{s},\mathrm{\text{out}},i}^{\left(n\right)}-{\dot{H}}_{\mathrm{s},\mathrm{\text{in}},i}^{\left(n\right)}\right)$$(9) are computed as(10) $${\dot{H}}_{\mathrm{s},\mathrm{\text{out}},i}={\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i}{c}_{p,\mathrm{s}}{\vartheta }_{\mathrm{s},i}$$(10) and(11) $${\dot{H}}_{\mathrm{s},\mathrm{\text{in}},i}=\{\begin{array}{l}0i=1\\ {\dot{H}}_{\mathrm{s},\mathrm{\text{out}},i-1}\mathit{\text{else}}\end{array}$$(11)

Furthermore, heat can be transported via convection from the bulk flow into the soil as according to(12) $${\dot{Q}}_i=h_{\vartheta ,i}{A}_{\mathrm{\text{int}},i}({\vartheta }_{\mathrm{f}}-{\vartheta }_{\mathrm{s},i})$$(12) where $h_{\vartheta ,i}$ is the heat transfer coefficient of segment $i\text{.}$ It must be determined individually for the problem at hand. The temperature of the bulk flow, ${\vartheta }_f\text{,}$ is assumed to be constant, since the mass of soil acting as a heat sink or source is small. The enthalpy is updated using(13) $$H_{\mathrm{s},i}^{\left(n\right)}={\tilde{H}}_{\mathrm{s},i}^{\left(n\right)}+\mathrm{\Delta }t{\dot{Q}}_i^{\left(n\right)}$$(13)

Computational algorithm

Assuming that hydrodynamic loads ${\tau }_{\mathit{\text{hyd}},i}\text{,}$ initial masses $m_{s,i}^0$ and enthalpies $H_{s,i}^0\text{,}$ are known, the following steps are computed within each time step:

1. Compute soil heights $h_{\mathrm{s},i}^{(n)}$ and temperatures ${\vartheta }_{\mathrm{s},i}^n$ from masses and energies.

2. Compute the shear rates ${\dot{\gamma }}_{\mathrm{s},i}^{(n)}$ (EquationEq. (2)(2) $${\tau }_{\mathrm{s}}=f({\dot{\gamma }}_{\mathrm{s}},{\vartheta }_{\mathrm{s}})$$(2) ) and the velocity profiles $u_{\mathrm{s},i}^{(n)}(y)$ (EquationEq. (3)(3) $$u_{\mathrm{s},i}\left(y_{\mathrm{s}}\right)={\dot{\gamma }}_{\mathrm{s},i}{y}_{\mathrm{s}}$$(3) ).

3. Propagate the mass to the next timestep using EquationEqs. (4)(4) $${\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i}={\rho }_{\mathrm{s}}{\int }_{A_{\mathrm{\text{out}},i}}^{}{u}_{\mathrm{s},i}\mathrm{d}A$$(4) , Equation(5)(5) $${\dot{m}}_{\mathrm{s},\mathrm{\text{in}},i}=\{\begin{array}{l}0i=1,\\ {\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i-1}\mathit{\text{else}}\end{array}$$(5) and Equation(7)(7) $$m_{\mathrm{s},i}^{(n+1)}=m_{\mathrm{s},i}^{\left(n\right)}-\mathrm{\Delta }t\left({\dot{m}}_{\mathrm{s},\mathrm{\text{out}},i}^{\left(n\right)}-{\dot{m}}_{\mathrm{s},\mathrm{\text{in}},i}^{\left(n\right)}\right)$$(7) .

4. Propagate the enthalpy to the next timestep using EquationEqs. (9–11).

5. Account for thermal energy transported by convection from the bulk flow into the soil using EquationEqs. (12)(12) $${\dot{Q}}_i=h_{\vartheta ,i}{A}_{\mathrm{\text{int}},i}({\vartheta }_{\mathrm{f}}-{\vartheta }_{\mathrm{s},i})$$(12) and Equation(13)(13) $$H_{\mathrm{s},i}^{\left(n\right)}={\tilde{H}}_{\mathrm{s},i}^{\left(n\right)}+\mathrm{\Delta }t{\dot{Q}}_i^{\left(n\right)}$$(13) .

If the investigated case is isothermal, the shear rates ${\dot{\gamma }}_{s,i}^{(n)}$ do not vary over time. Consequently, the shear rate only needs to be computed once at the beginning. Additionally, in an isothermal case, steps four and five are not necessary. The cleaning model was implemented in Python 3.9.

Setup

As a first validation, flushing processes of chocolates in a straight pipe under isothermal conditions at temperatures of $\vartheta \in \left\{40,52\right\}°C$ were investigated. The setup was taken from Liebmann etal. [Citation21]. In the reference, simulations of the flushing process using the volume of fluid method from the CFD library OpenFOAM were performed. The setup is depicted in . It features a pipe of length $L=1\hspace{0.17em}\mathrm{m}$ and radius $R=13\hspace{0.17em}\mathrm{\text{mm,}}$ which is initially filled with white chocolate. During the flushing process, dark chocolate pushes out the white chocolate. Applying the present wording, white chocolate corresponds to the soil, and dark chocolate to the cleaning fluid. The rheology of both fluids is described using the Windhab model [Citation43], reading(14) $${\tau }_{\mathrm{s}}\left(\dot{\gamma }\right)={\tau }_0+\left({\tau }_1-{\tau }_0\right)\left(1-\mathrm{\text{exp}}\hspace{0.17em}\left(-\frac{\dot{\gamma }}{{\dot{\gamma }}_{\mathrm{\text{ref}}}}\right)\right)+{\eta }_{\infty }\dot{\gamma }$$(14)

The rheological parameters ${\tau }_0,{\tau }_1,{\dot{\gamma }}_{\mathit{\text{ref}}},{\eta }_{\infty }$ are temperature dependent and their modeling is described in the appendix or in [Citation21].

The density of the white chocolate was determined to ${\rho }_s=1200\hspace{0.17em}\mathrm{\text{kg}}/{\mathrm{m}}^3$ and the density of the dark chocolate to ${\rho }_f=1260\hspace{0.17em}\mathrm{\text{kg}}/{\mathrm{m}}^3\text{.}$ In the reference simulation [Citation21], flushing fluid entered the pipe at the inlet with a developed velocity profile and a bulk velocity of $u_b=0.1\mathrm{m}/\mathrm{s}\text{.}$

In the reference, the whole flushing process was investigated. However, the present model can only be applied when the soil is film-like and the BCCM assumptions are valid. Hence, it was required (i) to identify a start time $t_{\mathit{\text{st}}}$ when this is the case, and (ii) to sample the soil height at $t_{\mathit{\text{st}}}$ to obtain suitable initial conditions for the cleaning simulation. The height of the soil layer $h_s$ was determined from the OpenFOAM v7 simulation by evaluating(15) $$h_{\mathrm{s}}\left(x,t\right)=R-\sqrt{\frac{1}{\pi }{\int }_{A_{\mathrm{\text{av}}}\left(x\right)}^{}\alpha \left(x,y,z\right)\mathrm{d}A}$$(15) where $\alpha$ is the phase indicator being equal to one in the cleaning fluid and $0$ in the soil. The definition of the coordinate directions and $A_{\mathit{\text{av}}}(x)$ can be found in . The second term in EquationEquation (15)(15) $$h_{\mathrm{s}}\left(x,t\right)=R-\sqrt{\frac{1}{\pi }{\int }_{A_{\mathrm{\text{av}}}\left(x\right)}^{}\alpha \left(x,y,z\right)\mathrm{d}A}$$(15) is the equivalent radius of the area occupied by the cleaning fluid [Citation21], and the integration was performed over cross-sections of the pipe, depicted in . A limit $h_{s,\mathrm{\text{max}}}$ was defined, for which the remaining layer of soil was sufficiently small to fulfill the assumptions of the present model. The first time, when $\forall x:h_s\left(x,t\right)\le h_{s,\mathrm{\text{max}}}$ was considered as the starting time $t_{\mathit{\text{st}}}$ for simulation using the present model. A value of $h_{s,\mathrm{\text{max}}}=0.1R$ was applied. Subsequently, the pipe was discretized in $N_{\mathit{\text{seg}}}$ equally sized segments of length $L_{\mathit{\text{seg}}}=L/N_{\mathit{\text{seg}}}\text{.}$ and in each segment, the soil height was averaged to obtain initial values of the soil thickness in each segment, $h_{s,i}^{(0)}\text{.}$

Finally, the hydrodynamic load ${\tau }_{\mathit{\text{hyd}}}$ was determined as input for the model. In the BCCM framework, the flow is calculated without consideration of the soil. In the case investigated, this resulted in computing the developed single-phase flow of a Windhab fluid in a straight pipe, where an analytical solution for the flow velocity and, hence, for the wall shear stress is available. This analytical solution, e.g., reported in the appendix, was used here to obtain the hydrodynamic load ${\tau }_{\mathit{\text{hyd}}}\text{.}$ Within the model framework, the domain occupied with soil corresponded to hollow cylinders where the areas depicted in are defined $A_{\mathit{\text{out}},i}=\pi \left(R^2-{\left(R-h_{s,i}\right)}^2\right)\text{,}$ $A_{\mathit{\text{int}},i}=2\pi \left(R-h_{s,i}\right)L_{\mathit{\text{seg}}}\text{.}$ The mass of each segment is computed as $m_{s,i}={{\rho }_{s}A}_{\mathit{\text{out}},i}{L}_{\mathit{\text{seg}}}\text{.}$ Cleaning simulations were performed with $N_{\mathit{\text{seg}}}=5$ and $\mathrm{\Delta }t=0.1\hspace{0.17em}\mathrm{s}\text{.}$

Results

displays the evolution of the soil height $h_{s,i}$ in the segments over time for ${\vartheta }_f={\vartheta }_s=40°C\text{.}$ The start time was determined to $t_{\mathit{\text{st}}}=16.2\mathrm{s}\text{.}$ The relative differences between the soil height from the present simulations and the reference [Citation21] are depicted in . The numbering of the segments is according to .

The model matches the OpenFOAM data in segments two, three, four and five. The finest cell near the wall used in the OpenFOAM simulations had a radial expansion of $\Delta r/R=0.01\text{.}$ Hence, the values below $h_{s,i}/R=0.01$ from the OpenFOAM simulation should be considered critically and cannot be employed for reasonable comparison. The reference simulation [Citation21] showed slight oscillations in segment five, indicating instabilities. Such instabilities are well-studied in the literature for core annular flows in straight pipes [Citation44–48]. Following [Citation48], there are three kinds of interfacial instabilities: (i) interfacial tension or capillary instabilities, which are dominant at very low Reynolds numbers, (ii) interfacial friction instabilities, caused by viscosity differences of the two fluids, dominant at low Reynolds numbers and (iii) Reynolds stress instabilities, caused by turbulence production in the bulk flow. Since surface tension was neglected and the Reynolds number is very small, these have to be interfacial friction instabilities. The cause of the instabilities may also be numerical. However, this has not yet been researched further and is beyond the scope of the present paper. The model at hand does not account for instabilities at all. However, the influence on the present results is negligible, and errors caused by this effect were minor. While the OpenFOAM simulation used as a reference case took around $20\hspace{0.17em}\mathrm{h}$ of computational time on 16 cores (AMD EPYC 7542 32-Core 2.9 GHz), the simulation with the present model took 0.2 s on four cores (Intel Core i5-6300, 2.4 GHz).

To obtain a single scalar value for the model performance, the relative root-mean-square error ${\epsilon }_{\mathit{\text{rms}}}$ was computed. This was done by comparing the height of the soil layer $h_{s,i}$ between the present simulation and the reference in each time step and averaging over the duration of the simulation, excluding values, where $h_{s,i,\mathit{\text{OF}}}/R<0.01\text{.}$ The error was normalized with the respective reference values $h_{s,i,\mathit{\text{OF}}}\text{.}$ Different timesteps $\mathrm{\Delta }t$ were investigated while holding the number of segments constant at $N_{\mathit{\text{seg}}}=5\text{.}$ The influence on the error was low if the stability criterion was fulfilled. The errors observed were around $6\%\text{.}$ The influence of the number of segments is shown in . The error is around $25\%$ for $N_{\mathit{\text{seg}}}=1$ and decreases toward a minimum value of $6\%$ at $N_{\mathit{\text{seg}}}=5\text{.}$ The error increases again for $N_{\mathit{\text{seg}}}>5\text{.}$ With a finer discretization, the local effect of interface instabilities on the soil height becomes more dominant.

The results obtained for $\vartheta =52°C$ are similar to $\vartheta =40°C$ and therefore provided in the appendix. While the minimum ${\epsilon }_{\mathit{\text{rms}}}$ error of $9\%$ is slightly higher compared to the $\vartheta =40°C$ case, the dependency of the error on the time step and the number of segments observed are the same.

Setup

A non-isothermal case from [Citation21] was investigated to validate the thermal modeling. The configuration is almost the same as described in the previous section. However, the temperature is not constant anymore and the walls are considered to be adiabatic. The soil was initially at $40°C\text{,}$ while the cleaning fluid entered the domain at $52°C\text{.}$ The thermal conductivity $\kappa$ of both fluids is assumed equal with $\kappa =0.3\hspace{0.17em}\mathrm{W}/(\mathrm{m}·\mathrm{K})$[Citation21]. Since the heat capacity varies only slightly with temperature, a constant value is assumed, with $c_{p,s}=1560\mathrm{J}/(\mathrm{\text{kg}}·\mathrm{K})\text{.}$

Additionally, the model required a value for the heat transfer coefficient $h_{\vartheta }\text{,}$ which was determined from the Nusselt number $\mathit{\text{Nu}}={(h_{\vartheta }L}_{\mathit{\text{chr}}})/\kappa \text{.}$ Chocolate is a highly viscous fluid, and the flow velocity was low here, which means the flow was laminar. Therefore, a Nusselt number of $\mathit{\text{Nu}}=4.36$ was assumed as an approximation, which is in the value for a laminar developed pipe flow of a Newtonian fluid with constant heat flux on the walls [Citation49]. Note that this approach does not account for developing a thermal boundary layer. This would be of particular importance in cases where an accurate prediction of a moving cleaning front is required. The characteristic length is $L_{\mathit{\text{chr}}}=2R$ in case of a pipe flow.

The initial condition for the soil height in each segment was determined the same way as in the previous section. However, an additional initial condition for the soil temperature was required here. It was obtained by averaging the temperature within the soil segments at $t_{\mathit{\text{st}}}\text{.}$ For the calculation of the hydrodynamic load ${\tau }_{\mathit{\text{hyd}}}\text{,}$ the analytical solution, reported in the appendix was used. To this end, it was assumed that the temperature of the flushing fluid remains constant at ${\vartheta }_f=52°C$ throughout the process.

Results

The simulations with the present model were repeated with the same discretization as above ($N_{\mathit{\text{seg}}}=5\text{,}$ $\mathrm{\Delta }t=0.1\hspace{0.17em}\mathrm{s}\text{,}$ $h_{s,\mathrm{\text{max}}}=0.1R$ and $t_{\mathit{\text{st}}}=19.0\mathrm{s}$) for the non-isothermal case. shows the evolution of soil height and soil temperature over time. Again, the present model matched the simulation results and the relative root-mean-square error ${\epsilon }_{\mathit{\text{rms}}}$ in terms of soil height is $6\%\text{,}$ while it is $1\%$ for the temperature. The present validation shows that the model can account for the temperature variation.

The simulation with the present model took 0.3 s on four cores (Intel Core i5-6300, 2.4 GHz), which is still negligible. It is worth noting that the thermal sub-model causes an increase of the duration by $50\%\text{,}$ which could become crucial in more complex cases. Step two of the algorithm listed above was the main reason for this increase. Adding the thermal sub-model made it necessary to recompute the shear rates ${\dot{\gamma }}_{s,i}^n$ using EquationEquation (2)(2) $${\tau }_{\mathrm{s}}=f({\dot{\gamma }}_{\mathrm{s}},{\vartheta }_{\mathrm{s}})$$(2) in each time step. In the present case, the Windhab model, which describes soil rheology, cannot be inverted analytically. Hence, a comparably expensive iterative procedure is required to invest EquationEquation (2)(2) $${\tau }_{\mathrm{s}}=f({\dot{\gamma }}_{\mathrm{s}},{\vartheta }_{\mathrm{s}})$$(2) .

Setup

As a third validation case, jet cleaning of a very thin oil layer is considered. The setup corresponds to the one presented in [Citation36]. A sketch of the setup with the definition of all relevant quantities is shown in . It features a vertical pipe with radius $r_0$ from which water with a temperature of $25°C$ is ejected to form a jet. The jet, which has a minimum radius $r_j<r_0\text{,}$ impinges perpendicularly on a flat plate. The flow can be characterized by the jet Reynolds number $R{e}_j=u_j{r}_j/{\nu }_j\text{,}$ which is computed with the maximum jet velocity $u_j\text{,}$ the minimum jet radius $r_j$ and the kinematic viscosity ${\nu }_j.$ A Newtonian oil layer covers the plate.

The time until the flow is established was negligible compared to the time required for cleaning. Hence, the development phase is neglected, and the flow is considered to be steady. In the region of impact, where $r<r_0\text{,}$ the pressure gradient and shear forces acting on the soil are of the same order of magnitude. The pressure drops rapidly with increasing radius and is negligible for $r>r_j$ [Citation50,Citation51]. Since the present modeling framework does not account for pressure forces, they are neglected entirely during the validation simulation. As in the previous validation simulation, an analytical solution of the single-phase flow is employed to obtain the hydrodynamic load ${\tau }_{\mathit{\text{hyd}}}\text{.}$ The analytical solution can be found in the appendix or in [Citation36,Citation52].

In the reference experiments [Citation36], the maximum jet velocity $u_j\text{,}$ the jet radius $r_0\text{,}$ the oil viscosity ${\eta }_s\text{,}$ the initial soil thickness $h_s^0\text{,}$ and the disk radius $R$ was varied. The cleaning fluid had a kinematic viscosity of the cleaning fluid, which is ${\nu }_j=8.9\cdot {10}^{-7}{\mathrm{m}}^2/\mathrm{s}\text{.}$ The quantity measured was the evolution of the soil height ${\overline{h}}_s$ over time. The overbar denotes an average over the plate. Furthermore, an analysis of the governing equations was performed in [Citation36] to derive a dimensionless time $t^{\ast }=0.089{\rho }_j{\nu }_j^{0.2}{u}_j^{1.8}{h}_s^{0}t/$(${\eta }_s{R}^{1.2})\text{.}$ Using this dimensionless time and rescaling the height using ${\overline{h}}_s^{\ast }={\overline{h}}_s/h_{s,0}$ results in a collapse of all measurements onto a single profile.

For validation, Run 5 of in [Citation36] was simulated, with $u_j=4.26\mathrm{m}/\mathrm{s}\text{,}$ $r_0=0.229\hspace{0.17em}\mathrm{\text{cm,}}$ ${\eta }_s=1\hspace{0.17em}\text{Pa s,}$ $h_s^0=h_{s,i}^0=60\hspace{0.17em}\mathrm{\mu }\mathrm{m}\text{,}$ $R=3.75\hspace{0.17em}\mathrm{\text{cm,}}$ $r_j=0.222\hspace{0.17em}\mathrm{\text{cm,}}$ $R{e}_j=10600$ using the present model. In that case $t^{\ast }=1$ corresponds to $t\approx 4\hspace{0.17em}\mathrm{s}\text{.}$ The geometry is described by $A_{\mathit{\text{out}},i}=2\pi i{L}_{\mathit{\text{seg}}}{h}_{s,i}\text{,}$ $A_{\mathit{\text{int}},i}=\pi L_{\mathit{\text{seg}}}^2(2i-1)\text{,}$ and $m_{s,i}={{\rho }_{s}A}_{\mathit{\text{int}},i}{h}_{s,i}\text{.}$ The length of a segment is $L_{\mathit{\text{seg}}}=R/N_{\mathit{\text{seg}}}\text{.}$

Results

Results obtained with the present model using $N_{\mathit{\text{seg}}}=5$ for the jet cleaning case are shown in and compared to experimental data and results obtained with data from the literature. A value of ${1/\overline{h}}_s^{\ast }=1$ means that the entire soil layer is present, ${1/\overline{h}}_s^{\ast }=100$ corresponds to $1\%$ of the initial soil being present.

The present model predicted almost the same results as the model in [Citation36]. For $t^{\ast }\le 4$ differences up to $40\%$ between the present model and the experiments were observed. The difference might be caused by soil being blasted off by the jet right at the footprint. For $t^{\ast }>4$ the present results agree better with the experimental data. The asymptotic behavior in the region $t^{\ast }>10$ is predicted almost perfectly. Apart from slight differences in the beginning, the present model predicts the same evolution of soil height over time as the model proposed in [Citation36]. These results validate, again, that the shear stress causes the removal of the thin soil layer. Since the exact data were not available, no further quantitative comparison was possible.