The constitutive model of Herschel-Bulkley fluids is similar to Bingham fluids or power-law fluids. The specific constitutive equation is given as follows: { τ = τ 0 + k γ ˙ n ,   | τ | > τ 0 γ ˙ = 0   ,   | τ | ≤ τ 0 (1) where τ is shear stress, k is viscosity coefficient, γ ˙ is the shear rate, n is the power-law index, and τ ₀ is the initial yield stress.

The spreading schematic of the advancing system of Newtonian or power-law fluids is shown in Figure 1. The derivation of the Herschel-Bulkley case was also started in Figure 1. In Figure 1, u is the horizontal velocity of the fluid, which is affected by time t and position x, z. U is the framing moving velocity of the contact line, θ(U) is the dynamic contact angle, and h is the height of film thickness of any position x.

With the proposed assumption (3), the pressure along the z-direction will meet the following equation. ∂ p ∂ z = − ρ g (2)

Applying the Young-Laplace equation to the free surface (z = h), there is p = p G + σ ( 1 R 1 + 1 R 2 ) ≈ p G - σ ∂ 2 h ∂ x 2 (3) where p _G represents the atmospheric pressure and σ is surface tension.

Integrating Eq. (2) with respect to the boundary condition listed in Eq. (3), the pressure can be further calculated as p = ρ g ( h − z ) + p G − σ ∂ 2 h ∂ x 2 (4)

When z = 0, the stress on the horizontal substrate can be obtained from Eq. (4) as τ s = p | z = 0 = ρ g h + p G − σ ∂ 2 h ∂ x 2 (5)

Consequently, when the initial yield stress τ ₀ is larger than τ _s , substantial flow will not occur. Thus the approximate condition for substantial flow is | τ 0 | < ρ g h + p G − σ ∂ 2 h ∂ x 2 (6)

Substantial flow will be discussed below. Owing to the initial yield stress, there must be a yield surface at z = h ₀. Thus the spreading zone will be divided into two zones, the significant shear exists only below the yield surface. The shear rate is zero at or above the yield surface. The viewpoint is consistent with the spreading of Bingham fluids proposed by Liu [20]. Therefore, the actual flow schematic of the advancing systems of Herschel-Bulkley fluids is shown in Figure 2. The dashed line separates the zone into sheared and yield zones respectively. The schematic figure is completely different from the Newtonian or power-law case.

For the Navier–Stokes equation, it can be written as ∂ p ∂ x = ∂ ∂ z ( μ ∂ u ∂ z ) (7) where p is the pressure and x is the distance of one point from the z-axis. According to Eq. (1), when the stress is larger than the yield initial stress, the viscosity can be expressed as μ = τ 0 | ∂ u / ∂ z | + k | ∂ u ∂ z | n − 1 (8)

Substituting Eq. (8) into Eq. (7), there is ∂ p ∂ x = ∂ ∂ z [ ( τ 0 | ∂ u / ∂ z | + k | ∂ u ∂ z | n − 1 ) ∂ u ∂ z ] (9)

Based on the boundary condition of no shear at the yield surface (h(x) = h ₀, ∂u/∂z = 0), integrating Eq. (9) with respect to z and the following equation can be acquired. ∂ u ∂ z = s i g n ( − ∂ p ∂ x ) [ | − 1 k ∂ p ∂ x | ( h 0 − z ) ] 1 n (10) With the boundary condition of no-slip at the solid surface (z = 0, u = 0), the expression of u can be obtained by integrating Eq. (10). u = n n + 1 ( 1 k ) 1 n s i g n ( − ∂ p ∂ x ) ⋅ | − ∂ p ∂ x | 1 n [ h 0 n + 1 n − ( h 0 − z ) n + 1 n ] (11)

Equation (11) shows that the maximum velocity u _p occurs at the yield surface z = h ₀ that can be obtained as u p = n n + 1 s i g n ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) ⋅ ( 1 k ) 1 n ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) 1 n h 0 n + 1 n (12)

