Gradient-enhanced visco-plastic models: Advancing localization and deformation analysis in toothpaste flow and hollow sphere deformation

Abstract

This study explores the incorporation of gradient-enhanced visco-plastic models to address the limitations of traditional visco-plastic formulations in capturing size-dependent and localized deformation behaviors. By introducing second-gradient terms, the proposed framework accounts for non-local effects and smoothens discontinuities in stress and deformation profiles, which are critical for accurately modeling materials with complex microstructures. Analytical solutions are derived for benchmark problems, including the deformation of a hollow sphere and toothpaste-like flow, demonstrating the enhanced predictive capabilities of the gradient-enhanced models. Numerical analyses further validate these formulations, highlighting their ability to stabilize computations and provide physically realistic stress distributions under highly localized loading conditions. The results reveal the significance of the gradient coefficient in influencing stress diffusion and plastic strain localization, emphasizing its role in material design and engineering applications. The study establishes robust theoretical and computational foundations for gradient-enhanced visco-plastic models, offering new insights into material behavior at micro and meso-scales, with implications for advanced manufacturing and processing technologies.

Keywords

Norton model numerical simulations gradient-enhanced model visco-plasticity strain localization

1. Introduction

The study of visco-elastic and visco-plastic materials has long been a focal point in mechanics, owing to their relevance across a wide array of applications, from geological processes to industrial material design. Visco-elasticity, describing materials that exhibit both viscosity and elasticity, and visco-plasticity, encompassing materials that yield under stress, offer essential frameworks for understanding and predicting the mechanical behavior of complex materials. In this work, we investigate the application of viscous models to two specific mechanical problems: the hydrostatic behavior of a hollow sphere and the flow dynamics of a material with toothpaste-like properties in a cylindrical tube.

1.1. Theoretical background and context

The analysis of materials under hydrostatic loading is a well-established area of study in continuum mechanics. Hollow spheres under hydrostatic stress have been extensively studied, as they provide simplified models for porous, granular, or damaged materials subject to compressive stress. Classic visco-plasticity theories (e.g., Hill [1], Drucker and Prager [2]) initially focused on perfectly plastic materials, while subsequent models integrated viscous effects to capture time-dependent deformation, such as the Norton model [3], which describes materials undergoing steady-state creep under constant stress. This model has been applied widely in geotechnics, metallurgy, and bioengineering to simulate long-term material behavior [4,5]. The Norton model, in particular, has been shown to accurately describe the visco-plastic flow of metals at high temperatures, where elastic effects can be negligible [6,7].

In the context of porous materials, the study of flow within a hollow sphere provides essential insight into the behavior of materials with microstructural voids [8]. Models focusing on porosity have advanced from initial elasticity-based approaches (e.g., Gurson [9]) to visco-plastic descriptions that account for both matrix flow and void growth at elevated temperatures [10]. Current research continues to examine these models under various boundary conditions, including both confined and unconfined compressive states.

1.2. The Norton viscous model

The Norton viscous model is a well-established framework for describing the time-dependent deformation of materials under stress, particularly in high-temperature creep and plasticity. It is especially useful for modeling steady-state deformation behavior in metals, polymers, and other viscoelastic materials. The model can be described mathematically as:

\overset{\cdot}{ϵ} = A σ^{n} \exp (- \frac{Q}{RT}),

where $\overset{\cdot}{ϵ}$ is the strain rate, σ is the applied stress, n is the stress exponent, A is a material constant, Q is the activation energy, R is the universal gas constant, and T is the absolute temperature. This equation captures the dependence of strain rate on stress and temperature, which is crucial for understanding high-temperature material behavior.

The model was initially proposed by Norton [11] to describe creep deformation in metals. Norton’s law states that the strain rate is proportional to the applied stress raised to a power n. While this simple power-law relationship adequately describes steady-state creep, subsequent studies incorporated the temperature dependence via an Arrhenius-type term to account for the thermal activation of deformation mechanisms [12].

One of the model’s critical features is its incorporation of thermal activation through the exponential term $\exp (- Q ∕ RT)$ , which reflects the reduction of creep resistance with increasing temperature. Orowan [12] demonstrated the relevance of thermal activation in creep processes, and this insight has been integral to the model’s success in predicting temperature-dependent deformation.

Several modifications have been made to the Norton model to address its limitations. For instance, McLean [13] proposed enhancements to account for the microstructural changes during deformation, such as grain boundary sliding and dislocation interactions. Multi-phase materials, such as composite systems [14], have also been modeled using adaptations of the Norton law.

In addition, Turner and McNeil [15] extended the model to incorporate stress relaxation behavior in polymers, thereby broadening its applicability beyond metals. Esposito and Bonera [16] demonstrated the use of the Norton model in finite element simulations to study creep in nickel-based superalloys, highlighting its potential for computational applications.

The increasing computational power in recent decades has facilitated the integration of the Norton model into numerical simulations. For instance, finite element modeling based on the Norton law has been widely used to simulate multi-axial creep behavior in structural components, such as pressure vessels and pipeline and keep [17].

Applied predominantly in the study of metals, ceramics, and geomechanics, the Norton model has been validated across various experimental and theoretical investigations [18]. The equation governing this model is of particular relevance in high-temperature applications where steady-state creep is prominent, as it reflects the flow characteristics of materials without a clear yield point [19]. The problem of a hollow sphere under hydrostatic loading within the Norton model is explored here to further understand the implications of such visco-plastic behavior on overall material stability and durability under high-temperature conditions [20].

1.3. The toothpaste flow problem

The study of semi-solid and semi-fluid materials, including colloids and pastes, introduces complexities beyond those encountered in purely fluid or solid mechanics. These materials exhibit yield behavior, transitioning from a solid-like state to a fluid state under sufficient stress [21]. This characteristic is exemplified in the toothpaste flow problem, where the material remains stationary until a yield stress threshold is reached, beyond which it flows as a visco-plastic medium [22]. Numerous studies have explored the flow characteristics of such materials in confined geometries, including tubes or channels, using models that incorporate yield stress and elasticity [23,24].

This behavior is particularly relevant for modeling the flow of highly viscous pastes and gels, where a balance between elasticity and flow defines the material response [25,26]. Recent work in visco-elastic flow has employed linear visco-elastic models such as the Norton model, as well as more complex approaches incorporating the Bingham or Herschel–Bulkley models for non-Newtonian fluids [27,28]. The stationary flow regime described in this study offers a precise solution to the toothpaste problem, providing comprehensive insights into stress distribution and flow behavior [29,30].

1.4. Gradient-enhanced visco-plastic models

As stated previously, visco-plasticity models have long been used to describe the time-dependent deformation of materials under stress, particularly in scenarios involving high-temperature creep and plasticity. While traditional visco-plastic models, such as the Norton model [11], are widely used, they often fail to capture the complex, size-dependent, and localized deformation behaviors observed in materials with intricate microstructures. These limitations arise because traditional models typically assume local stress and strain relationships that cannot fully account for non-local effects. This issue becomes particularly pronounced in materials with fine microstructures or under extreme loading conditions, where gradients in stress and strain play a crucial role.

In this paper, we introduce an enhanced framework by incorporating gradient terms into the visco-plasticity model, addressing these shortcomings. By introducing second-gradient terms, our model accounts for the non-local nature of the material response, smoothing discontinuities in the stress and deformation profiles. This gradient-enhanced formulation allows for the capture of size-dependent effects and provides a more accurate representation of materials under highly localized loading conditions, such as those observed in microscale deformations or materials with complex microstructures.

We apply this gradient-enhanced visco-plastic model to two benchmark problems: the deformation of a hollow sphere under hydrostatic loading and the flow dynamics of a material with toothpaste-like properties. The first scenario investigates the stress distribution and deformation behavior of a hollow sphere using the Norton viscous model, while the second scenario models the flow of a visco-plastic medium through a confined geometry, such as a cylindrical tube. For both problems, we derive analytical solutions that demonstrate the enhanced predictive capabilities of the gradient-enhanced model compared to traditional formulations.

This paper is structured as follows: First, we present the Norton model without gradient terms, providing a foundational understanding of its application to visco-plastic behavior. Next, we solve the toothpaste problem, first without and then with gradient terms, analyzing the impact of the gradient-enhanced framework on handling localization in the material’s flow. We then investigate the deformation of a hollow sphere, again considering both the standard Norton model and the gradient-enhanced model, and discuss how the inclusion of gradient terms improves the model’s ability to handle localized deformation. For each of these problems, we emphasize the role of the gradient framework in effectively addressing localization issues and capturing size-dependent behaviors.

The results of our study reveal that the inclusion of gradient terms significantly influences the material’s stress diffusion and plastic strain localization. The gradient coefficient, in particular, plays a critical role in controlling these effects, offering new insights into material behavior at micro- and meso-scales. These findings have important implications for engineering applications, particularly in the design and processing of materials where localized deformation and size-dependent behaviors are significant.

2. Elasto-visco-plastic response under uniaxial loading

The material’s viscous characteristics, particularly evident at elevated temperatures or when subjected to low deformation rates, can be effectively emphasized by employing diverse uniaxial experiments characterized by small deformations. Among these experiments, the simplest involves conducting a tensile test at a specified deformation rate $\overset{\cdot}{ϵ}$ , a procedure easily executed by adjusting the displacement speed of the grips on the tensile machine. It becomes apparent that as the imposed deformation rate $\overset{\cdot}{ϵ}$ increases, the resulting stress–strain curve $σ (ϵ)$ exhibits a corresponding increase, thereby illustrating the influence of $\overset{\cdot}{ϵ}$ on the material’s behavior.

At the extreme, when the rate of strain $\overset{\cdot}{ϵ}$ reaches exceedingly high values, the behavior exhibits quasi-elastic tendencies, leading to the disappearance of visco-plastic deformation due to insufficient time for its occurrence. On the converse end, in the limit as $\overset{\cdot}{ϵ} \to 0$ , two scenarios emerge. In the case of low temperature, the stress–strain curve $σ (ϵ)$ converges towards a non-zero limit, aligning with the characteristic elasto-plastic $σ (ϵ)$ curve. Consequently, no visco-plastic deformation is conceivable beneath this limit curve, denoting the existence of a threshold, which corresponds to the conventional plasticity threshold. Conversely, under high temperatures, no lower limit is discernible for the $σ (ϵ)$ curve as $\overset{\cdot}{ϵ} \to 0$ . As a result, the behavior lacks a threshold, with visco-plastic flow manifesting irrespective of the applied stress, even when it is exceptionally low.

Another enduring experiment revolves around the phenomenon of creep. In this particular trial, a tensile bar undergoes a consistent stress application—an easily attainable condition, achieved, for instance, by merely suspending a weight. In the context of elasto-plastic behavior, the relationship between ${\overset{\cdot}{ϵ}}^{p}$ and $\overset{\cdot}{σ}$ and the latter’s nullity affirm the existence of ${\overset{\cdot}{ϵ}}^{p}$ (similarly to ${\overset{\cdot}{ϵ}}^{e}$ ). Consequently, the total deformation promptly adopts the value $ϵ = ϵ^{e} + ϵ^{p} ”$ and remains constant thereafter. Conversely, in cases of viscous behavior, ϵ immediately assumes the value $ϵ^{e}$ but undergoes subsequent evolution. To elaborate further, the visco-plastic deformation rate $ϵ^{vp}$ experiences an initial decrease (primary creep), followed by an increase (tertiary creep). The ultimate phase corresponds to material damage culminating in final rupture, a facet generally omitted in visco-plastic behavior models, particularly those outlined in the subsequent section.

A modification of the uncomplicated creep test involves the two-stage creep experiment. Commencing with creep under stress $σ_{1}$ , one awaits the attainment of secondary creep (at a rate ${\overset{\cdot}{ϵ}}^{vp}$ ), stopping short of reaching tertiary creep. Subsequently, the stress undergoes an abrupt reduction to a value $σ_{2}$ . After a certain duration, a second phase of secondary creep emerges at a strain rate ${\overset{\cdot}{ϵ}}_{2}^{vp} < {\overset{\cdot}{ϵ}}_{1}^{vp}$ (an outcome expected due to the lower magnitude of $σ_{2}$ compared to $σ_{1}$ ). In the transitional stage, if $σ_{2}$ is adequately low, the occurrence of creep hesitation might manifest: rather than following a continuous decline from ${\overset{\cdot}{ϵ}}_{1}^{vp}$ to ${\overset{\cdot}{ϵ}}_{2}^{vp}$ (as might be intuitively anticipated), ${\overset{\cdot}{ϵ}}_{2}^{vp}$ abruptly adopts a negative value (even in the presence of $σ_{2}$ ), then subsequently ascends again until ${\overset{\cdot}{ϵ}}_{2}^{vp} > 0$ , passing through the zero point.

