On the Response of a Herschel–Bulkley Fluid Due to a Moving Plate

In this paper, we study the boundary-layer flow of a Herschel–Bulkley fluid due to a moving plate; this problem has been experimentally investigated by others, where the fluid was assumed to be Carbopol, which has similar properties to cement. The computational fluid dynamics finite volume method from the open-source toolbox/library OpenFOAM is used on structured quad grids to solve the mass and the linear momentum conservation equations using the solver “overInterDyMFoam” customized with non-Newtonian viscosity libraries. The governing equations are solved numerically by using regularization methods in the context of the overset meshing technique. The results indicate that there is a good comparison between the experimental data and the simulations. The boundary layer thicknesses are predicted within the uncertainties of the measurements. The simulations indicate strong sensitivities to the rheological properties of the fluid.


Introduction
Constitutive modeling of complex fluids [1], sometimes referred to as non-linear fluids or non-Newtonian fluids, has received much attention in the literature [2][3][4]. Most of the naturally occurring and synthetic fluids are non-linear fluids, for example, polymer melts, suspensions, blood, slurries, drilling fluids, mud, etc., [5][6][7]. There are many empirical or semi-empirical constitutive equations that have been suggested for these fluids. Many nonlinear constitutive relations have also been derived based on the techniques of continuum mechanics [8][9][10][11]. The non-linearities oftentimes appear due to higher gradient terms or time derivatives.
Cement and concrete are among two of the most interesting complex materials. In fact, at least since the publication of a paper by Rivlin & Ericksen [12], who discussed fluids of complexity n (see also Truesdell & Noll [13]), to the recently published book [14], the term 'complex fluid' refers, in general, to fluid-like materials whose response, namely the stress tensor, is 'non-linear' in some fashion. This non-linearity can manifest itself in a variety of forms such as memory effects, yield stress, creep or relaxation, normalstress differences, etc., [15,16]. Cement has many applications, and it has been used in the oil and gas industries, where a cement slurry is pumped in the annulus space between the well casing and the geological formations surrounding the wellbore. This is carried out primarily to isolate the wellbore to keep fluids from migrating to other layers of the formation and secondly to prevent the corrosion and the eventual damage to the casing for the life of the well [17,18]. In time, the cement begins to harden. If the fluids from the surrounding formations penetrate the well, then disasters, both financial and operational, can occur, causing a shutdown and replacement of the cement. This unwanted phenomenon is known as 'gas migration' (see, e.g., [19]). In their powder forms, cement and concrete behave as bulk solids (granular materials); when mixed with water, initially they act as flowing suspensions (slurry) [20,21]; with chemical reactions and hydration occurring inside the suspension, the cement becomes a paste-like material (viscoelastic or viscoplastic) exhibiting yield stress and thixotropy [22,23]; and eventually, when it has hardened, the cement behaves as a poro-elastic material. Thus, cement can behave and can respond differently depending on the application and the conditions. When measuring its viscosity or its yield stress, cement generally behaves like a viscoplastic material.
In polymers, we can see a similar diversity. Agassant et al. [24] (p. xix) indicate that, in general, in polymer processing applications, one can distinguish three different stages: (1) the plastification (molten) stage where the polymer goes from a solid-like material, for example, powders or granular, to a fluid-like material, followed by (2) the molten polymer (see also [25]) being pushed or forced into molds or dies, and finally, (3) the stage where a final shape is given to the material, usually carried out via cooling. An excellent and early reference where the fundamentals of the modeling of these various stages in polymer processing are considered is the book by Middleman [26]. The Herschel-Bulkley (H-B) fluid model has been used in a variety of applications. Some recent applications are mentioned here. Ziaee at al. [27] studied colloidal-gas-aphron (CGA)-based fluids in drilling applications by modeling the fluid as a H-B fluid. In their study of solid-free polymer drilling fluid (SFPDF) with natural gas hydrates (NGH), Wang et al. [28], used the Herschel-Bulkley model. A new promising area for the application and use of polymeric gels seems to be in CO 2 underground storage where supercritical gas tends to leak through microcracks in wellbores (see [29]), where in some cases cement slurries, used in oilfields, are too vicious and are not able to penetrate the cracks. Chauhan et al. [30] looked at the characteristics of gum karaya suspensions as a fracturing fluid and developed an empirical Herschel-Bulkley model capable of predicting the temperature and concentration sensitivity of the apparent viscosity. Zheng et al. [31] looked at the effects of temperature and the rheological impact of a commonly used drilling fluid polymer-treating agent used in the petroleum industries; they mention that the dispersion was reasonably described by a Herschel-Bulkley model. Millian et al. [32] studied the rheological behavior of gel polymer electrolytes (GPE) used as a suspending fluid in a zinc-slurry-air RFB by fitting their experimental data to the Herschel-Bulkley model.
In his pioneering work on flow of yield stress fluids, Oldroyd [33] proposed a plastic boundary-layer theory, defined as a region of sufficiently slow plastic flow characterized by a large Oldroyd (Od or Bingham) number and evolving in the limit of small Reynolds numbers. Oldroyd further argued that the thickness of such a plastic boundary layer is of the order of Od −1/2 d (where d is a characteristic length). The proposed theory relied upon a certain number of assumptions, especially at the boundary between the elastic and the plastic states of the material, in addition to the assumption of a constant positive sign of the velocity gradient inside the boundary layer. When such a slow steady plastic flow develops near infinite or semi-finite thin plates, Oldroyd discussed a case of constant thickness, as well as a case of a variable layer with a thickness ranging from zero at the leading edge to a finite value far away from the edge. He also derived expressions for both the velocity and the pressure distributions inside the boundary layer, which depend on the Oldroyd number.
Piau [34] discussed Oldroyd's theory by pointing out certain inconsistencies buried in the approach; he mentioned five points. Essentially, these points can be summarized as the Dirichlet and the Neumann velocity boundary condition issues at the outer limit of the boundary layer where the transition occurs from a flowing material (a fluid) to an elastic material and at the wall. In addition, Piau pointed out that the assumptions on the pressure gradient and the symmetry conditions were not satisfied. Balmforth et al. [35] also claim that "Oldroyd's analysis runs into difficulties when the boundary layer buffers a wall, being unable to satisfy all the boundary conditions and the continuity equation". Piau revisited the theory and, in contrast to Oldroyd's approach, the Bingham stress of the material's plastic behavior was supplemented with the Hooke model for the linear elastic behavior that prevails in the outer unyielded regions. Piau [34] further identified and derived constant and variable lens-shaped boundary layer thicknesses and velocity distributions; he looked at the lower and the upper bound solutions consistent with the outer elastic region. The (viscous) drag forces acting on the plate in the context of these solutions, as derived by Piau, increase linearly with the yield stress and exhibits a relatively weak dependence on the Oldroyd (or the Bingham) number, which was assumed to be 'large' while the dependence of the drag forces was proportional to Od −1/2 .
Piau & Debiane [36] extended this work to shear-thinning fluids in the context of a Herschel-Bulkley fluid with no-slip conditions at the walls. They showed that the boundary thickness along with the velocity distribution, as well as the (viscous) drag force acting on the plate, are explicitly functions of the power-law exponent. For instance, in the framework of the constant thickness model of the boundary layer, the slope of the velocity distribution is found to be determined by the power-law exponent. For the case of slip at the walls, Piau & Debiane [36] introduced a dimensionless number (ratio of the yield stress to the consistency index and the fluid velocity), which is a measure of the slip at the wall. In the context of the constant-boundary-thickness model, they found that slip at the wall reduces the viscous drag, while the fluid velocity inside the boundary layer increases with decreasing slip. Ahonguio et al. [37] experimentally investigated the influence of slip at the wall in the limit of non-inertial flow with relatively large Oldroyd numbers (varying from 16 to 40). These authors found that the slip velocity decreases with the Oldroyd number. The consequences are (i) a thinner boundary layer and (ii) a reduction in the drag, which is consistent with the slip at the wall described by Piau and Debiane [36], although the definition of the boundary layer thickness in Ahonguio et al. [37] is different from that used by Piau & Debiane. In another work, Ahonguio et al. [38] showed that their laboratory measurements of the drag coefficient for a Carbopol gel flowing past a thin fixed plate compared well with the Piau & Debiane model. Balmforth et al. [39] revisited the derivation of Oldroyd's theory by numerically studying flow past a thin plate. These authors also investigated a jet-like intrusion; these two examples, referred to as Oldroyd's canonical problems, were meant to illustrate their approach. Most importantly, Balmforth et al. discussed that the magnitude of the small parameter ( ), associated with the re-scaling of the flow in the normal direction that sets both the thickness of the boundary layer and the angular velocity of the rotating plug, must be Bi −1/2 (Bi being the Bingham number) in order to match the pressure within the viscoplastic boundary layer. These authors claim that there is "a missing ingredient in Piau's boundary-layer scaling argument".
In the context of a boundary layer developing away from rigid boundaries (referred to "remote boundary layers" by Oldroyd [33]), Chevalier et al. [40] discussed experiments of slow injections of a yield stress fluid into a stagnant fluid (the same fluid), by means of an extrusion syringe. They observed that the injected fluid penetrates as a solid-like block over the whole injection surface, while the large surrounding material remains at rest. Their study shows the existence of a thin boundary through which the motion of the injected material occurs. They also reported a decrease in that boundary layer with an increasing of the Bingham number. Chevalier et al. also looked at the data from Boujlel et al. [1] where a plate was slowly immersed into a bath of Carbopol gel at rest. As the fluid was stressed beyond the yield point, a thin layer around the plate developed. Boujlel et al. showed the measurements of the boundary layer size for various immersion velocities as well as the distribution of the velocity within the yielded region.
The response of a yield-stress fluid, for example, a Herschel-Bulkley fluid [41], due to the motion of a plate can provide useful information about the resistance (drag) to flow and how the plate can move in the yielded regions as opposed to the viscous regions of the flow. In a sense, this flow arrangement can be thought of as an idealization of the slow movement of a vane in a viscometer while measuring the yield stress of cement or a yield stress fluid [42]. Carbopol is a fluid which has been studied extensively and has similar properties to cement. In this paper, we present numerical solutions to the boundary-layer flow of a Herschel-Bulkley fluid showing a solid-like behavior away from the boundaries, which were reported by Boujlel et al. [1]. In Section 2, we provide a description of the problem while briefly mentioning the experimental investigation of Boujlel et al. [1], which is relevant to our work. This is followed by an overview of the mathematical model (in Sections 3 and 4) used to describe the fluid and the flow conditions. The numerical method is presented in Section 5, and the results which are obtained using regularization approach in the context of the overset meshing technique are compared and discussed against the measurements of Boujlel et al. [1]. Finally, some conclusions and interesting points for future work are provided.