When z > h ₀, the velocity is still equal to the maximum velocity u _p . According to mass conservation, there is ∂ h ∂ t + ∂ q ∂ x = 0 (13) where q is the flow of fluids.

For Herschel-Bulkley fluids, q is made up of two parts as q = ∫ 0 h 0 u d z + u p ( h − h 0 ) (14)

The first term of the right side represents the shear zone below the yield surface and its detailed expression is ∫ 0 h 0 u d z = n 2 n + 1 sign ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) ( 1 k ) 1 n ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) 1 n h 0 2 n + 1 n (15)

Then the flow q will be acquired as follows: q = n 2 n + 1 sign ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) ( 1 k ) 1 n ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) 1 n h 0 2 n + 1 n                               + n n + 1 sign ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) ( 1 k ) 1 n ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) 1 n h 0 n + 1 n ( h − h 0 ) (16)

In Fig. 2, the velocity along the x-direction is parabolic, which can also be explained by Eqs. (11)–(12). When it refers to the velocity, the case here is similar to classical Poiseuille flow. Both of them are driven by the pressure gradient, and they have the same initial and boundary conditions, thus the relationship between h and h ₀ can be obtained based on our previous work [21–22] as h 0 = h − τ 0 | ∂ p / ∂ z | = h − τ 0 ρ g (17)

Substituting Eq. (17) into Eq. (16), a new formula of flow q is generated as q = n 2 n + 1 sign ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) ( 1 k ) 1 n ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) 1 n h 0 2 n + 1 n                               + n n + 1 sign ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) ( 1 k ) 1 n ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) 1 n h 0 n + 1 n τ 0 ρ g (18)

Finally, substituting Eq. (18) into Eq. (13), the following equation is obtained as ∂ h ∂ t + n 2 n + 1 sign ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) ( 1 k ) 1 n ∂ ∂ x [ | σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x | 1 n h 0 2 n + 1 n ] + n n + 1 sign ( σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x ) ( 1 k ) 1 n ∂ ∂ x [ | σ ∂ 3 h ∂ x 3 − ρ g ∂ h ∂ x | 1 n h 0 n + 1 n τ 0 ρ g ] = 0 (19)

According to Eq. (16), the pressure mainly depends on the gravity and capillary forces. In the following part, two limiting regimes will be considered further. One is only the gravitational regime taken into consideration by ignoring the capillary force, ρgR ²>>σ. The other is only the capillary force regime considered, ignoring the gravity, σ>>ρgR ².

The capillary force is ignored here, so Eq. (19) can be simplified into the following equation. ∂ h ∂ t + sign ( − ρ g ∂ h ∂ x ) ( ρ g k ) 1 n ∂ ∂ x [ | − ∂ h ∂ x | 1 n h 0 n + 1 n ( n 2 n + 1 h 0 + n n + 1 τ 0 ρ g ) ] = 0 (20)

When the contact line moves at velocity U, the above equation can be converted into the following form by introducing a new variable ξ = x − U t . − U d h 0 d ξ + ( ρ g k ) 1 n d d ξ [ ( d h 0 d ξ ) 1 n h 0 n + 1 n ( n 2 n + 1 h 0 + n n + 1 τ 0 ρ g ) ] = 0 (21)

Integrating the above equation with respect to ξ, there is U n = ( n 2 n + 1 ) n ρ g k d h 0 d ξ h 0 ( h 0 + 2 n + 1 n + 1 τ 0 ρ g ) n (22) With the boundary condition: ξ = 0, h = 0, Eq. (22) can be solved as 1 n + 2 [ ( h + n n + 1 τ 0 ρ g ) n + 2 − ( n n + 1 τ 0 ρ g ) n + 2 ] − 2 n + 1 ( n + 1 ) 2 ⋅ τ 0 ρ g [ ( h + n n + 1 τ 0 ρ g ) n + 1 − ( n n + 1 τ 0 ρ g ) n + 1 ] = ( U 2 n + 1 n ) n k ρ g ξ (23)

Equation (23) addresses the fact that film thickness in capillary spreading is determined by initial yield stress τ ₀ and power index n.