In the exploration of uniaxial elasticity, a third experiment delves into the realm of relaxation. This particular experimental scenario involves the instantaneous imposition of a specific deformation, a task easily accomplished. Subsequently, the deformation is held constant, resulting in an immediate stress level dictated by elasticity. Over time, there is a gradual decline in stress, marking the transition from elastic deformation to visco-plastic deformation. The outcome hinges on whether the material exhibits behavior without a threshold, typically observed at high temperatures, or if a threshold is present, typically occurring at low temperatures. In the former scenario, the stress asymptotically approaches zero, while in the latter, it decreases to the threshold value without actually reaching zero.

3. Closed-form solution for toothpaste flow with an elasto-visco-plastic model

Leveraging the methodologies employed for modeling the elastic-plastic behavior of materials, we integrate the premise of strain rate partitioning, as elucidated in the ensuing equation. This strategic incorporation facilitates a more comprehensive understanding of the material’s mechanical response, accounting for the intricate interplay between strain rates and their respective contributions to the overall deformation process. Drawing upon the methods utilized in modeling the elastic-plastic behavior of materials, we introduce the concept of strain rate partitioning through the subsequent equation:

D = D^{e} + D^{vp} .

(1)

This addition enhances the sophistication of our analysis, providing a nuanced framework to discern the nuanced contributions of different strain rates to the overarching deformational dynamics of the material.

In this context, the rate of elastic deformation, denoted as $D^{e}$ , is intricately governed by a specific hypo-elasticity law. The current predicament lies in formulating an apt expression for the visco-plastic strain, $D^{vp}$ , that not only conforms to the stipulated conditions but also evolves as a progressively increasing function in response to stress.

Navigating this challenge requires a nuanced understanding of the intricate balance between elastic and visco-plastic behaviors. Crafting an expression for $D^{vp}$ that adheres meticulously to the requirement of a progressively increasing function necessitates a delicate interplay of material properties, demanding a comprehensive synthesis of theoretical foundations and empirical insights.

The simple visco-plastic model we consider here aligns with Norton’s behavior without a yield factor. In the uniaxial scenario, this law presents us with the following relationship:

D^{vp} = {\overset{\cdot}{ε}}_{0} {(\frac{σ}{σ_{0}})}^{n} .

(2)

In this equation, ${\overset{\cdot}{ε}}_{0}$ , $σ_{0}$ , and n are treated as the constants. It is important to note that the model, in actuality, is only reliant on two constants, n and ${\overset{\cdot}{ε}}_{0} σ_{0}^{- n}$ . We will, however, explain the necessity of introducing a third constant in the subsequent sections for comprehensive understanding.

The exponent n is typically greater than 1, and its value varies significantly, decreasing from high values at lower temperatures to figures closely approximating 1 near the melting point of the material. In many instances, n could be substantially high, on the order of 5 to 10, indicating that materials rarely exhibit linear viscous behavior, especially when contrasted with fluid mechanics.

It is noteworthy to mention that in this model, if the stress σ is kept constant, $D^{vp}$ will also remain unchanged. Consequently, this model is capable of replicating only the second, or stationary, stage of creep, but it fails to accurately simulate primary creep. Furthermore, it does not adequately represent the ’hesitation’ in creep initiation or the third stage of creep.

To generalize this law to the more complex tri-axial scenario, we shall assume that $D^{vp}$ is colinear with the deviatoric component of the stress, similar to the alignment observed in plasticity theory when the von Mises yield criterion is applied.

The expression of $D^{vp}$ can be written as follows:

D^{vp} = \frac{3}{2} {\overset{\cdot}{ε}}_{0} {(\frac{σ_{eq}}{σ_{0}})}^{n} \frac{s}{σ_{eq}}, σ_{eq} = {(\frac{3}{2} s : s)}^{1 ∕ 2} .

(3)

This equation demonstrates that $D^{vp}$ is proportional to the deviatoric tensor s and can be utilized to verify the behavior under uniaxial tension load using equation (3).

One intriguing property of the Norton visco-plastic model is its reduction to the rigid perfectly plastic von Mises model with a yield limit of $σ_{0}$ when we fix ${\overset{\cdot}{ε}}_{0}$ and allow n to approach infinity. This adjustment requires treating ${\overset{\cdot}{ε}}_{0}$ and $σ_{0}$ separately and incorporating them into the single parameter ${\overset{\cdot}{ε}}_{0} σ_{0}^{- n}$ . Consequently, as n tends to infinity, $D^{vp}$ becomes zero for $σ_{eq} < σ_{0}$ and infinite for $σ_{eq} > σ_{0}$ , resulting in unlimited flow when the material is produced. As a result, Norton’s model serves as a “visco-plastic regularizer” that mitigates numerical issues such as spurious mesh dependency effects inherent in the von Mises model.

A more sophisticated variant of the Norton model without yield is the Norton model with yield, which is written as:

v^{vp} \equiv v = \frac{3}{2} {\overset{\cdot}{ε}}_{0} {⟨ \frac{σ_{eq}}{σ_{0}} - 1 ⟩}^{n} \frac{s}{σ_{eq}},

(4)

where $⟨ ⟩$ denote the “positive part”: $⟨ x ⟩ = x$ if $x \geq 0$ , 0 if $x < 0$ . Thus, $v^{vp}$ vanishes if $σ_{eq} < σ_{0}$ (existence of a yield) but non-zero whenever $σ_{eq} > σ_{0}$ (note that, although a “yield exists”, there is no “criterion” the material should obey, $σ_{eq}$ could go beyond the yield $σ_{0}$ ). Note that this model depends on three parameters ${\overset{\cdot}{ε}}_{0}, σ_{0}$ , and n. Despite the refinement introduced by the existence of a yield, this model does not reproduce the secondary creep (because $σ = cst$ implies that $D^{vp}$ is also constant). It also does not reproduce hesitation in creep. In addition, it reduces to the perfect rigid plastic von Mises model in the limit ${\overset{\cdot}{ε}}_{0} \to \infty$ , $σ_{0}$ and n fixed (since $v^{vp}$ ) vanishes for $σ_{eq} < σ_{0}$ and very large for $σ_{eq} > σ_{0}$ .

There exist other models, more complicated. Among these, we shall only mention here the Chaboche model (at least in one of its versions), particularly interesting since it is a relatively simple model which can reproduce (contrary to the Norton models) the primary creep phenomenon and the “hesitation” to creep phenomenon. The main equation of the Chaboche model reads:

\begin{matrix} {\begin{matrix} D^{vp} & = & \frac{3}{2} {\overset{\cdot}{ε}}_{0} {(\frac{σ_{eq}}{σ_{0}})}^{n} \frac{s - α}{σ_{eq}} \\ \tilde{α} & = & \frac{2}{3} h D^{vp} - c α, \end{matrix} \end{matrix}

(5)

where $σ_{eq} = {(\frac{3}{2} (s - α) : (s - α))}^{1 ∕ 2}$ is the generalized equivalent stress.

This model introduces an internal (deviatoric) tensorial parameter α, as in plasticity model accounting for kinematics hardening. It depends on four constants, ${\overset{\cdot}{ε}}_{0} σ_{0}^{- n}$ , n, h, and c. The term $- c α$ in the expression of $\tilde{α}$ (the Jaumann derivative of α) represents a viscous relaxation of the parameter α, of characteristic time $1 / c$ .

To confirm that this model accurately represents both the primary creep phase and the “hesitation” to creep phenomena, we need to formulate the equations for the uniaxial case, considering the assumption of small deformations.

Given the configuration of $D^{vp} \equiv {\overset{\cdot}{ε}}^{vp}$ , σ, and s in such a scenario (as outlined previously), and taking into account that $α = α e_{1} \otimes e_{1} - \frac{α}{2} e_{2} \otimes e_{2} - \frac{α}{2} e_{3} \otimes e_{3}$ and $σ_{eq} = | | (σ - \frac{3}{2} α | |$ , we can reformulate the equations as follows:

\begin{matrix} {\begin{matrix} {\overset{\cdot}{ε}}^{vp} & = & {\overset{\cdot}{ε}}_{0} {(\frac{| | (σ - \frac{3}{2} α | |}{σ_{0}})}^{n} \frac{σ - \frac{3}{2} α}{| | (σ - \frac{3}{2} α | |} = {\overset{\cdot}{ε}}_{0} {(\frac{| | (σ - \frac{3}{2} α | |}{σ_{0}})}^{n} sgn (σ - \frac{3}{2} α) \\ \overset{\cdot}{α} & = & \frac{2}{3} h {\overset{\cdot}{ε}}^{vp} - c α . \end{matrix} \end{matrix}

(6)

Initially, ( $α = 0$ ), and therefore, we have ${\overset{\cdot}{ε}}^{vp} = {\overset{\cdot}{ε}}_{0}^{vp} = {\overset{\cdot}{ε}}_{0} {(σ / σ_{0})}^{n}$ . But later on, α increases, and thus ${\overset{\cdot}{ε}}^{vp}$ decreases (primary creep). The secondary creep is reached when α a stationary value $α_{\infty}$ ( $\overset{\cdot}{α}$ ), ${\overset{\cdot}{ε}}^{vp} = {\overset{\cdot}{ε}}_{1 \infty}^{vp}$ becomes constant, less than ${\overset{\cdot}{ε}}_{0}^{vp}$ .

To analyze the “hesitation” to creep phenomena, let us consider that we have arrived at the secondary creep stage, corresponding to the stress $σ_{1}$ . In this case, α and ${\overset{\cdot}{ε}}^{vp}$ are the constant and equivalent to $α_{1 \infty}$ and ${\overset{\cdot}{ε}}_{1 \infty}^{vp}$ , respectively. Now, assume that a secondary stress, denoted by $σ_{2}$ and less than $\frac{3}{2} α_{1 \infty}$ , is applied. Under such circumstances, ${\overset{\cdot}{ε}}_{\infty}^{vp}$ turns negative, illustrating the hesitation or resistance to the creep phenomenon. However, as time passes by α goes toward a second stationary value $α_{2 \infty}$ , and ${\overset{\cdot}{ε}}^{vp}$ goes toward ${\overset{\cdot}{ε}}_{2 \infty}^{vp}$ , which is positive but less than ${\overset{\cdot}{ε}}_{1 \infty}^{vp}$ (second secondary creep).

To delve into more details, we aim to derive the precise equations representing primary and secondary creep stages, specifically in scenarios where $n = 1$ . Not only does this value align well with physical reality, but it also facilitates analytical computation. In such cases, the term $σ - \frac{3}{2} α$ consistently stays positive, thereby allowing us to express the flow rule as follows:

{\overset{\cdot}{ε}}^{vp} = \frac{{\overset{\cdot}{ε}}_{0}}{σ_{0}} (σ - \frac{3}{2} α) .

(7)

Using the expression of ${\overset{\cdot}{ε}}^{vp}$ in the evolution equation of α, we get:

\overset{\cdot}{α} = \frac{2}{3} h \frac{{\overset{\cdot}{ε}}_{0}}{σ_{0}} (σ - \frac{3}{2} α) - c α = \frac{2}{3} h {\overset{\cdot}{ε}}_{0} \frac{σ}{σ_{0}} - (c + h \frac{{\overset{\cdot}{ε}}_{0}}{σ_{0}}) α .

(8)

The stationary solution of this equation can be written as:

α = α_{\infty} \equiv \frac{2}{3} \frac{h {\overset{\cdot}{ε}}_{0} / σ_{0}}{h {\overset{\cdot}{ε}}_{0} / σ_{0} + c} σ .

(9)

Using the initial condition $α (0) = 0$ , the intermediate solution can be written as:

α = α_{\infty} (1 - e^{- (h {\overset{\cdot}{ε}}_{0} / σ_{0} + c) t}) = \frac{2}{3} \frac{h {\overset{\cdot}{ε}}_{0} / σ_{0}}{h {\overset{\cdot}{ε}}_{0} / σ_{0} + c} σ (1 - e^{- (h {\overset{\cdot}{ε}}_{0} / σ_{0} + c) t}) .

(10)

The value of ${\overset{\cdot}{ε}}^{vp}$ is obtained using the flow rule:

{\overset{\cdot}{ε}}_{0} \frac{σ}{σ_{0}} \frac{(h {\overset{\cdot}{ε}}_{0} / σ_{0}) e^{- (h {\overset{\cdot}{ε}}_{0} / σ_{0} + c) t} + c}{h {\overset{\cdot}{ε}}_{0} / σ_{0} + c} .

(11)

From the above discussion, we observe that ${\overset{\cdot}{ε}}^{vp}$ varies with time, transitioning from an initial value of ${\overset{\cdot}{ε}}_{0}^{vp} = {\overset{\cdot}{ε}}_{0} \frac{σ}{σ_{0}}$ to a final value of ${\overset{\cdot}{ε}}_{\infty}^{vp} = {\overset{\cdot}{ε}}_{0} \frac{σ}{σ_{0}} \frac{c}{h {\overset{\cdot}{ε}}_{0} / σ_{0} + c}$ , which is less than ${\overset{\cdot}{ε}}_{0}^{vp}$ .

This analysis highlights the role of the viscous relaxation of α in producing the secondary creep phenomena. In the absence of this factor (i.e., $c = 0$ ), the equation would result in ${\overset{\cdot}{ε}}_{\infty}^{vp} = 0$ , suggesting that the model could only account for the primary creep stage and not the secondary one.

4. Norton viscous model applied to cylindrical toothpaste flow

We are investigating the flow of toothpaste within a cylindrical tube, aligned along the Oz axis, with a radius R (see Figure 1). This material, as its name suggests, possesses a pasty consistency, which means it falls between a fluid and a solid. It is adequately described by a rigid, linear viscoelastic Norton law (with $n = 1$ ) that includes yield stress. The flow under consideration is stationary in nature, and the only non-zero component of the velocity is $v_{z} = v (x, z)$ . This component is zero on the inner wall of the tube, reflecting the behavior of a viscous fluid. We assume that $v > 0$ , which means that the fluid flows to the right.

Figure 1.

Geometry of the toothpaste model problem.

In the context of incompressible materials, where elasticity is neglected and the incompressible visco-elasticity rule is applied, the dynamics are governed by certain constraints. One crucial implication of these constraints is the incompressibility of the material, which has significant consequences on the behavior of the velocity variable v.

In the scenario described, the incompressibility condition imposes restrictions on how the velocity variable v can vary. Specifically, due to the incompressibility assumption, the velocity variable v becomes solely dependent on the radial coordinate r. This restriction is essential in understanding the material’s response to deformation, as it simplifies the representation of velocity in the system.

Digging deeper into the mechanics, when considering the strain tensor, a fundamental quantity characterizing deformation, the incompressibility assumption further narrows down the non-zero components. In this case, the only non-zero component of the strain tensor is $d_{rz} = \frac{1}{2} \frac{dv}{dr}$ . This expression signifies the rate of change of velocity with respect to the radial coordinate, emphasizing the axial-radial coupling inherent in the system.

By recognizing and analyzing these constraints imposed by incompressibility, researchers and engineers can gain valuable insights into the behavior of the material and design more accurate models for predicting its response under various conditions. The simplifications introduced by neglecting elasticity and applying incompressible viscoelasticity rules pave the way for a focused understanding of the mechanics at play, facilitating the development of effective solutions and optimizations in diverse engineering and scientific applications.

First, let us consider the visco-plastic zone (where $σ_{eq} > σ_{0}$ ). In this zone, $⟨ \frac{σ_{eq}}{σ_{0}} - 1 ⟩ \neq 0$ , and therefore, the flow rule implied that the only non-zero component of v is $s_{rz} \equiv σ_{rz} < 0$ (since it is clear that from the physics point of view $d_{zr} = \frac{1}{2} \frac{dv}{dr} < 0$ , see Figure 1). Thus, $σ_{eq} = {(\frac{3}{2} 2 σ_{rz}^{2})}^{1 ∕ 2} = - \sqrt{3} σ_{rz}$ so that the component “rz” of the flow rule can be written as:

d_{rz} = \frac{1}{2} \frac{dv}{dr} = \frac{\sqrt{3}}{2} {\overset{\cdot}{ε}}_{0} (\sqrt{3} \frac{σ_{rz}}{σ_{0}} + 1),

(12)

after substitution and simplification. From equation (12), we deduce that:

σ_{rz} = σ_{0} (\frac{1}{3 {\overset{\cdot}{ε}}_{0}} \frac{dv}{dr} - \frac{1}{\sqrt{3}}),

(13)

which shows that $σ_{rz}$ depends only on the coordinate r

Next, let us find the pressure $p = - \frac{1}{3} tr (σ)$ and the velocity fields in the visco-plastic zone. Due to the cylindrical properties involved in the problem, the pressure p depends only on r and z. The radial and axial balance equations (note that the orthoradial equation is automatically satisfied) are as follows:

\begin{matrix} {\begin{matrix} \frac{\partial σ_{rr}}{∂r} + \frac{\partial σ_{rz}}{∂z} + \frac{σ_{rr} - σ_{θθ}}{r} & = & 0 \\ \frac{\partial σ_{rz}}{∂r} + \frac{\partial σ_{rz}}{∂z} + \frac{σ_{rz}}{r} & = & 0, \end{matrix} \end{matrix}

(14)

which implies

\begin{matrix} {\begin{matrix} \frac{∂p}{∂r} & = & 0 \\ \frac{\partial σ_{rz}}{∂r} - \frac{d p}{d z} + \frac{σ_{rz} (r)}{r} & = & 0, \end{matrix} \end{matrix}

(15)

given that $σ = - p 1 + s$ . Thus, the pressure depends only on the coordinate z. In addition, the second equation of equation (25) is written as:

(\frac{\partial σ_{rz}}{∂r} + \frac{σ_{rz} (r)}{r}) (r) = \frac{d p}{d z} (z),

(16)

which implies that the left and right hand sides of these equations equal the same constant, say $- P$ which represent the drop of the pressure (due to the viscosity) by unit length of the tube. Clearly, $P > 0$ for a flow towards the right. Given the expression of $σ_{rz}$ as a function of v, we get:

σ_{0} (\frac{1}{3 {\overset{\cdot}{ε}}_{0}} \frac{d^{2} v}{d r^{2}} + \frac{1}{3 {\overset{\cdot}{ε}}_{0}} \frac{1}{3} \frac{d v}{d r} - \frac{1}{\sqrt{3} r}) = - P,

(17)

which implies that:

\frac{d^{2} v}{d r^{2}} + \frac{1}{r} \frac{d v}{d r} = - 3 {\overset{\cdot}{ε}}_{0} \frac{P}{σ_{0}} + \sqrt{3} \frac{{\overset{\cdot}{ε}}_{0}}{r} .

(18)

Let us find a quadratic solution in r (by analogy with Couette flow of a viscous fluid) which vanishes on $r = R$ , i.e., $v = - \frac{a}{2} {(r - R)}^{2} - b (r - R)$ . We then get $\frac{d v}{d r} = - a (r - R) - b,$ $\frac{d^{2} v}{d r^{2}} = - a$ and the previous equation becomes $- a - a (1 - \frac{R}{r}) - \frac{b}{r} = - 2 a + \frac{aR - b}{r} = - 3 {\overset{\cdot}{ε}}_{0} \frac{P}{σ_{0}} + \frac{\sqrt{3} {\overset{\cdot}{ε}}_{0}}{r}$ which yields:

\begin{matrix} {\begin{matrix} a & = & \frac{3}{2} {\overset{\cdot}{ε}}_{0} \frac{P}{σ_{0}} \\ b & = & aR - \sqrt{3} {\overset{\cdot}{ε}}_{0} = {\overset{\cdot}{ε}}_{0} (\frac{3}{2} \frac{PR}{σ_{0}} - \sqrt{3}) . \end{matrix} \end{matrix}

(19)

Therefore, in the visco-plastic zone, $v (r)$ becomes:

v (r) = - \frac{3}{4} {\overset{\cdot}{ε}}_{0} \frac{P}{σ_{0}} {(r - R)}^{2} - {\overset{\cdot}{ε}}_{0} (\frac{3}{2} \frac{PR}{σ_{0}} - \sqrt{3}) (r - R) .

(20)

This velocity profile is illustrated in Figure 2 for various values of the strain rates.

Figure 2.

Velocity field as a function of the radius for different strain regimes of the material.

We can now find the internal radius of the visco-plastic zone. This radius corresponds to a point where $σ_{eq} = σ_{0}$ , which implied that $d_{rz} = 0$ and $\frac{dv}{dr} (R_{0}) = 0$ which as it turn gives the equality $- a (R_{0} - R) - b = 0$ . With this we have:

R_{0} = R - \frac{b}{a} = R + \frac{{\overset{\cdot}{ε}}_{0} (\frac{3}{2} \frac{PR}{σ_{0}} - \sqrt{3})}{\frac{3}{2} {\overset{\cdot}{ε}}_{0} \frac{P}{σ_{0}}} = \frac{2}{\sqrt{3}} \frac{σ_{0}}{P} .

(21)

Note that $R_{0} = \frac{2}{\sqrt{3}} \frac{σ_{0}}{P}$ is indeed the internal radius since $\frac{d v}{d r}$ is a decreasing function of r, which vanishes for $r = R_{0}$ . We then have for $r > R_{0}$ , $\frac{d v}{d r} < 0$ , this implies that:

σ_{rz} = σ_{0} (\frac{1}{3 {\overset{\cdot}{ε}}_{0}} \frac{d v}{d r} - \frac{1}{\sqrt{3}}) < \frac{σ_{0}}{\sqrt{3}} .

(22)

With this $σ_{eq} = - \sqrt{3} σ_{rz} > σ_{0}$ . The region $r < R_{0}$ is in the contrary rigid: it is where visco-plastic we would have $\frac{d v}{d r} > 0$ which should $σ_{eq}$ by the same reasoning.

In order for the material to be visco-plastic, the internal radius $R_{0}$ of the visco-plastic zone must be smaller than the radius of the tube R, i.e., $\frac{2}{\sqrt{3}} \frac{σ_{0}}{P} < R$ which implies $P > \frac{2}{\sqrt{3}} \frac{σ_{0}}{R}$ . There is no material flow whenever $P < \frac{2}{\sqrt{3}} \frac{σ_{0}}{R}$ : all the material remains, see Figure 3. In the case where there is a material flow ( $P > \frac{2}{\sqrt{3}} \frac{σ_{0}}{R})$ , the velocity profile looks like. Note that we recovered the parabolic velocity profile of the well-known Couette flow for a viscous fluid by letting the yield limit $σ_{0}$ tending toward 0, the “viscosity coefficient” $σ_{0} ∕ {\overset{\cdot}{ε}}_{0}$ being kept constant.

σ_{rz} = σ_{0} (\frac{1}{3 {\overset{\cdot}{ε}}_{0}} \frac{d v}{d r} - \frac{1}{\sqrt{3}}) = σ_{0} (\frac{1}{3 {\overset{\cdot}{ε}}_{0}} (- a (r - R) - b) - \frac{1}{\sqrt{3}}) = - \frac{\Pr}{2} .

(23)

Figure 3.

Velocity profile of the toothpaste model, showing various regions within the tube.

In the rigid region, see Figure 3, we introduce the assumption that only the non-zero component of s is $σ_{sz}$ only depends on r, as in the visco-plastic region. The balance equations read:

\begin{matrix} {\begin{matrix} \frac{∂p}{∂r} & = & 0 \\ \frac{\partial σ_{rz}}{∂r} - \frac{d p}{d z} + \frac{σ_{rz} (r)}{r} & = & 0, \end{matrix} \end{matrix}

(24)

which implies that:

\begin{matrix} {\begin{matrix} \frac{∂p}{∂z} & = & P^{'} \\ \frac{\partial σ_{rz}}{∂r} + \frac{σ_{rz} (r)}{r} & = & P^{'}, \end{matrix} \end{matrix}

(25)

where $P^{'}$ is a new constant (“a priori” different P). Integrating the second equation, we get:

\frac{\partial σ_{rz}}{∂r} + σ_{rz} (r) = \frac{d}{d r} (r σ_{rz}) = - P r^{'},

(26)

which implies that:

r σ_{rz} = - P^{'} \frac{r^{2}}{2} + A

(27)

and

σ_{rz} = - P^{'} \frac{r}{2} + \frac{A}{r},

(28)

As $σ_{rz}$ cannot vanish at infinity at 0, we necessary have $A = 0$ , and thus, $σ_{rz} = - P^{'} \frac{r}{2}$ . The continuity of $σ_{rz}$ in $r = R_{0}$ then implies $P = P^{'}$ . The formula:

σ_{rz} = - P \frac{r}{2},

(29)

is then valid everywhere in the tube.

Note that $σ_{rz}$ is independent of the yield $σ_{0}$ : the distribution of the stresses is exactly the same as for a Newtonian viscous fluid.

5. Gradient-enhanced toothpaste flow analysis

The flow of semi-solid materials, such as toothpaste, is traditionally modeled using visco-plastic models that capture the transition from solid-like to fluid-like behavior. However, these models may fail to accurately represent localized phenomena near rigid-visco-plastic boundaries. This work extends the model by incorporating second-gradient effects, providing additional terms that account for higher-order stress and strain variations. The benefits of this extension include improved numerical stability and a more realistic representation of material behavior.

This section presents an extended derivation and analysis of visco-plastic flow in a cylindrical tube, incorporating second-gradient effects to model localized stress behavior near the rigid-visco-plastic boundary. The primary focus is on refining the velocity profile and stress distributions by considering higher-order stress terms and validating the solutions numerically. The application discussed uses toothpaste flow in a tube as a case study.

5.1. Governing equations

5.1.1. Stress decomposition

The stress tensor σ is expressed as a sum of Cauchy stress and second-gradient contributions:

σ = σ^{Cauchy} + σ^{SG},

where:

$σ^{Cauchy}$ is the classical Cauchy stress tensor.

$σ^{SG} = η \frac{\partial^{2} d}{\partial r^{2}}$ , with η being a material parameter representing second-gradient effects.

In cylindrical coordinates, for axial symmetry, we focus on the shear stress $σ_{rz}$ :

σ_{rz} = σ_{rz}^{Cauchy} + η \frac{\partial^{2} d_{rz}}{\partial r^{2}} .

Appendix 4 provides a variational-oriented justification of this decomposition.

5.1.2. Linear momentum balance

The balance of linear momentum in the axial-radial plane is given by:

\frac{\partial σ_{rz}}{∂r} + \frac{σ_{rz}}{r} = \frac{dp}{dz} .

Substituting the stress decomposition:

\frac{\partial}{∂r} (σ_{rz}^{Cauchy} + η \frac{\partial^{2} d_{rz}}{\partial r^{2}}) + \frac{σ_{rz}^{Cauchy} + η \frac{\partial^{2} d_{rz}}{\partial r^{2}}}{r} = \frac{dp}{dz} .

Expanding the derivatives:

\frac{\partial σ_{rz}^{Cauchy}}{∂r} + \frac{σ_{rz}^{Cauchy}}{r} + η \frac{\partial^{3} d_{rz}}{\partial r^{3}} + \frac{η}{r} \frac{\partial^{2} d_{rz}}{\partial r^{2}} = \frac{dp}{dz} .

5.1.3. Strain rate relationship

The shear strain rate $d_{rz}$ is related to the velocity gradient by:

d_{rz} = \frac{1}{2} \frac{dv}{dr} .

Thus, its higher derivatives are given by:

\frac{\partial^{2} d_{rz}}{\partial r^{2}} = \frac{1}{2} \frac{d^{3} v}{d r^{3}}, \frac{\partial^{3} d_{rz}}{\partial r^{3}} = \frac{1}{2} \frac{d^{4} v}{d r^{4}} .

Substituting into the momentum balance equation:

\frac{\partial σ_{rz}^{Cauchy}}{∂r} + \frac{σ_{rz}^{Cauchy}}{r} + \frac{η}{2} \frac{d^{4} v}{d r^{4}} + \frac{η}{2 r} \frac{d^{3} v}{d r^{3}} = \frac{dp}{dz} .

5.1.4. Stress in the visco-plastic zone

In the visco-plastic zone ( $r > R_{0}$ ), the Cauchy stress $σ_{rz}^{Cauchy}$ is governed by the visco-plastic flow rule:

\frac{dv}{dr} = \frac{\sqrt{3}}{2} ϵ_{0} (\frac{\sqrt{3} σ_{rz}^{Cauchy}}{σ_{0}} + 1) .

Rearranging:

σ_{rz}^{Cauchy} = σ_{0} (\frac{1}{3 ϵ_{0}} \frac{dv}{dr} - \frac{1}{\sqrt{3}}) .

Substitute this expression into the momentum balance:

\frac{\partial}{∂r} (σ_{0} (\frac{1}{3 ϵ_{0}} \frac{dv}{dr} - \frac{1}{\sqrt{3}})) + \frac{σ_{0}}{r} (\frac{1}{3 ϵ_{0}} \frac{dv}{dr} - \frac{1}{\sqrt{3}}) + \frac{η}{2} \frac{d^{4} v}{d r^{4}} + \frac{η}{2 r} \frac{d^{3} v}{d r^{3}} = \frac{dp}{dz} .

Simplifying:

\frac{σ_{0}}{3 ϵ_{0}} \frac{d^{2} v}{d r^{2}} + \frac{σ_{0}}{r} (\frac{1}{3 ϵ_{0}} \frac{dv}{dr} - \frac{1}{\sqrt{3}}) + \frac{η}{2} \frac{d^{4} v}{d r^{4}} + \frac{η}{2 r} \frac{d^{3} v}{d r^{3}} = \frac{dp}{dz} .

5.1.5. Boundary layer correction near $r = R_{0}$

Near the rigid-visco-plastic boundary, the second-gradient contributions dominate. Assuming a small perturbation around $r = R_{0}$ :

\frac{d^{2} v}{d r^{2}} + α \frac{d^{4} v}{d r^{4}} = \frac{3 ϵ_{0}}{σ_{0}} P, α = \frac{3 ϵ_{0} η}{2 σ_{0}} .

5.2. Solution and trial function

To solve this equation, we use a trial solution:

v (r) = c_{1} {(r - R_{0})}^{2} + c_{2} {(r - R_{0})}^{4} .

Taking the derivative of the trial function and substituting these derivatives into the governing equation, we get:

2 c_{1} + 12 c_{2} {(r - R_{0})}^{2} + α (24 c_{2}) = \frac{3 ϵ_{0}}{σ_{0}} P .

At $r = R_{0}$ (where $(r - R_{0}) = 0$ ):

2 c_{1} + 24 α c_{2} = \frac{3 ϵ_{0}}{σ_{0}} P .

For a small perturbation $(r - R_{0} = δ)$ :

2 c_{1} + 12 c_{2} δ^{2} + 24 α c_{2} = \frac{3 ϵ_{0}}{σ_{0}} P .

From the first equation:

c_{1} = \frac{\frac{3 ϵ_{0}}{σ_{0}} P - 24 α c_{2}}{2} .

From the second equation:

c_{2} = \frac{\frac{3 ϵ_{0}}{σ_{0}} P - 2 c_{1}}{12 δ^{2} + 24 α} .

5.3. Numerical implementation

To validate the proposed model, we numerically solve the governing equation:

\frac{d^{2} v}{d r^{2}} + α \frac{d^{4} v}{d r^{4}} = \frac{3 ϵ_{0}}{σ_{0}} P,

with the boundary conditions:

v (R_{0}) = 0, \frac{dv}{dr} (R_{0}) = 0 .

The solution is implemented using finite differences to approximate the derivatives, discretizing the spatial domain $[R_{0}, R_{\max}]$ into N equally spaced points with step size $Δ r$ .

5.3.1. Finite difference discretization

The second derivative is approximated using a central difference scheme:

\frac{d^{2} v}{d r^{2}} \approx \frac{v_{i + 1} - 2 v_{i} + v_{i - 1}}{Δ r^{2}},

and the fourth derivative is discretized as:

\frac{d^{4} v}{d r^{4}} \approx \frac{v_{i + 2} - 4 v_{i + 1} + 6 v_{i} - 4 v_{i - 1} + v_{i - 2}}{Δ r^{4}} .

The governing equation at each interior point i becomes:

\frac{v_{i + 1} - 2 v_{i} + v_{i - 1}}{Δ r^{2}} + α \frac{v_{i + 2} - 4 v_{i + 1} + 6 v_{i} - 4 v_{i - 1} + v_{i - 2}}{Δ r^{4}} = \frac{3 ϵ_{0}}{σ_{0}} P .

Boundary conditions are imposed explicitly:

At $R_{0}$ : The Dirichlet boundary condition $v (R_{0}) = 0$ is enforced by setting $v_{1} = 0$ .

At $R_{0}$ : The Neumann boundary condition $\frac{dv}{dr} (R_{0}) = 0$ is discretized as:

\frac{v_{2} - v_{1}}{Δ r} = 0 .

This leads to $v_{2} = v_{1}$ , which is incorporated into the system matrix.

At $R_{\max}$ : The Dirichlet condition $v (R_{\max}) = 0$ is enforced by setting $v_{N} = 0$ .

5.3.2. Matrix formulation

The system is represented as a linear system:

A v = b,

where:

A is the finite difference coefficient matrix derived from the discretized equations.

v is the vector of unknowns ${[v_{1}, v_{2}, \dots, v_{N}]}^{⊤}$ .

b is the right-hand side vector containing the source term $\frac{3 ϵ_{0}}{σ_{0}} P$ .

To improve numerical stability, a small regularization term ( $ϵ = 1 0^{- 8}$ ) is added to the diagonal of A, ensuring that the matrix remains well-conditioned for inversion.

5.3.3. Numerical solution procedure

The steps for solving the system are as follows:

Define the spatial grid by discretizing the radial coordinate r with a step size $Δ r = \frac{R_{\max} - R_{0}}{N - 1}$ , where $R_{0}$ and $R_{\max}$ represent the minimum and maximum values of the radial domain, respectively, and N is the total number of grid points.

Construct the system matrix A and the right-hand side vector b by applying finite difference approximations to the governing differential equations, ensuring that the discretization reflects the specific problem being solved.

Enforce the boundary conditions by explicitly modifying the first and last rows of the matrix A and updating the corresponding entries in b. This step ensures that the solution satisfies the required conditions at the boundaries of the domain.

Solve the resulting linear system $A x = b$ using MATLAB’s backslash operator (∖), which efficiently computes the solution vector x.

5.4. Results and discussion

The computed solution $v (r)$ is plotted over the domain $[R_{0}, R_{\max}]$ , providing insights into the deformation profile. The accuracy of the numerical solution is validated by ensuring that it satisfies the boundary conditions and the original differential equation within the discretization error (see Figures 4 and 5).

Figure 4.

Velocity profile $V (r)$ , corrected by incorporating second-gradient effects near the rigid-visco-plastic boundary $r = R_{0}$ , with different values for α and a fixed value of P.

Figure 5.

Velocity profile $V (r)$ , corrected by incorporating second-gradient effects near the rigid-visco-plastic boundary $r = R_{0}$ for different values of the pressure P and $α = 2$ .

Next, we present a detailed analysis of the numerical results obtained for the velocity profile $v (r)$ , corrected by incorporating second-gradient effects near the rigid-visco-plastic boundary $r = R_{0}$ . The goal is to understand the impact of these effects on the material behavior and verify the model’s consistency with physical expectations.

5.4.1. Velocity profile near $r = R_{0}$

The velocity profile $v (r)$ starts at zero at $r = R_{0}$ , consistent with the no-slip boundary condition. In addition, the initial gradient $\frac{dv}{dr}$ is also zero at $r = R_{0}$ , reflecting the transition from rigid to visco-plastic behavior in the material. The inclusion of second-gradient terms introduces curvature corrections near $R_{0}$ , resulting in a more gradual increase in velocity beyond $R_{0}$ compared to classical models. This effect highlights the role of higher-order stress terms in smoothing the velocity gradient.

The second-gradient parameter $α = \frac{3 ϵ_{0} η}{2 σ_{0}}$ significantly affects the velocity profile:

Higher α: Larger values of α result in a smoother profile with reduced velocity gradients near $R_{0}$ . This behavior indicates increased resistance to shear deformation, emphasizing the importance of size-dependent effects.

Lower α: Smaller values of α lead to a sharper profile, approximating the classical quadratic form. This corresponds to weaker size effects and a more localized deformation response.

5.5. Comparison to classical visco-plastic models

For $α = 0$ (i.e., no second-gradient effects), the velocity profile reduces to the quadratic form given by:

v (r) = - \frac{3}{4} ϵ_{0} \frac{P}{σ_{0}} {(r - R)}^{2} - ϵ_{0} (\frac{3}{2} \frac{PR}{σ_{0}} - \sqrt{3}) (r - R) .

(30)

The addition of second-gradient terms refines this profile by introducing corrections proportional to ${(r - R)}^{4}$ , enhancing the model’s ability to capture localized microstructural behavior.

The inclusion of second-gradient effects has several key implications for material behavior:

Shear localization: Second-gradient terms distribute shear more evenly across the visco-plastic zone, preventing sharp velocity gradients near the boundary $R_{0}$ .

Material behavior: The corrected model aligns with the observed behavior of semi-solid materials, such as toothpaste, where microstructural resistance smoothens localized stress and deformation near boundaries.

The numerical solution demonstrates the following characteristics:

Stability: The inclusion of higher-order derivatives improves numerical stability, avoiding oscillations or mesh dependency in the velocity profile (see Appendix 1).

Accuracy: The numerical results closely match theoretical predictions, validating the correctness of the finite difference scheme used in the computations (here also see Appendix 1).

5.6. Velocity profile across the entire domain

The velocity profile $v (r)$ was solved numerically for the entire domain $r \in [0, R]$ , incorporating both the rigid and visco-plastic zones. The numerical solution reveals a distinct behavior in each zone, reflecting the underlying physics of the toothpaste flow model with second-gradient effects.

5.6.1. Rigid zone ( $r R_{0}$ )

In the rigid zone, where $r < R_{0}$ , the material exhibits no flow, consistent with the assumption of rigidity. The velocity remains zero throughout this region:

v (r) = 0, for r \in [0, R_{0}] .

(31)

This behavior confirms that the material is unable to deform or flow below the yield stress threshold.

5.6.2. Transition at $r = R_{0}$

At the interface between the rigid and visco-plastic zones, located at $r = R_{0}$ , the velocity profile transitions smoothly from zero to non-zero values. The boundary conditions enforce:

v (R_{0}) = 0, \frac{dv}{dr} |_{r = R_{0}} = 0 .

(32)

The inclusion of second-gradient terms results in a smooth, gradual transition across the boundary, reducing sharp gradients and providing a more realistic representation of the material’s mechanical response near the interface.

5.6.3. Visco-plastic zone ( $r \geq R_{0}$ )

In the visco-plastic zone, the governing equation for the velocity profile is given by:

\frac{d^{2} v}{d r^{2}} + α \frac{d^{4} v}{d r^{4}} = \frac{3 ϵ_{0}}{σ_{0}} P, α = \frac{3 ϵ_{0} η}{2 σ_{0}},

(33)

with appropriate boundary conditions at $r = R_{0}$ and $r = R$ . The numerical solution shows a parabolic increase in velocity, modified by higher-order terms due to the second-gradient effects:

Smooth profile: The velocity increases gradually from zero at $r = R_{0}$ , reflecting the material’s visco-plastic flow behavior (see Figures 4 and 5).

Impact of second-gradient terms: The second-gradient contributions lead to curvature corrections near the boundaries, smoothing the transition between the rigid and visco-plastic zones (here also see Figures 4 and 5).

5.6.4. Influence of the second-gradient parameter α

The parameter α controls the influence of second-gradient effects. Larger values of α produce a smoother profile with reduced velocity gradients near $r = R_{0}$ . Conversely, smaller values of α result in a sharper transition, approximating classical visco-plastic behavior without higher-order corrections. This demonstrates the role of α in capturing size-dependent and localized stress effects within the material.

5.7. Comparison with classical visco-plastic models

When second-gradient effects are neglected ( $α = 0$ ), the velocity profile reduces to the classical quadratic form observed in standard visco-plastic models:

v (r) = - \frac{3}{4} ϵ_{0} \frac{P}{σ_{0}} {(r - R)}^{2} - ϵ_{0} (\frac{3}{2} \frac{PR}{σ_{0}} - \sqrt{3}) (r - R) .

(34)

The inclusion of second-gradient terms provides a more accurate representation of the material behavior by incorporating microstructural resistance and distributing shear more evenly across the visco-plastic zone.

5.8. Physical implications

The numerical results highlight several important physical characteristics:

Shear localization: The second-gradient terms distribute shear stresses more evenly, reducing sharp gradients and enhancing numerical stability.

Realistic material behavior: The model aligns well with the observed behavior of semi-solid materials, such as toothpaste, where localized stress and deformation near boundaries are smoothed due to internal microstructural resistance.

5.9. Numerical stability and accuracy

The numerical solution demonstrates stability and accuracy in capturing both the primary velocity profile and the higher-order corrections, see Appendix 1. The use of fourth-order finite differences avoids numerical artifacts, such as oscillations, and ensures consistent results across different mesh resolutions.

6. Additional numerical results and discussion

This section presents some additional numerical results of the gradient-based model, obtained using the finite difference method to discretize the governing equations. The discussion includes the simulation parameters, the condition number of the matrix system, the velocity profiles, and the residual trends, with relevance to the context of gradient models. Figures 6 –9 summarize the results.

Figure 6.

Raw velocity profiles $v (r)$ for varying P. The profiles scale linearly with P, highlighting proportional system response.

Figure 7.

Normalized velocity profiles $ṽ (r)$ for varying P. Profiles collapse into a single curve, confirming consistent scaling and trends.

Figure 8.

Residuals of the governing equation for varying P. Residuals increase for higher P, highlighting the need for additional refinement at larger loads.

Figure 9.

Normalized residuals $Residual ∕ P$ for varying P. Consistent scaling across all P values indicates robust numerical stability.

6.1. Simulation setup and matrix conditioning

The governing equation is derived from a gradient-enhanced model designed to include higher-order terms that account for material gradient effects. The equation is written as:

α \frac{d^{4} v}{d r^{4}} + β \frac{d^{2} v}{d r^{2}} + γv = P,

(35)

where:

$v (r)$ : velocity profile;

$α = 1 \times 1 0^{- 7}$ : gradient stiffness coefficient;

$β = 5$ : damping coefficient;

$γ = 0$ : higher-order term neglected for simplicity;

P: load term.

The problem is defined over a non-dimensionalized domain $r \in [0, 1]$ , where:

\begin{matrix} v (1) & = 0 (Dirichlet boundary condition at r = 1), \end{matrix}

(36)

\begin{matrix} \frac{d v}{d r} |_{r = 0} & = 0 (Neumann boundary condition at r = 0) . \end{matrix}

(37)

The finite difference discretization uses a uniform grid with $N = 800$ points, resulting in a spatial step size of $Δ r = 0.00125$ . A small regularization term ( $1 \times 1 0^{- 5}$ ) is added to the diagonal of the matrix to enhance numerical stability. The linear system is solved for multiple P values ( $P = 1, 10, 50, 100, 500$ ).

The matrix condition number is a critical measure of numerical stability and accuracy. Across all simulations, the condition number remains consistent and well-controlled, with values on the order of $1 0^{6}$ to $1 0^{7}$ . For $P = 500$ , the condition number is approximately $9.05 \times 1 0^{6}$ , ensuring that the matrix system is not ill-conditioned. This stability is crucial for accurately resolving higher-order gradient terms.

6.2. Velocity profiles and residuals

The velocity profiles and residuals are analyzed in Figures 6 –9.

6.2.1. Raw velocity profiles

Figure 6 illustrates the raw velocity profiles $v (r)$ for different P. The profiles increase proportionally with the load P, reflecting the linear dependence of the solution on the forcing term. For $P = 500$ , the maximum velocity reaches approximately $250$ , while for $P = 1$ , it remains near zero. The smooth decay to $v (1) = 0$ aligns with the Dirichlet boundary condition. The matrix condition number’s consistency across simulations ensures that the velocity profiles remain stable and free from numerical artifacts.

6.2.2. Normalized velocity profiles

The normalized velocity profiles are presented in Figure 7. These profiles collapse into a single curve, confirming consistent scaling and numerical behavior across all P. The smooth decay from 1 at $r = 0$ to 0 at $r = 1$ validates the enforcement of boundary conditions. The normalization also highlights the relative trends in $v (r)$ , independent of the load magnitude.

6.2.3. Residuals and normalized residuals

Figure 8 shows the residuals for varying P. For $P = 1$ , the residuals remain small and positive, demonstrating good numerical balance. As P increases, the residuals grow, peaking at $\sim 600$ for $P = 500$ . This trend suggests that higher loads pose greater challenges for balancing the governing equation, despite the well-conditioned matrix.

Normalized residuals are plotted in Figure 9, where residuals are scaled by P. These results reveal consistent behavior across all loads, with normalized residuals peaking at approximately $1.5$ near $r = 0$ . The consistent scaling of residuals further demonstrates the stability of the numerical scheme.

6.3. Discussion on gradient models

Gradient-enhanced models include higher-order terms $(α \frac{d^{4} v}{d r^{4}})$ that capture spatial localization effects. For lower P, the gradient stiffness parameter α dominates, resulting in smooth velocity profiles. However, at higher P, the damping term $(β \frac{d^{2} v}{d r^{2}})$ and forcing term (P) overshadow the gradient effects, leading to higher residuals.

The matrix condition number remains stable even as the gradient term contributes to the stiffness. This stability highlights the robustness of the finite difference scheme. Future work could explore adaptive meshing or higher-order schemes to improve residual control for larger loads.

6.3.1. Contour plot of velocity profiles

The contour plot of the velocity profiles is presented in Figure 10. This plot visualizes the spatial variation of the velocity $v (r)$ along the non-dimensional radial coordinate r (horizontal axis) for different pressure values P (vertical axis). The smooth and continuous variation of $v (r)$ across the domain reflects the stability and effectiveness of the numerical solution. Also:

Scaling with pressure: The velocity profiles scale proportionally with P, with the maximum velocity near $r = 0$ increasing from near-zero values for $P = 1$ to approximately $250$ for $P = 500$ . This is consistent with the raw velocity profiles shown in Figure 6.

Boundary condition enforcement: The velocity decreases monotonically along r and reaches $v (1) = 0$ at the domain boundary, satisfying the Dirichlet boundary condition. Near $r = 0$ , the Neumann condition ensures a maximum velocity, where $\frac{d v}{d r} \to 0$ .

Localized effects: The smooth transition from high velocity near $r = 0$ to zero velocity at $r = 1$ highlights the importance of the gradient model in regularizing the solution. Without the inclusion of the higher-order gradient term ( $α \frac{d^{4} v}{d r^{4}}$ ), sharp transitions or oscillations might occur, particularly for high P.

Figure 10.

Contour plot of velocity profiles $v (r)$ for different pressure values P. The plot illustrates the spatial variation of velocity across the non-dimensional domain, demonstrating smooth scaling with increasing P. The gradient term ensures stability and regularization, particularly near the boundaries.

The contour plot also provides a visual link to the residuals discussed in Figure 8. The smoothness of $v (r)$ correlates with the small residuals for lower pressures ( $P = 1, 10$ ). However, as P increases, the forcing term dominates over the regularizing effect of the gradient term, leading to larger residuals at higher pressures ( $P = 100, 500$ ).

From a numerical perspective, the contour plot demonstrates the robustness of the finite difference scheme and the stability of the matrix system. The condition number remains consistently below $1 0^{7}$ for all simulations, even as the gradient term introduces additional stiffness. This ensures that the velocity profiles remain free of numerical artifacts and converge smoothly across the domain.

In conclusion, the results validate the gradient-enhanced model in capturing smooth and proportional velocity profiles across varying load conditions. The matrix condition number ( $\sim 1 0^{6}$ ) ensures stability and accuracy, while normalized residuals reveal consistent scaling. These findings support the model’s applicability for problems involving higher-order spatial effects.

6.4. Analysis of velocity profiles and gradient effects

In this section, we analyze the effect of the gradient term parameter α on the velocity profile $v (r)$ , the energy-like term $\int v^{2} (r) d r$ , and the numerical stability of the solution. The governing equation incorporates a fourth-order derivative term scaled by α, which introduces regularization effects, alongside a second-order damping term scaled by β.

6.4.1. Velocity profiles for varying α

Figure 11 shows the velocity profiles $v (r)$ for fixed pressure $P = 100$ and varying values of α. For small values of α (e.g., $1 0^{- 12}$ and $1 0^{- 10}$ ), the velocity profile remains nearly flat across the interior of the domain, forming a plateau, with a steep drop near the boundary at $r = R_{\max}$ . This behavior reflects the dominance of the second-order term scaled by β, which governs the primary damping mechanism in the system.

Figure 11.

Stabilized velocity profiles for fixed $P = 100$ and varying α. The profiles remain stable and smooth across all values of α.

The plateau in the interior arises because the second-order term counterbalances the forcing term P, resulting in a nearly uniform velocity distribution. The steep gradient at the boundary is due to the Dirichlet condition $v (R_{\max}) = 0$ , which enforces dissipation at the domain edge.

As α increases (e.g., $1 0^{- 6}$ ), the profile becomes smoother, reflecting the regularizing effect of the fourth-order term. However, for the chosen range of α, the velocity profiles exhibit minimal variation, suggesting that the second-order term and forcing term P dominate the solution.

6.4.2. Energy-like term vs. α

To quantify the effect of α on the overall solution, we compute the energy-like term $E (α) = \int v^{2} (r) d r$ . Figure 12 shows the variation of $E (α)$ with α on a log-log scale.

Figure 12.

Energy-like term $E (α) = \int v^{2} (r) d r$ vs. α. The decrease in $E (α)$ for larger α reflects the regularizing effect of the fourth-order term.

For small values of α, the energy-like term remains nearly constant, indicating that the fourth-order term has a negligible impact on the solution. As α increases, $E (α)$ decreases significantly, demonstrating the regularizing effect of the gradient term. This trend is consistent with the expected behavior of gradient models, where higher-order terms smooth the solution and reduce sharp variations.

6.4.3. Numerical stability

The numerical stability of the solution is validated through the condition number of the system matrix A and the residual norm. For small α, the condition number remains moderate ( $\sim 1 0^{11}$ ), and the residual norm is small ( $< 1 0^{- 3}$ ), indicating a stable solution. As α increases, the condition number grows significantly (e.g., $1 0^{13}$ for $α = 1 0^{- 4}$ ), highlighting the need for careful regularization.

Table 1 summarizes the condition numbers for the system matrix A and the corresponding residual norms for varying α.

Table 1.

Condition numbers and residual norms for the system matrix A for varying α.

α	Condition number of A	Residual norm
$1 0^{- 12}$	$1.01 \times 1 0^{11}$	$4.86 \times 1 0^{- 4}$
$1 0^{- 10}$	$1.01 \times 1 0^{11}$	$4.99 \times 1 0^{- 4}$
$1 0^{- 8}$	$9.21 \times 1 0^{10}$	$5.25 \times 1 0^{- 4}$
$1 0^{- 6}$	$2.62 \times 1 0^{13}$	$2.70 \times 1 0^{- 2}$

The stabilization strategy employed includes a regularization term added to the diagonal of A and an increased damping coefficient β. These adjustments ensure smooth and stable velocity profiles, even for larger α.

6.4.4. Implications for gradient models

The results demonstrate the role of the gradient term in regularizing the velocity distribution. For small α, the second-order term dominates, leading to nearly flat profiles in the interior of the domain, forming plateaux. As α increases, the fourth-order term suppresses sharp gradients, reducing the energy-like term and introducing smoothness.

The steep drop near $r = R_{\max}$ highlights the boundary effect imposed by the Dirichlet condition $v (R_{\max}) = 0$ . This boundary behavior reflects the dissipation or termination of velocity at the edge of the domain, consistent with physical systems requiring zero velocity at boundaries.

This behavior is typical of gradient-enhanced models, such as those used in strain-gradient elasticity, phase-field modeling, and regularized plasticity. The minimal variation in velocity profiles for the chosen α range suggests that gradient effects are secondary in this setup, dominated by the forcing term P and the damping term β.

In conclusion, the analysis highlights the stabilizing influence of the gradient term α on the velocity profiles and the energy-like term. Numerical stability is achieved through careful parameter selection, including regularization and damping. These results provide insights into the interplay between second-order and fourth-order terms in gradient models and their implications for physical and numerical behavior.

7. Norton viscous model applied to the problem of the flow in a hollow sphere under hydrostatic load

In this section, we address the scenario of a hollow sphere subject to hydrostatic loading conditions. We assume that the sphere’s material matrix follows the Norton visco-plastic rigid behavior model without yield, effectively neglecting elasticity. This problem’s solution is of interest because it provides insights that can be valuable in defining a homogenized plastic porous material at high temperatures.

Just like with plastic materials, deriving an exact solution for arbitrary loads (where D is arbitrary) poses significant challenges. While it is possible to find an approximate solution, it first necessitates the development of an approximation technique in visco-plasticity.

In this work, we will focus exclusively on the specific case of purely hydrostatic loading, represented by $D = D_{m} 1$ , where $D_{m} > 0$ indicates the sphere’s expansion. In such a case, with loads exhibiting spherical symmetry, it is feasible to arrive at an analytical solution.

Due to the incompressibility and the spherical symmetry of the problem, the velocity field can be written as:

v (r) = \frac{α}{r^{2}} e_{r}, α = b^{3} D_{m} .

(38)

As a result, the non-zero components of the strain rate and the equivalent strain rate can be found as:

\begin{matrix} {\begin{matrix} D_{rr} = \frac{d v_{r}}{d r} = - 2 D_{m} \frac{b^{3}}{r^{3}} \\ D_{θθ} = D_{ϕϕ} = \frac{v_{r}}{r} = D_{m} \frac{b^{3}}{r^{3}} . \end{matrix} \end{matrix}

(39)

The equivalent strain rate is then defined as:

D_{eq} = {(D_{rr}^{2} + D_{θθ}^{2} + D_{ϕϕ}^{2})}^{\frac{1}{2}} = | D_{rr} | = 2 D_{m} \frac{b^{3}}{r^{3}} .

(40)

By inverting the flow rule, we get:

D = \frac{3}{2} {\overset{\cdot}{ε}}_{0} {(\frac{σ_{eq}}{σ_{0}})}^{n} \frac{s}{σ_{eq}},

(41)

which implies that:

D_{eq} = {\overset{\cdot}{ε}}_{0} {(\frac{σ_{eq}}{σ_{0}})}^{n},

(42)

the verification of the latter formula is easy as the flow rule was developed such that this relation is automatically satisfied and the colinearity between D and s.

As a consequence, we get:

σ_{eq} = σ_{0} {(\frac{d_{eq}}{{\overset{\cdot}{ε}}_{0}})}^{1 / n},

(43)

which yields:

s = \frac{2}{3} \frac{d}{{\overset{\cdot}{ε}}_{0}} σ_{0} {(\frac{d_{eq}}{{\overset{\cdot}{ε}}_{0}})}^{\frac{1 - n}{n}},

(44)

which finally gives:

s = \frac{2}{3} σ_{0} {(\frac{d_{eq}}{{\overset{\cdot}{ε}}_{0}})}^{\frac{1}{n}} \frac{s}{d_{eq}} .

(45)

From there, we can deduce the components of the deviatoric tensor s read:

\begin{matrix} {\begin{matrix} s_{rr} = \frac{2}{3} σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r^{3}})}^{1 / n} \frac{d_{rr}}{| d_{rr} |} = - \frac{2}{3} σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r^{3}})}^{1 / n} \\ s_{θθ} = s_{ϕϕ} = - s_{rr} / 2 = {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r^{3}})}^{1 / n} . \end{matrix} \end{matrix}

(46)

Given the spherical symmetry of the system, the macroscopic stress tensor adopts a hydrostatic form: $Σ = Σ_{m} 1$ . Our next step is to derive an expression for $Σ_{m}$ , leveraging Hill [1]’s theory and Mandel [31]’s lemma. Notably, the virtual velocity field, denoted by $v^{⋆}$ , in this scenario aligns with the real velocity field:

v^{⋆} (x) = D_{m}^{⋆} \frac{b^{3}}{r^{3}},

(47)

which implies:

D_{rr}^{⋆} = - 2 D_{m}^{⋆} \frac{b^{3}}{r^{3}}, D_{θθ}^{⋆} = D_{ϕϕ}^{⋆} = D_{m}^{⋆} \frac{b^{3}}{r^{3}} .

(48)

As a result, we get:

{\begin{matrix} Σ : D^{⋆} = Σ_{m} D_{m}^{⋆} 1 : 1 = 3 Σ_{m} D_{m}^{⋆} \\ = (1 - f) {〈 σ : D^{⋆} (x) 〉}_{Ω - ω} \\ = (1 - f) {〈 σ_{rr} (x) D_{rr}^{⋆} (x) + 2 σ_{θθ} (x) D_{θθ}^{⋆} (x) 〉}_{Ω - ω} \\ = (1 - f) 2 D_{m}^{⋆} {〈 \frac{(σ_{rr} - σ_{θθ}) (x)}{r^{3}} 〉}_{Ω - ω} \\ = 2 (1 - \frac{a^{3}}{b^{3}}) D_{m}^{⋆} b^{3} \frac{1}{\frac{4}{3} π (b^{3} - a^{3})} \int_{a}^{b} \frac{(σ_{rr} - σ_{θθ}) (r)}{r^{3}} 4 π r^{2} d r \\ = 3 D_{m}^{⋆} \int_{a}^{b} 2 (σ_{rr} - σ_{θθ}) (r) \frac{d r}{r} \end{matrix}

(49)

and the latter equation implies that:

Σ_{m} = \int_{a}^{b} 2 (σ_{rr} - σ_{θθ}) (r) \frac{d r}{r} .

(50)

As in the case with a cylindrical shape, this result can be written as:

Σ_{m} = \int_{a}^{b} \frac{d σ_{rr}}{dr} (r) = σ_{rr} (b),

(51)

since $σ_{a} = 0$ . Note that this result is intuitively of satisfaction (even though it is not of entire satisfaction for the effective calculation of $Σ_{m}$ ). From our previous calculations:

σ_{θθ} = σ_{rr} = s_{θθ} - s_{rr} = σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r^{3}})}^{1 / n},

(52)

and thus:

\begin{matrix} {\begin{matrix} Σ_{m} = 2 σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} b^{3})}^{1 / n} \int_{a}^{b} r^{- \frac{3}{n}} \frac{d r}{r} = \frac{2}{3} σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} b^{3})}^{1 / n} \int_{a^{3}}^{b^{3}} \frac{d r^{3}}{{(r^{3})}^{\frac{1}{n} + 1}} \\ = \frac{2}{3} σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} b^{3})}^{1 / n} (- n) (b^{- \frac{3}{n}} - a^{- \frac{3}{n}}), \end{matrix} \end{matrix}