Problem Statement
The experiments of Boujlel et al. [1] of a plate slowly being immersed in a yield stress fluid are numerically investigated here. The yield stress fluid, which is at rest in a parallelpiped-shaped container 10 cm wide, 25 cm high and 16 cm deep, is a solution of Carbopol in water with a concentration of 0.5%. The fluid is assumed to behave as a Herschel-Bulkley fluid, and its rheological properties are obtained from fitting of the rheometrical flow curves: the yield stress (τ 0 ), the consistency (k), and the power-law exponent (n) are approximated as 59.5 Pa, 23.6 Pa.s n , and 0.38, respectively, for shear rates ranging between 10 −2 s −1 and 10 2 s −1 . The plate is 25 cm long (l), 7 cm wide (w) with a thickness (d) 1.5 mm. The flow domain along with the immersed plate are sketched in Figure 1. which were reported by Boujlel et al. [1]. In Section 2, we provide a description of the problem while briefly mentioning the experimental investigation of Boujlel et al. [1], which is relevant to our work. This is followed by an overview of the mathematical model (in Sections 3 and 4) used to describe the fluid and the flow conditions. The numerical method is presented in Section 5, and the results which are obtained using regularization approach in the context of the overset meshing technique are compared and discussed against the measurements of Boujlel et al. [1]. Finally, some conclusions and interesting points for future work are provided.