In the condition, the gravitational action is ignored, then Eq. (19) can be simplified as follows ∂ h ∂ t + sign ( σ ∂ 3 h ∂ x 3 ) ( 1 k ) 1 n ∂ ∂ x { ( σ ∂ 3 h ∂ x 3 ) 1 n h 0 n + 1 n [ n 2 n + 1 h 0 + n n + 1 τ 0 ρ g ] } = 0 (24)

Introducing the above variable ξ and h ₀ into Eq. (24), − U d h 0 d ξ + sign ( σ d 3 h 0 d ξ 3 ) ( 1 k ) 1 n d d ξ { ( σ ∂ 3 h 0 ∂ ξ 3 ) 1 n h 0 n + 1 n [ n 2 n + 1 h 0 + n n + 1 τ 0 ρ g ] } = 0 (25)

Integrating the above equation with respect to ξ, − U + n 2 n + 1 ( 1 k ) 1 n ( σ d 3 h 0 d ξ 3 ) 1 n h 0 1 n [ h 0 + 2 n + 1 n + 1 τ 0 ρ g ] = 0 (26)

Then it is supposed that τ 0 ρ g = λ h 0 , so Eq. (26) can be translated into U n = [ n 2 n + 1 ⋅ ( 2 λ + 1 ) n + ( λ + 1 ) n + 1 ] n σ k d 3 h 0 d ξ 3 h 0 n + 1 (27) With the boundary condition: ξ = 0, h = 0, Eq. (27) can be solved as h 0 n + 2 = k σ [ U ( 2 n + 1 ) ( n + 1 ) n ( n ( 2 λ + 1 ) + ( λ + 1 ) ) ] n [ ( n + 2 ) 3 3 ( 2 n + 1 ) | n − 1 | ] ξ 3           ( n ≠ 1 ) (28)

Further, the film thickness equation can be obtained by replacing h _0, ( h − τ 0 ρ g ) n + 2 = k σ [ U ( 2 n + 1 ) ( n + 1 ) n ( n ( 2 λ + 1 ) + ( λ + 1 ) ) ] n [ ( n + 2 ) 3 3 ( 2 n + 1 ) | n − 1 | ] ξ 3               ( n ≠ 1 ) (29)

Based on Eq. (17), the variable λ can be calculated as λ = τ 0 / ρ g h − τ 0 / ρ g (30)

In this way, the equation of film thickness h can be obtained, which depends on the power index n and initial yield stress τ ₀.

The inclination angle at x = ξ can be calculated by differentiating film thickness h as shown in Eq. (29) for the capillary spreading regime. tan ⁡ θ = d h d ξ = k σ U n [ ( n + 2 ) 3 ( 2 n + 1 ) | n − 1 | ] ( h − τ 0 / ρ g ) [ n ( n + 1 ) h + n 2 τ 0 / ρ g ( 2 n + 1 ) ( n + 1 ) ] n − 1 { ( n + 2 ) [ n ( n + 1 ) h + n 2 τ 0 / ρ g ( 2 n + 1 ) ( n + 1 ) ] − n 2 τ 0 / ρ g n + 1 } ξ 2   ( n ≠ 1 ) (31)

In previous studies, many researchers found that the local microscopic contact angle can not be measured directly for boundary conditions. Thus dynamic contact angle θ is taken as the replacement, which is equal to the inclination angle at x = x _m . Furthermore, the dynamic contact angle can be obtained as θ ≈ tan ⁡ θ = k σ U n [ ( n + 2 ) 3 ( 2 n + 1 ) | n − 1 | ] ( h − τ 0 / ρ g ) [ n ( n + 1 ) h + n 2 τ 0 / ρ g ( 2 n + 1 ) ( n + 1 ) ] n − 1 { ( n + 2 ) [ n ( n + 1 ) h + n 2 τ 0 / ρ g ( 2 n + 1 ) ( n + 1 ) ] − n 2 τ 0 / ρ g n + 1 } x m 2     ( n ≠ 1 ) (32)

To obtain Eq. (32), it is considered that when the angle θ is small enough, it is approximately equal to sin θ or tan θ.