(53)

or by inverting, we get:

D_{m} = \frac{{\overset{\cdot}{ε}}_{0}}{2} {(\frac{3}{2 n} \frac{Σ_{m} / σ_{0}}{f^{- 1 / n} - 1})}^{n} .

(54)

Above, it was observed that the visco-plastic Norton model, devoid of yield, converges to the plastic von Mises model as the parameter n approaches infinity. In this asymptotic scenario, we are compelled to rediscover the expression for the yield limit of the plastic hollow sphere, which is given by $2 σ_{0} \ln \frac{a}{b} = - \frac{2}{3} σ_{0} \ln f$ .

It is noteworthy that this expression was originally derived for an internal pressure applied at $r = a$ , in contrast to the current scenario where an external tension is applied at $r = b$ . Nevertheless, our findings effortlessly demonstrate that the limit loads remain consistent for both loading conditions.

Subsequently, let us analyze the limit of the expression $n (f^{- 1 ∕ n} - 1)$ as n tends towards infinity:

n (f^{- 1 ∕ n} - 1) = n (e^{- \frac{\ln f}{n}} - 1) = n (1 - \frac{1}{n} \ln f - 1 + O (\frac{1}{n^{2}})) \to - \ln f .

In the same limit, as $n \to \infty$ , we observe that ${(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}})}^{1 ∕ n} \to 1$ . Consequently, for the limit $n \to \infty$ , we find:

Σ_{m} \to - \frac{2}{3} σ_{0} \ln f .

This result aligns with our expectations.

8. Solution of the hollow sphere problem including gradient effects

In this section, we extend the analysis of the hollow sphere problem under hydrostatic loading by incorporating high-gradient effects. This extension is particularly relevant in capturing localized deformation phenomena, size-dependent behavior, and enhanced mechanical stability, which are not accounted for in traditional visco-plastic models. Our goal is to derive the modified expressions for stress and strain distributions and analyze the physical implications of these gradient-enhanced contributions.

8.1. Incorporating high-gradient terms in stress components

The radial and hoop stresses, originally derived using a Norton visco-plastic model, are modified to include high-gradient contributions as follows:

\begin{matrix} σ_{rr} & = σ_{rr}^{visco-plastic} + η \frac{d^{3} u}{d r^{3}}, \\ σ_{θθ} & = σ_{θθ}^{visco-plastic} + η (\frac{1}{r} \frac{d^{2} u}{d r^{2}} - \frac{2}{r^{2}} \frac{du}{dr} + \frac{2 u}{r^{3}}) . \end{matrix}

(55)

where η is a material parameter representing the influence of high-gradient effects and $u (r)$ is the radial displacement. Here, $σ_{rr}^{visco-plastic}$ and $σ_{θθ}^{visco-plastic}$ denote the visco-plastic stress components derived from the original model. Further justifications of this “ad hoc” choice of stress modeling with physical and mechanical insights are provided in Appendix 2.