Problem Statement
The experiments of Boujlel et al. [1] of a plate slowly being immersed in a yield stress fluid are numerically investigated here. The yield stress fluid, which is at rest in a parallelpiped-shaped container 10 cm wide, 25 cm high and 16 cm deep, is a solution of Carbopol in water with a concentration of 0.5%. The fluid is assumed to behave as a Herschel-Bulkley fluid, and its rheological properties are obtained from fitting of the rheometrical flow curves: the yield stress ( 0 ), the consistency ( ), and the power-law exponent ( ) are approximated as 59.5 Pa, 23.6 Pa.s n , and 0.38, respectively, for shear rates ranging between 10 −2 s −1 and 10 2 s −1 . The plate is 25 cm long ( ), 7 cm wide ( ) with a thickness ( ) 1.5 mm. The flow domain along with the immersed plate are sketched in Figure 1. The flow conditions simulated in this work are summarized in Table 1. A total of three cases with different plate immersion velocity ( ) are considered. As shown in Table  1, the flow conditions correspond to very small Reynolds numbers, usually associated with the Stokes regime ( ≪ 1). In the limit of such (creeping) flow conditions where the inertial effects are negligible, the flow depends on the Bingham number, which for these experimental conditions, indicates yield stress effects which dominate the viscous ones (by about 2 to 3 times). While the limit of small Reynolds numbers is attained, one may notice that the Bingham numbers are not as large as they are prescribed in the theories (see [33,34,36,39]). The Reynolds ( ) and the Bingham ( ) numbers are defined below as (see the dimensionless form of the equation): where 0 the yield stress, the consistency, the power-law exponent, the density of the fluid, U p is a reference velocity, and is a reference length. In their experiments, Boujlel et al. defined an observation window (of 5 cm × 6.5 cm) 5 cm below the Carbopol bath surface, where successive pictures were taken once the leading edge of the plate appeared in the window until it was immersed to a depth of 20 cm. The average velocity profile of the Carbopol as well as the boundary layer thickness discussed in the result The flow conditions simulated in this work are summarized in Table 1. A total of three cases with different plate immersion velocity (U p ) are considered. As shown in Table 1, the flow conditions correspond to very small Reynolds numbers, usually associated with the Stokes regime (Re 1). In the limit of such (creeping) flow conditions where the inertial effects are negligible, the flow depends on the Bingham number, which for these experimental conditions, indicates yield stress effects which dominate the viscous ones (by about 2 to 3 times). While the limit of small Reynolds numbers is attained, one may notice that the Bingham numbers are not as large as they are prescribed in the theories (see [33,34,36,39]). The Reynolds (Re) and the Bingham (Bi) numbers are defined below as (see the dimensionless form of the equation): where τ 0 is the yield stress, k the consistency, n the power-law exponent, ρ the density of the fluid, U p is a reference velocity, and d is a reference length. In their experiments, Boujlel et al. defined an observation window (of 5 cm × 6.5 cm) 5 cm below the Carbopol bath surface, where successive pictures were taken once the leading edge of the plate appeared in the window until it was immersed to a depth of 20 cm. The average velocity profile of the Carbopol as well as the boundary layer thickness discussed in the result sections below are extracted from these measurements at 15 cm above the leading edge of the plate. The boundary layer thickness is estimated as the flow region extent delineated in its upper bound by the constant fluid velocity. As mentioned in the introduction, Piau & Debiane [36] showed that the boundary layer thickness (δ) can, in the context of the Herschel-Bulkley fluid, be approximated by a function which depends on the Oldroyd (or the Bingham) number:

Governing Equations
In this problem, we do not consider thermo-chemical or electromagnetic effects. Therefore, the governing equations of motion for a single component fluid include the conservation equations for mass, linear momentum, and angular momentum (see, e.g., [43]).

Conservation of Mass
where ∂/∂t is the partial derivative with respect to time, div is the divergence operator, v is the velocity vector, and ρ is the density of the fluid. If the fluid is assumed to be incompressible, then it can only undergo isochoric (i.e., volume preserving) motions, so: div v = 0 (5)

Conservation of Linear Momentum
where d/dt is the total time derivative given by d(.)/dt = ∂(.)/∂t + [grad(.)]v and grad is the gradient operator, b is the body force vector, and T is the Cauchy stress tensor.

Conservation of Angular Momentum
The conservation of the angular momentum indicates that the stress tensor is symmetric when there are no couple stresses, that is: Looking at the above equations, we can see that before we can solve any problems, we need a constitutive relation for T. In the next section, we provide a brief discussion of the stress tensor T used in this paper.
These conservation equations are supplemented with boundary and initial conditions, both at the walls and at the free-surface boundaries. The no-slip BC is prescribed at the plate and the walls such that: v = U p , at the plate boundaries (8) v = 0, at the container s walls (9) Polymers 2022, 14, 3890 where U p is the plate immersion velocity. The Neumann zero gradient condition is imposed at the free surface for the velocity vector. For the initial conditions, since the fluid is initially at rest, we use: v(x, 0) = 0 (10)