8.2. Modified mechanical equilibrium equation

The mechanical equilibrium equation for a hollow sphere under spherical symmetry, incorporating high-gradient contributions, is given by:

\begin{matrix} \frac{d}{dr} (σ_{rr}^{visco-plastic} + η \frac{d^{3} u}{d r^{3}}) & + \frac{2}{r} [σ_{rr}^{visco-plastic} + η \frac{d^{3} u}{d r^{3}} \\ - (σ_{θθ}^{visco-plastic} + η (\frac{1}{r} \frac{d^{2} u}{d r^{2}} - \frac{2}{r^{2}} \frac{du}{dr} + \frac{2 u}{r^{3}}))] = 0 . \end{matrix}

This equation ensures mechanical equilibrium by balancing the modified radial and hoop stresses, while accounting for the influence of high-gradient terms.

8.3. Simplified expressions for $Σ_{m}$

To derive the macroscopic stress $Σ_{m}$ , we use Hill’s averaging method, leading to:

\begin{matrix} Σ_{m} & = \int_{a}^{b} 2 (σ_{rr} - σ_{θθ}) \frac{dr}{r} . \end{matrix}

Substituting the modified stress components yields:

\begin{matrix} Σ_{m} & = \int_{a}^{b} 2 [(σ_{rr}^{visco-plastic} + η \frac{d^{3} u}{d r^{3}}) \\ - (σ_{θθ}^{visco-plastic} + η (\frac{1}{r} \frac{d^{2} u}{d r^{2}} - \frac{2}{r^{2}} \frac{du}{dr} + \frac{2 u}{r^{3}}))] \frac{dr}{r} . \end{matrix}

(56)

This expression can be further simplified by separating the visco-plastic and gradient terms, yielding:

\begin{matrix} Σ_{m} \approx & - 2 σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r_{avg}^{3}})}^{\frac{1}{n}} \ln (\frac{b}{a}) \\ + η \int_{a}^{b} [\frac{d^{2} u}{d r^{2}} (\frac{1}{r^{2}} - \frac{1}{r} \frac{d}{dr}) + \frac{2}{r^{3}} (\frac{du}{dr} - 2 u)] dr \\ + {η \frac{d^{2} u}{d r^{2}} \frac{1}{r} |}_{a}^{b} . \end{matrix}

(57)

The high-gradient terms provide enhanced mechanical stability by smoothing out stress and strain gradients, particularly near boundaries. This regularization effect reduces localization phenomena such as shear bands and promotes a more distributed flow of visco-plastic deformation. In addition, the gradient contributions introduce size-dependent behavior, making the model particularly relevant for small-scale or thin-walled structures. The dependence of the stress distribution on the radial coordinate highlights the importance of considering high-order terms in applications involving complex loading and geometrical constraints.

8.4. Simplified higher-order terms and physical interpretation

To further analyze the impact of high-gradient effects on the macroscopic stress $Σ_{m}$ in the hollow sphere model, we proceed by incorporating the simplified higher-order terms into our formulation and exploring their physical implications.

8.4.1. Simplified expression for $Σ_{m}$

Recall that the macroscopic stress $Σ_{m}$ is given by:

\begin{matrix} Σ_{m} \approx & - 2 σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r_{avg}^{3}})}^{\frac{1}{n}} \ln (\frac{b}{a}) \\ + \int_{a}^{b} η [\frac{d^{2} u}{d r^{2}} (\frac{1}{r^{2}} - \frac{1}{r} \frac{d}{dr}) + \frac{2}{r^{3}} (\frac{du}{dr} - 2 u)] dr \\ + {η \frac{d^{2} u}{d r^{2}} \frac{1}{r} |}_{a}^{b} . \end{matrix}

(58)

8.4.2. Simplification of the high-gradient integral term

We first consider the high-gradient integral term:

\int_{a}^{b} η [\frac{d^{2} u}{d r^{2}} (\frac{1}{r^{2}} - \frac{1}{r} \frac{d}{dr}) + \frac{2}{r^{3}} (\frac{du}{dr} - 2 u)] dr .

(59)

To simplify this term, we assume that the displacement field $u (r)$ varies smoothly, such that higher derivatives are small near the boundaries. In addition, we approximate the slowly varying term $(\frac{1}{r^{2}} - \frac{1}{r} \frac{d}{dr})$ over the interval $[a, b]$ .

By integrating by parts and considering boundary terms, we have:

\int_{a}^{b} η \frac{d^{3} u}{d r^{3}} \frac{dr}{r} \approx {η \frac{d^{2} u}{d r^{2}} \frac{1}{r} |}_{a}^{b} + \int_{a}^{b} η \frac{d^{2} u}{d r^{2}} \frac{dr}{r^{2}} .

(60)

8.4.3. Final approximate expression for $Σ_{m}$

Combining the contributions, the approximate expression for $Σ_{m}$ becomes:

\begin{matrix} Σ_{m} \approx & - 2 σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r_{avg}^{3}})}^{\frac{1}{n}} \ln (\frac{b}{a}) \\ + η {\frac{d^{2} u}{d r^{2}} |}_{avg} [\int_{a}^{b} (\frac{1}{r^{2}} - \frac{1}{r} \frac{d}{dr}) dr] \\ + {η \frac{d^{2} u}{d r^{2}} \frac{1}{r} |}_{a}^{b} . \end{matrix}

(61)

8.5. Asymptotic analysis of the modified stress equation

In this section, we investigate the asymptotic behavior of the modified stress parameter, $Σ_{m}$ , derived from the governing equations. The expression for $Σ_{m}$ is given as in the previous equation.

This equation combines contributions from logarithmic terms, grading effects encapsulated in the integral, and boundary terms. Here, we analyze the asymptotic behavior under the assumption of $b ≫ a$ , which is typical for many practical applications.

8.5.1. Analysis of dominant terms

8.5.2. Logarithmic contribution

The first term,

- 2 σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r_{avg}^{3}})}^{\frac{1}{n}} \ln (\frac{b}{a}),

dominates as $b ≫ a$ due to the logarithmic growth of $\ln (b ∕ a)$ . The prefactor, which includes the material parameters $D_{m}$ , ${\overset{\cdot}{ε}}_{0}$ , and $r_{avg}$ , modulates this growth. For large n, the scaling factor ${(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r_{avg}^{3}})}^{1 ∕ n}$ approaches unity, simplifying the term.

8.5.3. Integral contribution

The integral term in $Σ_{m}$ involves a grading term,

\int_{a}^{b} (\frac{1}{r^{2}} - \frac{1}{r} \frac{d}{dr}) dr .

Under the assumption of a smooth displacement field $u (r)$ , the second term $\frac{1}{r} \frac{d}{dr}$ becomes negligible compared to $\frac{1}{r^{2}}$ . Evaluating the integral gives:

\int_{a}^{b} \frac{1}{r^{2}} dr \approx \frac{1}{a} - \frac{1}{b},

which simplifies further to $\frac{1}{a}$ as $b ≫ a$ . Consequently, this contribution scales with:

η {\frac{d^{2} u}{d r^{2}} |}_{avg} \frac{1}{a},

where ${\frac{d^{2} u}{d r^{2}} |}_{avg}$ represents the average curvature of the displacement field.

8.5.4. Boundary contributions

The boundary term,

{η \frac{d^{2} u}{d r^{2}} \frac{1}{r} |}_{a}^{b},

captures the effects of grading at the boundaries. Explicitly, this term evaluates to:

η (\frac{1}{b} \frac{d^{2} u}{d r^{2}} |_{r = b} - \frac{1}{a} \frac{d^{2} u}{d r^{2}} |_{r = a}) .

For $b ≫ a$ , the contribution from $r = a$ dominates, simplifying this term to:

- η \frac{1}{a} \frac{d^{2} u}{d r^{2}} |_{r = a} .

8.6. Asymptotic expression

Combining the dominant contributions, the asymptotic expression for $Σ_{m}$ as $b ≫ a$ is:

\begin{matrix} Σ_{m} \sim & - 2 σ_{0} {(2 \frac{D_{m}}{{\overset{\cdot}{ε}}_{0}} \frac{b^{3}}{r_{avg}^{3}})}^{\frac{1}{n}} \ln (\frac{b}{a}) \\ + \frac{η}{a} ({\frac{d^{2} u}{d r^{2}} |}_{avg} - \frac{d^{2} u}{d r^{2}} |_{r = a}) . \end{matrix}

(62)

8.7. Impact of the grading term

The grading term, encapsulated in the integral and boundary contributions, adds a correction to the logarithmic behavior. Its influence depends on:

The material property η, which represents the effect of viscosity or grading.

The curvature of the displacement field ( $\frac{d^{2} u}{d r^{2}}$ ) at the boundaries and over the interval $[a, b]$ .

The scaling of a and b, with the grading term becoming relatively less significant as $b ∕ a \to \infty$ .

This analysis highlights the interplay between the dominant logarithmic term and the grading correction, emphasizing the importance of boundary effects and material properties in the asymptotic behavior of $Σ_{m}$ .

Figure 13 correctly captures the logarithmic behavior of $Σ_{m}$ as a function of $b ∕ a$ . The linear decrease with respect to $\ln (a ∕ b)$ indicates that the logarithmic term is indeed dominating the grading contribution.

Figure 13.

Logarithmic behavior of $Σ_{m}$ as a function of $b ∕ a$ .

We analyze the boundary effects arising from the gradient term in a stress equation, particularly focusing on the contribution:

Boundary Effect = η \frac{d^{2} u}{d r^{2}} \frac{1}{r},

evaluated at $r = a$ . The objective is to examine how the boundary term behaves with varying a, which represents the inner radius of a cylindrical domain, while keeping the outer radius b fixed. This analysis provides insights into the localized effects of gradients near boundaries in materials with graded properties.