Constitutive Relation for the Stress Tensor
For many complex fluids, yield stress is an important rheological parameter [44][45][46][47]. In the oil and gas industries, predicting the yield stress for cement slurries is also important [48]. Here, a cement slurry is pumped in the well and then it begins to hydrate quickly and develop strength [49]. The difficulties related to yield stress measurements have been discussed, for example in [50][51][52][53]. Coussot [53] and Coussot et al. [54] reviewed different methods for measuring the yield stress for thixotropic non-Newtonian fluids. Experimental measurements for the yield stress are usually conducted either by direct rheometric techniques or indirect techniques. For a detailed discussion related to cement applications, see the report by Tao et al. [55]. One of the disadvantages of the direct techniques is the wall-slip effects, which cause under-estimation of the yield stress [56,57]. One of the most widely used techniques to measure the yield stress is the vane method since there is no wall slip during the shearing process within the material [56,[58][59][60]. Using the vane method [58,61], we can measure the peak torque-time response by rotating the vane immersed in the fluid. As mentioned in the Introduction section, the motion of a vane in the paste/suspension is of interest to us here.
In general, it can be assumed that the (Cauchy) stress tensor T for yield stress fluids, such as cement, can be defined as where T y is the yield stress tensor and T v is the viscous stress tensor. In general, for cement, the yield stress can be a function of many parameters, such as the volume fraction, w/c, etc.
where φ is the volume fraction, and w/c is the water-to-cement ratio. In a recent review article, Tao et al. [62], proposed a very general constitutive relationship for T v : where the kinematical tensors A 1 and A 2 are defined as: where p is the pressure, λ(t) is the structural parameter describing the degree of flocculation or aggregation. They used Krieger's idea [63] for the volume fraction dependence of the viscosity, where µ 0 is the (reference) coefficient of viscosity, tr is the trace operator, and m is the power law exponent, a measure of non-linearity of the fluid related to the shear-thinning effects (when m < 0) or shear-thickening effects (when m > 0). This model potentially is capable of exhibiting normal stress effects through the terms α 1 and α 2 , thixotropy effects because of the presence of the structural parameter λ, shear-rate-dependent effects of the viscosity through the two parameters α and m (showing shear-thinning or shear-thickening effects), and the concentration dependency of viscosity through the two parameters φ m and β. A simplified version Equation (13) was used in our earlier study [62]. For the yield stress part, historically, Oldroyd [33] derived a proper (frame invariant) 3D form for the Bingham fluid [64] by assuming that the material behaves as a linear elastic solid below the yield stress; he used the von Mises criterion for the yield surface. Thus: where G is the shear modulus, indicating that below the yield stress, the material behaves as a linear elastic solid, obeying the Hooke's Law, and where E is the strain tensor and the second invariant of the tensor A 1 is: As Denn [65] indicates, if the material is assumed to be inelastic prior to yielding, then G → ∞ , and Equation (18) is replaced by: Macosko [66] (p. 96) mentions that for many fluids with a yield stress, there is a lower Newtonian regime rather than a Hookean one, and thus one can use a two-viscosity (bi-viscous) model, such as: where . γ c is the critical shear rate. In this paper, we use the Herschel-Bulkley model, and thixotropy is not considered. Thus, the stress in the fluid is given by: where p is the pressure (the mean value of the stress tensor), I is the identity tensor, and τ is the stress tensor: in which τ 0 is the yield stress, k is the consistency index, and n is the power-law exponent, which measures of non-linearity of the fluid and is related to the shear-thinning effects (when n < 1) or shear-thickening effects (when n > 1). II τ and II A 1 are the second invariants of the stress tensor and of the kinematical tensor A 1 . In Equation (24), the total contribution in the brackets, which defines the viscosity of the fluid, is the sum of the shear (viscous) (µ v = k|II A 1 | (n−1)/2 ) and the apparent (µ app = τ 0 /|II A 1 | 1/2 ) viscosities.
In this paper, we ignore the micro-structure of the cement, i.e., the size and the shape of the particles, and the impact of the porosity and how the volume fraction affects the motion of the fluid. Thus, we represent the cement suspension as a viscoplastic fluid modeled as a Herschel-Bulkley fluid, given by Equations (23)-(25).

Numerical Approach
The computational fluid dynamics finite volume method from the open-source toolbox/library, OpenFOAM [67], is used on structured quad grids to solve the mass and the linear momentum conservation equations, using the solver "overInterDyMFoam" customized with non-Newtonian viscosity libraries. In our numerical scheme, we use the regularization methods. Indeed, to avoid the numerical implementation challenges associated with the discontinuity (singularity) in the stress tensor field between the unyielded and the yielded regions, the regularization method is used where the stress tensor τ is replaced with an -dependent small parameter such that: where the -dependent viscosity η is approximated in this work according to Papanastasiou [68] by: Two other commonly used regularization methods are also employed to study the sensitivity of the solution to such viscosity approximations; these two methods are the "simple" algebraic procedure (see e.g., [69]) and the approximation suggested by Bercovier & Engelman [70], given below, respectively: A detailed examination of the convergence challenges and issues associated with the regularized solutions are discussed, for example, in Frigaard & Nouar [71] and Saramito & Wachs [72].
Substituting Equations (23)- (25) in Equation (6), we have the basic equations, which need to be solved numerically: div And the boundary conditions are, • at the moving plate: • at the container's walls: and the initial conditions are: v(x, 0) = 0, p(x, 0) = 0 The dimensionless forms of the equations are presented in Appendix A.
The grid of the computational domain is generated relying upon the overset mesh technique, which in this work consists of a uniform grid (in each direction) of the background mesh of the entire flow domain, supplemented by a fine grid around the downward moving plate. With this fine (or overset) mesh, the grid in the normal direction to the plate is clustered using nonuniform spacings according to a geometric series with a rational stretching factor; this would better capture the gradients as the plate moves through the fluid. Figure 2 shows the background mesh, as well as the overset mesh which covers a region that extends over 3 cm at either side of the plate. Shown in Table 2 is the summary of the overset grid densities used to study the sensitivity of the solution to the mesh refinement.
• grad p = 0 (35) and the initial conditions are: The dimensionless forms of the equations are presented in Appendix A.
The grid of the computational domain is generated relying upon the overset mesh technique, which in this work consists of a uniform grid (in each direction) of the background mesh of the entire flow domain, supplemented by a fine grid around the downward moving plate. With this fine (or overset) mesh, the grid in the normal direction to the plate is clustered using nonuniform spacings according to a geometric series with a rational stretching factor; this would better capture the gradients as the plate moves through the fluid. Figure 2 shows the background mesh, as well as the overset mesh which covers a region that extends over 3 cm at either side of the plate. Shown in Table 2 is the summary of the overset grid densities used to study the sensitivity of the solution to the mesh refinement.  The convective term in the momentum equations is discretized using the second order "linear" scheme. Spatial gradients are also discretized using the second order "linear" scheme (central differences with linear interpolation). The simulations are performed while discretizing the unsteady terms with a backward Eulerian scheme. The coupling  The convective term in the momentum equations is discretized using the second order "linear" scheme. Spatial gradients are also discretized using the second order "linear" scheme (central differences with linear interpolation). The simulations are performed while discretizing the unsteady terms with a backward Eulerian scheme. The coupling between the background and the overset grids at these mesh boundaries is accomplished for the resolved fields (v, p) through a cell-volume-weighted interpolation scheme.
In the next section we discuss the results of our numerical simulations. Figures 3a, 4a and 5a show the contour plots of the instantaneous (vertical) velocity of the fluid as the plate is being immersed. As seen, the fluid mainly moves within a narrow area around the plate. There are two recirculation zones below the plate's leading edge, and the upward displacements at either side of the plate are readily apparent beyond this narrow area in the vicinity of the plate, where the fluid is dragged down by the plate (see Figures 3b, 4b and 5b). This flow pattern, which is consistent with the experimental observations, remains unchanged with an increase in the immersion velocity.

Flow Visualization
Error! Reference source not found.a, 4a and 5a show the contour plots of the instantaneous (vertical) velocity of the fluid as the plate is being immersed. As seen, the fluid mainly moves within a narrow area around the plate. There are two recirculation zones below the plate's leading edge, and the upward displacements at either side of the plate are readily apparent beyond this narrow area in the vicinity of the plate, where the fluid is dragged down by the plate (see Figures 3b, 4b and 5b). This flow pattern, which is consistent with the experimental observations, remains unchanged with an increase in the immersion velocity.  The square root of the second invariant of the viscous tensor (i.e., |II τ | 1/2 ) plotted in Figures 3c, 4c and 5c, is compared against the yield stress τ 0 . It appears that the larger values of the second invariants of the viscous stress tensor are localized just around the plate. The fluid region, which is identified as the region where the second invariant exceeds the yield stress τ 0 , following the von Mises yield criterion, exhibits an anchor-like shape around the plate. This indicates that, in addition to the material, which behaves as a fluid within a small envelope surrounding the plate, there exists also a fluid-like region attached to the leading edge of the plate which extends at either side of the leading edge of the plate. This could have originated from the material deformation caused and sustained by the continuous penetration of the leading edge of the plate. Boujlel et al. [1] reported that although some material was liquefied just below the leading edge of the plate, most of the material was liquefied around the leading edge. Furthermore, two narrow fluid regions also appear at the walls of the container.      [73] reported from their experimental observations that such a transition is characterized by a competition between destruction and reformation of the gel.
The boundary-layer thickness, measured as the extent of the yielded region adjacent The boundary-layer thickness, measured as the extent of the yielded region adjacent to the plate, is plotted for the three velocities in Figure 7. The plot shows that the variation in the yielded region with the increase in the plate velocity is adequately predicted both in trend and quantitatively; the maximum relative error over these three conditions is only about 1.07%. The predictions are indeed within the uncertainties of the experimental data as shown by the error bars.
in the yielded region with the increase in the plate velocity is adequately predicted both in trend and quantitatively; the maximum relative error over these three conditions is only about 1.07%. The predictions are indeed within the uncertainties of the experimental data as shown by the error bars. We can see from Figure 7 that an increase in the Bingham number results in a decrease in the boundary layer thickness. This suggests that for creeping flows, where the yield stress dominates over the viscous stress, the boundary layer developing around the moving plate will exhibit a negative correlation dependence on the Bingham number. This is consistent with the boundary layer theory proposed by Piau & Debiane [36], which is an improvement and extension to the Herschel-Bulkley model of Oldroyd's theory [33]. Indeed, Piau & Debiane showed that the boundary layer thickness scales as , which with the current rheological properties results in −0.7246 (given that = 0.38). However, the experimental data (or at least the subset investigated in this work), indicate that the thickness scaled as −0.1213 .

Mesh and Regularization Sensitivities
The sensitivity of the solution to the mesh and the regularization method are examined below for the plate velocity of 1 mm/s. The sensitivity to the mesh resolution is performed using Papanastasiou's approximation (with = 10 −5 ). The mesh densities along with the resolutions are summarized in Table 2. Figure 8 compares the velocity profiles for the three grids against the measurements. It is apparent that the medium-grid solution does quite well with the fine-grid prediction, while the coarse-grid one slightly deviates We can see from Figure 7 that an increase in the Bingham number results in a decrease in the boundary layer thickness. This suggests that for creeping flows, where the yield stress dominates over the viscous stress, the boundary layer developing around the moving plate will exhibit a negative correlation dependence on the Bingham number. This is consistent with the boundary layer theory proposed by Piau & Debiane [36], which is an improvement and extension to the Herschel-Bulkley model of Oldroyd's theory [33]. Indeed, Piau & Debiane showed that the boundary layer thickness scales as Bi −1/(1+n) , which with the current rheological properties results in Bi −0.7246 (given that n = 0.38). However, the experimental data (or at least the subset investigated in this work), indicate that the thickness scaled as Bi −0.1213 .

Mesh and Regularization Sensitivities
The sensitivity of the solution to the mesh and the regularization method are examined below for the plate velocity of 1 mm/s. The sensitivity to the mesh resolution is performed using Papanastasiou's approximation (with = 10 −5 ). The mesh densities along with the resolutions are summarized in Table 2. Figure 8 compares the velocity profiles for the three grids against the measurements. It is apparent that the medium-grid solution does quite well with the fine-grid prediction, while the coarse-grid one slightly deviates from those two solutions. Furthermore, it can be seen from the fine-grid solution that the slight mismatch remains despite the refinement, especially in the vicinity of the plate. As summarized in Table 3, the predicted boundary layer thickness with respect to different grids are within the measurement uncertainties. The prediction with the medium grid compares well with the fine grid predictions and the slight mismatch in the velocity profile marginally reflects grid resolution issues.

Mesh and Regularization Sensitivities
The sensitivity of the solution to the mesh and the regularization method are examined below for the plate velocity of 1 mm/s. The sensitivity to the mesh resolution is performed using Papanastasiou's approximation (with = 10 −5 ). The mesh densities along with the resolutions are summarized in Table 2. Figure 8 compares the velocity profiles for the three grids against the measurements. It is apparent that the medium-grid solution does quite well with the fine-grid prediction, while the coarse-grid one slightly deviates from those two solutions. Furthermore, it can be seen from the fine-grid solution that the slight mismatch remains despite the refinement, especially in the vicinity of the plate. As summarized in Table 3, the predicted boundary layer thickness with respect to different grids are within the measurement uncertainties. The prediction with the medium grid compares well with the fine grid predictions and the slight mismatch in the velocity profile marginally reflects grid resolution issues.  The influence of the regularization method is shown in Figure 9. This plot shows that there is little difference between the Bercovier-Engelman and the Papanastasiou methods. The boundary layer thickness prediction compared in Table 3 shows that the difference is very small. As for the "simple" regularization approach, the velocity profile seems to be  The influence of the regularization method is shown in Figure 9. This plot shows that there is little difference between the Bercovier-Engelman and the Papanastasiou methods. The boundary layer thickness prediction compared in Table 3 shows that the difference is very small. As for the "simple" regularization approach, the velocity profile seems to be less accurately predicted; the same observation holds for the thickness of the boundary layer (see Table 3). less accurately predicted; the same observation holds for the thickness of the boundary layer (see Table 3). Figure 9. Sensitivity of the velocity distribution to the regularization model.

Results for Different Fluid Properties
Here, we look at the difference between the numerical predictions of the velocity inside the boundary layer obtained by prescribing the thickness of the boundary layer ( ), and the velocity distribution measured by Boujlel et al. (2012). This velocity distribution is a function of the boundary layer thickness ( ) (see [1,36]): where is the plate velocity and is the power-law exponent of the Herschel-Bulkley fluid model.
In Figure 10, we plot the velocity distribution using the measured thickness of = 9.01 mm and = 14.4 mm proposed by Boujlel et al. The plot prediction using the measured thickness exhibits a less satisfactory agreement contrary to the larger boundary layer thickness for which a very good agreement is found. Since the only difference is due to a different value of the boundary layer thickness, which depends on the Bingham number, we think it would be interesting to look at different values for the properties of the fluid and see the impact on the results. Although Boujlel et al. [1] carefully prepared and performed the measurements, they did not mention uncertainties in their measurements. Piau [74] argued from an extensive literature review that often the rheometry of Carbopol aqueous gels is not without measurement difficulties.

Results for Different Fluid Properties
Here, we look at the difference between the numerical predictions of the velocity inside the boundary layer obtained by prescribing the thickness of the boundary layer (δ), and the velocity distribution measured by Boujlel et al. (2012). This velocity distribution is a function of the boundary layer thickness (δ) (see [1,36]): where U p is the plate velocity and n is the power-law exponent of the Herschel-Bulkley fluid model. In Figure 10, we plot the velocity distribution using the measured thickness of δ = 9.01 mm and δ = 14.4 mm proposed by Boujlel et al. The plot prediction using the measured thickness exhibits a less satisfactory agreement contrary to the larger boundary layer thickness for which a very good agreement is found. Since the only difference is due to a different value of the boundary layer thickness, which depends on the Bingham number, we think it would be interesting to look at different values for the properties of the fluid and see the impact on the results. less accurately predicted; the same observation holds for the thickness of the boundary layer (see Table 3). Figure 9. Sensitivity of the velocity distribution to the regularization model.

Results for Different Fluid Properties
Here, we look at the difference between the numerical predictions of the velocity inside the boundary layer obtained by prescribing the thickness of the boundary layer ( ), and the velocity distribution measured by Boujlel et al. (2012). This velocity distribution is a function of the boundary layer thickness ( ) (see [1,36]): where is the plate velocity and is the power-law exponent of the Herschel-Bulkley fluid model.
In Figure 10, we plot the velocity distribution using the measured thickness of = 9.01 mm and = 14.4 mm proposed by Boujlel et al. The plot prediction using the measured thickness exhibits a less satisfactory agreement contrary to the larger boundary layer thickness for which a very good agreement is found. Since the only difference is due to a different value of the boundary layer thickness, which depends on the Bingham number, we think it would be interesting to look at different values for the properties of the fluid and see the impact on the results. Although Boujlel et al. [1] carefully prepared and performed the measurements, they did not mention uncertainties in their measurements. Piau [74] argued from an extensive literature review that often the rheometry of Carbopol aqueous gels is not without measurement difficulties. Although Boujlel et al. [1] carefully prepared and performed the measurements, they did not mention uncertainties in their measurements. Piau [74] argued from an extensive literature review that often the rheometry of Carbopol aqueous gels is not without measurement difficulties.
We now look at the boundary-layer flow for different properties of a yield-stress fluid modeled as a Herschel-Bulkley fluid. The results for the velocity profiles show a systematic slight mismatch of the predictions against the measurements in the yielded region behaving as a solid-like material. Table 4 lists the different rheological properties which we use in our simulations. In Figures 11-14, we plot the variations in the mean velocity (U), the kinematic viscosity (µ/ρ), and the mean shear stress (τ xy ) for different material properties.     Figure 14, shows that an increase in the consistency index results in a steeper profile, thus reducing the   Figure 14, shows that an increase in the consistency index results in a steeper profile, thus reducing the Figure 13. Influence of the power-law exponent on the viscous (left plot) and the apparent (right plot) viscosities for the plate immersion velocity of 1 mm/s. discrepancy. Typically, the change in the consistency index appears to show a similar change in the viscosity and the mean shear stress. This is illustrated by an increase in the viscosity and in the mean shear stress, especially in the wall region (see Figure 14b,c), causing a flatter profile as a result of the increased friction at the wall.

Conclusions
The boundary layer experiments of Boujlel et al. [1], where a plate is slowly immersed into a Carbopol yield-stress fluid, are numerically investigated. The fluid is modeled as a Herschel-Bulkley fluid. The numerical approach uses the Papanastasiou regularization of the apparent viscosity associated with the Herschel-Bulkley constitutive model in the context of the overset meshing technique. The physics of the flow appears to be adequately captured for the three flow conditions investigated in this work. It is seen that the boundary layer develops in the vicinity of the plate as the leading edge steadily advances into the fluid. The velocity distribution seems to be in good agreement with the measurements. Specifically, the upward velocity distributions of the unyielded material are reasonably  For smaller values of the yield stress, we see a better agreement with the experimental data. We also notice similar behavior for the viscosity and the mean shear stress shown in Figure 11b,c. The tendency of the fluid to resist shearing motion as conveyed by the viscosity, clearly increases with the yield stress as shown in Figure 11b. This plot further exhibits a steeper increase in the viscosity within the transition region between the fluid-like and the solid-like behaviors. Figure 11c also shows how the mean shear stress changes with the tangential shearing forces between the fluid and the plate, from near the plate to the far-field regions. These clearly indicate an increase in the friction at the plate with the yield stress, which result in slower velocities.
For different values of the power-law exponent (n), we could see similar tendencies (see Figure 12). However, the variations are not that accentuated as is the case for the yield stress discussed above. Only slight increases are apparent in both the viscosity and the mean shear stress (see Figure 12b,c), which result in a slight change in the velocity, especially for n = 0.25, where the mean shear stress is relatively higher compared to n = 0.38 and n = 0.50. The viscous and the apparent viscosities plotted in Figure 13 show that the apparent viscosity is one order of magnitude larger than the regular (shear) viscosity for each of the power-law exponents used in this study. Therefore, even though some shear-thinning effects takes place, the yield stress effects appear to be dominant. Piau & Debiane [36] and later Boujlel et al. [1] showed that the slope is determined by the power law exponent. In other words, the shear-thinning effects strongly influence the velocity field in the boundary layer.
The velocity distribution in the vicinity of the plate plotted in Figure 14, shows that an increase in the consistency index results in a steeper profile, thus reducing the discrepancy. Typically, the change in the consistency index appears to show a similar change in the viscosity and the mean shear stress. This is illustrated by an increase in the viscosity and in the mean shear stress, especially in the wall region (see Figure 14b,c), causing a flatter profile as a result of the increased friction at the wall.

Conclusions
The boundary layer experiments of Boujlel et al. [1], where a plate is slowly immersed into a Carbopol yield-stress fluid, are numerically investigated. The fluid is modeled as a Herschel-Bulkley fluid. The numerical approach uses the Papanastasiou regularization of the apparent viscosity associated with the Herschel-Bulkley constitutive model in the context of the overset meshing technique. The physics of the flow appears to be adequately captured for the three flow conditions investigated in this work. It is seen that the boundary layer develops in the vicinity of the plate as the leading edge steadily advances into the fluid. The velocity distribution seems to be in good agreement with the measurements. Specifically, the upward velocity distributions of the unyielded material are reasonably captured. However, a slight mismatch is noticeable within the yielded region near the solid-plug flow. It is shown that the transition between the fluid-like and the solid-like states, is very sensitive to the rheological properties of the fluid. Subsequent simulations show a strong dependence of the mean shear stress on these rheological properties that could cause different values for the friction at the wall. This behavior is consistent with the existing theories [36]. In our numerical investigations, we assume a no-slip condition at the walls. The influence of the slip at the walls (see [36,38]) is left for future work.