The boundary effect was implemented in MATLAB with the following setup:

The parameter η represents the grading coefficient and scales the contribution of the gradient term.

The displacement curvature, $\frac{d^{2} u}{d r^{2}}$ , was modeled as a smooth function of $r = a$ :

\frac{d^{2} u}{d r^{2}} |_{r = a} = 0.1 + 0.05 \cos (a),

allowing us to simulate spatial variations in the curvature.

The boundary effect was computed for a range of a values, ensuring $a < b$ to maintain physical consistency.

The fixed outer radius was set to $b = 10$ , and the inner radius a was varied from 1 to $b - 0.1$ .

The behavior of the boundary effect was analyzed with respect to $r = a$ , producing the plot shown in Figure 14.

Figure 14.

Boundary effects arising from the gradient term in a stress equation.

For small values of a, the term $\frac{1}{a}$ dominates, leading to a sharp increase in the boundary effect. This reflects the mathematical nature of $\frac{1}{r}$ , which grows rapidly as $r \to 0$ . Physically, this suggests that the boundary effect is highly localized and significant near small radii. In gradient-enhanced material models, this corresponds to a higher contribution of curvature effects near smaller boundaries.

As a increases, the boundary effect diminishes due to the $\frac{1}{a}$ term. This is consistent with gradient models, where the contribution of curvature becomes less pronounced farther from small-scale boundaries. The grading term essentially smoothens the transition.

For larger values of a, small oscillations in the boundary effect are observed. These are due to the smooth variation in the modeled curvature, $\frac{d^{2} u}{d r^{2}}$ , which was chosen as $0.1 + 0.05 \cos (a)$ . Physically, this suggests that the material properties or displacement curvature still have subtle effects even at larger radii, reflecting inherent variations in the material’s grading.

The magnitude of the boundary effect is directly proportional to the grading coefficient η. Larger values of η would amplify the boundary effect, making it a more significant contributor to the overall stress response. Conversely, smaller values of η reduce the boundary contribution, emphasizing the dominance of other terms in the stress model.

8.8. Physical interpretation of the high-gradient terms

8.8.1. Regularization and stability

The inclusion of high-gradient terms, represented by the parameter η, introduces additional stiffness into the system. This regularization effect smoothens abrupt changes in the stress and strain fields, particularly near the boundaries of the hollow sphere. Such smoothing is critical in preventing localized deformation phenomena, such as the formation of sharp stress discontinuities or shear bands.

8.8.2. Size-dependent behavior

The dependence of high-gradient terms on the radial coordinate r and the shell thickness $(b - a)$ indicates a size-dependent mechanical response. Smaller spheres or thin-walled shells exhibit more pronounced gradient effects, altering the stress distribution compared to larger structures. This behavior aligns with observations in micro- and nano-scale applications, where traditional continuum mechanics models may inadequately capture such size-dependent effects.

9. Conclusion

This paper introduces a second-gradient visco-plastic framework to address limitations in classical visco-plastic models for size-dependent and localized deformation behaviors. The main contributions of this work are as follows:

Theoretical advancements: We incorporate second-gradient terms into the governing equations of visco-plasticity, enabling the capture of localized phenomena such as stress diffusion and shear banding. This approach smoothens stress and velocity discontinuities near rigid-visco-plastic boundaries, which classical models fail to address.

Numerical stability: Our model demonstrates improved numerical stability, avoiding oscillations and mesh dependency through higher-order gradient corrections. The implementation is validated using finite difference schemes, ensuring robustness and consistency.

Practical applications: We provide analytical and numerical solutions for two benchmark problems:

– The toothpaste flow in a cylindrical tube.

– The hollow sphere under hydrostatic loading.

The inclusion of second-gradient effects refines the stress and velocity profiles, yielding new insights into size-dependent material behaviors.

Comparison with classical models: The proposed framework is shown to outperform classical visco-plastic models in capturing smooth transitions and localized phenomena. In addition, the model reduces to classical behavior when second-gradient effects are negligible ( $α = 0$ ).

Material design implications: The second-gradient parameter (α) emerges as a key factor in influencing stress localization and material response, offering a powerful tool for designing materials with tailored mechanical properties.

Our work bridges the gap between classical visco-plasticity and modern requirements for capturing size-dependent phenomena, with implications for advanced manufacturing, material engineering, and microstructural analysis.

Footnotes

Appendix 1

Appendix 2

Appendix 3

Appendix 4 Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Koffi Enakoutsa

References

Hill

The mathematical theory of plasticity. Oxford: Oxford University Press, 1950.

Drucker

Prager

Soil mechanics and plastic analysis or limit design. Q Appl Math 1952; 10(2): 157–165.

Norton

FH.

The creep of steel at high temperatures. New York: McGraw-Hill, 1929.

Kachanov

LM.

Time of the rupture process under creep conditions. Izv Akad Nauk SSSR 1958; 8: 26–31.

Lee

EH.

Visco-plasticity. In: Flugge

(ed.) Handbook of engineering mechanics, p. 29. New York: McGraw-Hill, 1976.

Chen

High temperature creep of nickel-based superalloys. Mater Sci Eng A 1989; 110(1): 83–92.

Hodge

IM.

The viscosity of metals at elevated temperatures. J Appl Phys 1955; 26(7): 844–848.

Cocks

ACF.

A model for the behavior of hollow spheres under hydrostatic loading. Int J Solids Struct 1991; 27(4): 461–473.

Gurson

AL.

Continuum theory of ductile rupture by void nucleation and growth. J Eng Mater Technol 1977; 99: 2–15.

10.

Rice

Tracey

DM.

On the ductile enlargement of voids in triaxial stress fields. J Mech Phys Solids 1969; 17: 201–217.

11.

Norton

FH.

The creep of steel at high temperatures. Trans ASME 1929; 51: 607–635.

12.

Orowan

Creep of metals under stress. J Inst Met 1948; 74: 122–133.

13.

McLean

MJ.

Viscous flow in materials with non-ideal microstructures. Metall Trans A 1982; 13: 1137–1145.

14.

Abd El-Rehim

. Effect of grain size on the primary and secondary creep behavior of Sn–3 wt.% Bi alloy. J Mater Sci 2008; 43: 1444–1450.

15.

Turner

McNeil

Viscoelastic models for polymers under high temperature loading. Polym Eng Sci 2001; 41: 1333–1340.

16.

Esposito

Bonora

Constitutive modeling and computational analysis of creep behavior in single crystal nickel-based superalloys. In: Proceedings of 12th International Conference on Creep and Fracture of Engineering Materials and Structures, Kyoto, Japan, 2012.

17.

Run-Zi

Lei

Takamoto

Hideo

Xian-Cheng

Shan-Tung

TU.

On multiaxial creep-fatigue considering the non-proportional loading effect: Constitutive modeling, deformation mechanism, and life prediction. Int J Plast 2022; 155: 103337.

18.

Riedel

Fracture at high temperatures. Cham: Springer, 2005.

19.

Argon

AS.

Creep and anelasticity of aluminum. J Appl Phys 1975; 46(2): 303–314.

20.

Byrne

Ling

FF.

On the creep behavior of porous materials. New York: Pergamon, 1977.

21.

Barnes

Hutton

Walters

An introduction to rheology. Amsterdam: Elsevier, 1989.

22.

Denn

MM.

Extrusion instabilities and wall slip. Annu Rev Fluid Mech 2001; 33: 265–287.

23.

Bird

Armstrong

Hassager

Dynamics of polymeric liquids. Hoboken, NJ: Wiley, 1976.

24.

Bingham

EC.

Fluidity and plasticity. New York: McGraw-Hill, 1922.

25.

Macosko

CW.

Rheology: Principles, measurements, and applications. Hoboken, NJ: Wiley, 1994.

26.

Mewis

Wagner

NJ.

Colloidal suspension rheology. Cambridge: Cambridge University Press, 2011.

27.

Chhabra

Richardson

JF.

Non-Newtonian flow and applied rheology. Oxford: Butterworth-Heinemann, 2006.

28.

Nguyen

Huang

HD.

Flow of non-Newtonian fluids in pipes. Chem Eng Sci 1983; 38(9): 1357–1364.

29.

Laun

HM.

The viscosity of gels and pseudoplastic fluids. J Non Newtonian Fluid Mech 1986; 20(3): 263–275.

30.

Oldroyd

JG.

On the formulation of non-Newtonian viscosity. Proc R Soc Lond A Math Phys Sci 1950; 200: 523–541.

31.

Mandel

Contribution théorique à l’étude de l’écrouissage et des lois d’écoulement plastique. In: Proceedings of the eleventh international congress of applied mechanics, Munich, 30 August–5 September, 1964, pp. 502–509. Cham: Springer.

Gradient-enhanced visco-plastic models: Advancing localization and deformation analysis in toothpaste flow and hollow sphere deformation

Abstract

Keywords

1. Introduction

1.1. Theoretical background and context

1.2. The Norton viscous model

1.3. The toothpaste flow problem

1.4. Gradient-enhanced visco-plastic models

2. Elasto-visco-plastic response under uniaxial loading

3. Closed-form solution for toothpaste flow with an elasto-visco-plastic model

4. Norton viscous model applied to cylindrical toothpaste flow

5. Gradient-enhanced toothpaste flow analysis

5.1. Governing equations

5.1.1. Stress decomposition

5.1.2. Linear momentum balance

5.1.3. Strain rate relationship

5.1.4. Stress in the visco-plastic zone

5.1.5. Boundary layer correction near r = R 0

5.2. Solution and trial function

5.3. Numerical implementation

5.3.1. Finite difference discretization

5.3.2. Matrix formulation

5.3.3. Numerical solution procedure

5.4. Results and discussion

5.4.1. Velocity profile near r = R 0

5.5. Comparison to classical visco-plastic models

5.6. Velocity profile across the entire domain

5.6.1. Rigid zone ( r R 0 )

5.6.2. Transition at r = R 0

5.6.3. Visco-plastic zone ( r ≥ R 0 )

5.6.4. Influence of the second-gradient parameter α

5.7. Comparison with classical visco-plastic models

5.8. Physical implications

5.9. Numerical stability and accuracy

6. Additional numerical results and discussion

6.1. Simulation setup and matrix conditioning

6.2. Velocity profiles and residuals

6.2.1. Raw velocity profiles

6.2.2. Normalized velocity profiles

6.2.3. Residuals and normalized residuals

6.3. Discussion on gradient models

6.3.1. Contour plot of velocity profiles

6.4. Analysis of velocity profiles and gradient effects

6.4.1. Velocity profiles for varying α

6.4.2. Energy-like term vs. α

6.4.3. Numerical stability

6.4.4. Implications for gradient models

7. Norton viscous model applied to the problem of the flow in a hollow sphere under hydrostatic load

8. Solution of the hollow sphere problem including gradient effects

8.1. Incorporating high-gradient terms in stress components

8.2. Modified mechanical equilibrium equation

8.3. Simplified expressions for Σ m

8.4. Simplified higher-order terms and physical interpretation

8.4.1. Simplified expression for Σ m

8.4.2. Simplification of the high-gradient integral term

8.4.3. Final approximate expression for Σ m

8.5. Asymptotic analysis of the modified stress equation

8.5.1. Analysis of dominant terms

8.5.2. Logarithmic contribution

8.5.3. Integral contribution

8.5.4. Boundary contributions

8.6. Asymptotic expression

8.7. Impact of the grading term

8.8. Physical interpretation of the high-gradient terms

8.8.1. Regularization and stability

8.8.2. Size-dependent behavior

9. Conclusion

Footnotes

Appendix 1

Appendix 2

Appendix 3

Appendix 4

Funding

ORCID iD

References

5.1.5. Boundary layer correction near $r = R_{0}$

5.4.1. Velocity profile near $r = R_{0}$

5.6.1. Rigid zone ( $r R_{0}$ )

5.6.2. Transition at $r = R_{0}$

5.6.3. Visco-plastic zone ( $r \geq R_{0}$ )

8.3. Simplified expressions for $Σ_{m}$

8.4.1. Simplified expression for $Σ_{m}$

8.4.3. Final approximate expression for $Σ_{m}$