A thermodynamically consistent variational framework for non-Newtonian fluids with evolving internal variables

Non-Newtonian fluids exhibit complex rheological behavior characterized by an effective viscosity that evolves both in space and time (Damanik et al., 2010; Toulopoulos, 2023; Janela et al., 2010; Galdi, 2008; Hron et al., ,2003; Lee et al., ,2018; McClure and Kang, 2017; Galdi et al., ,2008). Unlike Newtonian fluids, whose viscosity remains constant, the effective viscosity of non-Newtonian fluids depends on local flow conditions, such as the shear rate, accumulated strain and flow history. Classical constitutive models often express these effects using empirical algebraic stress–strain relationships, such as power-law or Carreau-type models (Galdi, 2008; Carreau, 1972; Málek et al., 2002;  Hirn 2013). However, these approaches typically define viscosity as a function of spatial variables only, without accounting for temporal variations. Furthermore, defining a distinct algebraic relation restricts the function space of the velocity, limiting its ability to adapt to evolving conditions. As a result, such models are inherently limited in accurately capturing the behavior of non-Newtonian fluids whose viscosity evolves over time.

Recent numerical studies on non-Newtonian flow problems have proposed various strategies to accurately handle the nonlinear and time-dependent behavior of viscosity. For instance, Deteix and Yakoubi (2019) introduced a shear-rate projection technique that enhances the numerical stability of nonlinear viscous terms, while Barrenechea et al. (2024) developed stabilized implicit–explicit (IMEX) schemes for incompressible flows with variable viscosity, achieving both robustness and computational efficiency. Moreover, Schussnig et al. (2022) proposed a time-splitting framework for flows with nonlinear viscosity that preserves consistent boundary conditions while allowing the use of equal-order finite element pairs with improved computational efficiency. In addition, recent works have extended numerical formulations to capture complex constitutive responses in practical settings, including anisotropic viscosity and transient evolution.

In a more application-oriented context, Kumar and Mondal (2024) investigated the mixing of inelastic non-Newtonian fluids under inlet swirl conditions, demonstrating how anisotropic viscous structures interact with flow instabilities and mixing efficiency in high-fidelity simulations. Furthermore, in Garg et al. (2025), a perspective on the current state of non-Newtonian fluid dynamics is done. As highlighted in their review, such current approaches inherently lack the ability to incorporate time-dependent effects, internal-variable evolution or structural remodeling mechanisms in a thermodynamically consistent manner. An interesting class of models are of thixotropic type. Traditionally, they incorporate a structure parameter to mimic the time-dependent rebuilding and breakdown of microstructure under shear (Barnes, 1997; Mewis and Wagner, 2009; Larson and Wei, 2019). In these formulations, viscosity is expressed as a function of this structure parameter, whose evolution is typically prescribed through empirical rate equations; see e.g. Cunha et al., 2020; Siqueiraand et al.,2020. While such models provide a phenomenological way to include temporal effects, their evolution laws are not derived from thermodynamic principles, and thus, they do not offer a fully consistent energy-based description of viscosity changes.

To address (some of) these limitations, a variational modeling framework based on Hamilton’s principle has recently been introduced by Junker and Wick (Junker and Wick, 2025). In their work, an internal variable was incorporated into the space–time Hamiltonian formulation to describe the evolution of viscosity. The resulting model, referred to as the Type 1 model, captures both shear-thinning and shear-thickening behavior through a scalar internal variable that interacts with the displacement and velocity gradients (Maugin and Muschik, 1994; Berezovski, 2022). While the Type 1 model has successfully captured the behavior of non-Newtonian fluids, general non-Newtonian fluids can exhibit a wide variety of flow patterns and complex viscosity variations.

To describe the previously mentioned diverse behaviors in a broader context, this study introduces a second model, denoted as Type 2. In this model, the internal variable evolves based on the local strain energy, enabling the description of viscosity evolution through a physical mechanism distinct from that of Type 1. Consequently, the aim of this work is to present for both types a thermodynamically consistent and extensible modeling framework for non-Newtonian fluid flow, incorporating the effects of time-dependent and spatially varying viscosity through internal variables. By systematically deriving the governing equations via Hamilton’s principle of stationary action (Junker and Wick, 2025; Junker and Balzani, 2021), we obtain naturally a coupled space–time formulation. Here, it should be emphasized that the term “space–time formulation” in this context refers to the variational structure of the governing equations derived from Hamilton’s principle, which accounts for physical processes evolving both in space and time. In the present paper, it does not denote a numerical space–time finite element discretization in which the shape functions depend on both spatial and temporal coordinates; see e.g. Bazilevs et al., 2013; Tezduyar and Takizawa, 2025; Tezduyar et al., ,2006; Hughes and Hulbert, 1988; Langer and Steinbach, 2019; Schafelner, 2021; Endtmayer et al., ,2024; Tezduyar et al., ,1992; Behr, 2008. Rather in our work, the classical Rothe method is used: finite differences in time and finite elements in space. However, we notice that a Galerkin discretization with discontinuous finite elements in time results into sequential time-stepping related to classical finite difference time-stepping schemes. But still, this has the advantage from a theoretical viewpoint to have at hand the entire Galerkin technology for numerical analyses, a posteriori error estimation, and adaptivity. The proposed framework is validated through benchmark simulations in channel and cylinder flow. Therein, in extension to prior work, besides using the newly developed model Type 2, also two- and three-dimensional spatial configurations are considered.

The outline of this paper is as follows. In Section 2, our non-Newtonian fluid flow model is recapitulated and extended to Type 1 and Type 2. Then, the space–time weak formulation is derived, and the discretization is discussed. Next, in Section 3, several numerical simulations for both channel flow and flow around a cylinder are carried out. Finally, in Section 4, we conclude our discussions and provide conclusions.

In this section, we introduce our models, governing discretization and solution schemes. The key extension to prior work is the introduction of a Type 2 model.

To model fluids with a viscosity that evolves both over space and time, we adopt a thermodynamically motivated variational framework based on Hamilton’s principle of stationary action. Central to this approach is the introduction of an internal variable $γ:(0,T)→R$⁠, where $T>0$ being the end time value, which governs the dynamic evolution of the effective viscosity. Instead of prescribing an algebraic relation between stress and strain rate, the proposed method derives a differential equation for $γ$ from the stationarity of the Hamiltonian functional.

The viscosity evolution is captured through a nonlinear scaling function $f(γ)$ defined as:

which satisfies $f(γ)→1$ as $γ→−∞$ (corresponding to shear-thickening behavior), and $f(γ)→0$ as $γ→+∞$ (shear-thinning behavior). The reference state $f(0)=0.5$ provides a baseline viscosity for unaltered flow. The function $f(γ)$ used in this study is chosen to ensure both thermodynamic consistency and physical interpretability. The expression $f(γ)=11+exp(γ)$ represents a monotonically decreasing sigmoidal function with the following characteristics. First, since $0<f(γ)<1$⁠, it naturally guarantees the physical upper and lower bounds of the viscosity. Second, it enables a unified and continuous representation of both shear-thinning and shear-thickening behaviors within a single functional form. Third, the sign of its derivative, $f′(γ)<0$⁠, ensures that the evolution of viscosity always proceeds in the direction of energy dissipation, thereby satisfying the second law of thermodynamics. Moreover, due to its smooth and sigmoidal nature, the viscosity response to variations in the internal variable $γ$ remains continuous and stable, providing a numerically robust and physically consistent description of viscosity evolution. Consequently, the adopted form of $f(γ)$ can be regarded as the simplest thermodynamically consistent formulation that simultaneously satisfies the physical bounds of viscosity, thermodynamic irreversibility, unified representation of shear-dependent behavior and mathematical stability.

The primary objective of this study is to establish a unified variational framework that derives the spatial and temporal evolution of viscosity directly from physical state variables, rather than relying on conventional phenomenological non-Newtonian fluid models such as the Carreau, power-law or Bingham formulations. To this end, we build upon the internal-variable-based fluid model proposed by Junker and Wick (Junker and Wick, 2025), which is formulated within the space–time Hamiltonian principle. Within this theoretical foundation, we design distinct forms of free energy potentials according to the underlying physical mechanisms driving the viscosity evolution. Accordingly, we consider two prototype models, each based on a distinct form of free energy potential: Type 1 and Type 2.

Type 1: Interaction-based viscosity model.

In this formulation, the internal variable evolves in response to the interaction between the symmetric gradients of the displacement $u$ and velocity $v$ (Junker and Wick, 2025). The free energy potential is defined as:

and the evolution law for $γ$ reads as non-linear ordinary differential equation:

Here, $∇symw:=12(∇w+(∇w)⊤)$ denotes the symmetric gradient of a vector field $w$⁠. The parameter $μ˜$ is a viscous scaling factor that controls the magnitude of the stress contribution. The constant c represents a linear relaxation coefficient for the internal variable $γ$⁠, while $η$ denotes a viscosity-like parameter governing the rate of $γ$’s temporal evolution.

Type 2: Strain-energy-based viscosity model.

The main objective of this work is the derivation of a new model. Here, we propose that the viscosity evolution is primarily influenced by the magnitude of the strain in the displacement field, reflecting the elastic strain energy stored in the material at a given state. Then, the free energy is given by:

where $||·||:=||·||L2$⁠. The corresponding evolution equation is:

Here, $E˜$ denotes a scaling elastic modulus that weights the contribution of the strain energy $||∇symu||2$ to the free energy potential.

Differences between the two models.

These two models share a common structure in which the internal variable $γ$ modifies the effective viscosity through $f(γ)$⁠, while a linear damping term regulates its temporal variation. However, the physical quantity that drives the evolution differs between them: Type 1 emphasizes the interaction between velocity and displacement, whereas Type 2 focuses on strain energy, which seems a trivial change, but results in completely different behavior in theory and numerical simulations. More specifically, this formulation provides a strain-driven mechanism for the evolution of viscosity and is particularly suitable for modeling materials where viscosity is sensitive to the strain energy. The entire proposed framework with both models represent a wide range of non-Newtonian fluid behaviors observed in unsteady flows within a unified variational structure.

The governing equations of the non-Newtonian fluid models are derived by applying Hamilton’s principle of stationary action to a space–time functional (Junker and Wick, 2025; Junker and Balzani, 2021). The total Hamiltonian is defined as:

where $I=(0,T)$ denotes the time interval of interest. For the following, let $Ω ⊂ Rd,d=2,3$ be the spatial domain. The product $Ω×I$ is known as the space–time domain. Then:

• $K$ denotes the kinetic energy:

• $G$ denotes the Gibbs (potential) energy:

where $b★$ and $−p★n$ are, respectively, the body force vector and the traction vector for the fluid.

• $D$ represents the energy dissipation (from fluid viscosity and internal variable evolution):

where the dissipation functions are defined as:

• $C$ denotes the constraint enforcing incompressibility:

with the Lagrange multiplier $p˜$⁠.

Hamilton’s principle requires the variation of the Hamiltonian to vanish:

which implies that each individual variation must vanish:

We now derive the weak forms for each prototype model. Importantly, these weak forms are automatically in a space–time format, thanks to the Hamilton principle.

The weak form with respect to the displacement $u$ reads as follows after integration by parts:

The variation with respect to the velocity $v$ yields after integration by parts the weak form:

The pressure term imposes the incompressibility constraint and follows from a variation with respect to the Lagrange multiplier as after integration by parts:

Finally, the weak form associated with the internal variable $γ$ governs its evolution and reads:

For completeness, we also state the strong form, which has the advantage that the boundary and initial conditions are explicitly visible. The resulting system represents the local form of the Type 1 model and reads as follows:

Type 1 governing equations.

Find $v:Ω×I→Rd,u : Ω×I→Rd,p˜ : Ω×I→R$ and $γ:I→R$ such that:

Boundary conditions.

On the space–time boundary $∂Ω×I$ with the decomposition $∂Ω=∂ΩD∪∂ΩN$⁠, the following conditions are imposed:

Initial conditions.

At the initial time $t=0$⁠, all fields are assumed to be at rest:

corresponding to the reference viscosity state $f(γ(0))=0.5$⁠. These equations form a coupled space–time system that simultaneously captures the nonlinear fluid dynamics and the evolution of the internal viscosity state (Junker and Wick, 2025).

The weak form of the Hamiltonian functional with respect to the displacement field $u$ is obtained after variation and integration by parts as follows:

The weak form of the Hamiltonian functional with respect to the velocity field $v$ is obtained after variation and integration by parts as follows:

The weak form with respect to the pressure variable $p˜$⁠, which enforces the incompressibility constraint, remains unchanged from the previous formulation:

Finally, the evolution of the internal variable $γ$ is described by the following weak form after variation, which includes contributions from strain energy and a linear relaxation term:

As before, for completeness, we state the strong form of the Type 2 model:

Type 2 governing equations.

Find $v : Ω×I→Rd,u : Ω×I→Rd,p˜ : Ω×I→R$ and $γ:I→R$ such that:

Boundary conditions.

On the space–time boundary $∂Ω×I$ with the decomposition $∂Ω=∂ΩD∪∂ΩN$⁠, the following conditions are imposed:

Initial conditions.

At the initial time $t=0$⁠, all fields are assumed to be at rest:

corresponding to the reference viscosity state $f(γ(0))=0.5$⁠.

The numerical solution starts from the space–time variational formulation introduced in the previous subsections. First, we notice that we closely follow Junker and Wick (2025). The principle idea is a decoupled scheme where $(v,u,p˜)$ are solved monolithically and then iteratively coupled to the solution of $γ$⁠. A classical Rothe method with finite differences in time and finite elements in space is used. The resulting global system at each time point is solved using a monolithic Newton method. Residuals and tangent matrices are computed via automatic differentiation.

In more detail, the spatial discretization uses the Taylor–Hood finite element pair, where the velocity $v$ and displacement $u$ fields are approximated by continuous quadratic quadrilateral elements (⁠$Qc2$⁠), while the pressure $p˜˙$ is approximated by continuous linear quadrilateral elements (⁠$Qc1$⁠). This choice satisfies the Ladyzhenskaya–Babuška–Brezzi (LBB) condition, ensuring numerical stability for the incompressibility constraint Link to the cited article. Furthermore, during the time integration procedure, the Newton–Raphson iterations typically converged within four iterations per time step, confirming that the proposed numerical framework operates in a stable and efficient manner. Although the governing equations of the present formulation are derived within a space–time variational framework based on Hamilton’s principle, the numerical implementation does not employ a fully space–time finite element method in which shape functions depend on both spatial and temporal coordinates. Instead, a semi-discrete approach is adopted: the Taylor–Hood shape functions depend only on the spatial variables, while the temporal discretization is carried out independently by means of implicit time-integration schemes, such as the backward Euler or Crank–Nicolson methods.

In addition to the spatial and temporal discretization described above, the treatment of the internal variable $γ$ and the chosen element integration strategy play an essential role in the numerical stability and computational efficiency of the proposed framework. Element integration is performed using standard Gaussian quadrature, and the global degrees of freedom for the velocity, displacement and pressure fields are stored at the nodal points. By contrast, the internal variable $γ$ is stored and updated locally at each Gauss integration point. This design maintains consistency of the constitutive formulation while avoiding unnecessary coupling terms that would arise if $γ$ were included as a global unknown.

A one-way coupling strategy is adopted (Junker and Wick, 2025). Specifically, the internal variable $γ$ from the previous time step is used to compute the global degrees of freedom for velocity $v$⁠, displacement $u$ and pressure $p˜˙$⁠. After the global fields have been updated, the internal variable $γ$ is independently updated at each integration point based on the newly obtained state variables. This integration-point-based storage strategy enables an accurate evaluation of the constitutive relations without introducing the internal variables as global degrees of freedom, thereby offering an implementation that is advantageous in terms of accuracy while simultaneously avoiding unnecessary global memory consumption. As a result, the present approach constitutes a well-balanced computational strategy that preserves numerical stability and faithfully reflects the physical meaning of the internal variables.

Moreover, the focus is given to the accurate reproduction of the internal-variable evolution and the physical fidelity of the constitutive formulation. For large-scale three-dimensional problems using high-performance computing (HPC), storing internal variables at the nodal points and interpolating them to the Gauss points when needed represents a viable alternative. In the present work, where quadrilateral (2D) and hexahedral (3D) elements are predominantly used and massively parallel computations are not targeted, an integration-point-based storage strategy is preferred, as it more clearly preserves the locality of the internal-variable evolution and the thermodynamic consistency of the constitutive model.

This section presents numerical simulations for the proposed non-Newtonian fluid models with internal variables. Two benchmark problems are considered, which are investigated each in $d=2$ and $d=3$⁠, thus resulting into four different configurations. The first one is a laminar flow in a simple rectangular channel, and the second one involves flow past a circular cylinder embedded within the channel (Schäfer and Turek, 1996). These examples are used to evaluate how the evolution of the internal variable influences the velocity field and viscosity distribution. Comparisons with Newtonian fluid results are provided to assess the nonlinear effects.

Boundary conditions.

Since a fully developed laminar flow is assumed, a parabolic velocity profile is prescribed at the inlet boundary. To prevent numerical instabilities that may arise from a sudden increase in inlet velocity at the initial time, a time-dependent scaling function is applied to gradually ramp up the inlet velocity. The complete inlet velocity profile $d=2$ is given by:

where H is the height of the channel. This function starts at zero when $t=0$ and asymptotically approaches the original velocity profile as $t→2$⁠, ensuring a smooth transition in the inlet velocity and improving convergence behavior in the initial transient regime of the nonlinear problem. Dirichlet conditions (non-slip conditions) are applied to the velocity and displacement on the top and bottom walls, while a zero traction condition is imposed at the outlet boundary. For $d=3$⁠, the following inlet velocity distribution is prescribed:

where H and D represent the channel height and depth, respectively.

Initial conditions.

At the initial time $t=0$⁠, all fields are assumed to be at rest, that is,

The internal variable is initialized as $γ|t=0=0$⁠, corresponding to the reference viscosity state $f(γ)|t=0=0.5$⁠. This choice represents a quiescent configuration with the nominal viscosity and ensures a stable start of the transient simulation.

In this first example, the flow behavior is analyzed by considering two representative Reynolds number regimes. The first corresponds to the low-Reynolds-number regime, for which we select $Re<100$ to examine the viscosity-dominated behavior. In this regime, the characteristics of the proposed internal-variable-based model are assessed through direct comparisons with classical non-Newtonian models such as the Carreau and power-law formulations.

The second regime represents a moderately high Reynolds number, where we adopt the representative interval of $Re>100$⁠. This case is used to investigate conditions in which inertial effects become more pronounced.

Through these two complementary regimes, we aim to consistently verify that the proposed model reproduces physically meaningful behavior across both the low-Reynolds-number (viscosity-dominated) and moderately high-Reynolds-number (inertia-influenced) conditions.

This simulation considers two-dimensional channel flow in a rectangular domain with a length of $2.5m$ and a height of $0.41m$⁠; see Figure 1. The fluid is assumed to be incompressible, with a density of $ρ=1.0kg/m3$ and a kinematic viscosity of $μ=0.005m2/s$⁠. To capture non-Newtonian behavior, an internal variable $γ$ is introduced. The evolution of $γ$ is governed by two physical constants: $c=4.412×10−8kg/(m·s2)$ and $η=10.0kg/(m·s)$⁠.

To provide a reference against classical empirical formulations, we compare the proposed internal-variable-based model with the well-established Carreau model (Carreau, 1972; Bird et al.,1986; Tanner, 2000). The Carreau model is widely used to describe non-Newtonian fluids whose viscosity varies nonlinearly with the shear rate, and it is capable of reproducing both shear-thinning and shear-thickening behaviors depending on the choice of model parameters. In this study, we adopt the commonly used simplified form of the Carreau viscosity law (with $a=2$⁠), given by:

The model parameters have the following physical interpretation:

• $μ0$⁠: zero-shear viscosity;

• $μ∞$⁠: infinite-shear viscosity;

• $λ$⁠: characteristic relaxation time governing the transition between viscosity regimes;

• n: power-law index (⁠$n<1$ for shear-thinning, $n>1$ for shear-thickening); and

• $D$⁠: $∇symv$⁠, the symmetrized velocity gradient.

In this study, the Carreau parameters used for qualitative comparison were chosen as follows. The zero-shear and infinite-shear viscosities were set to $μ0=2$ and $μ∞=0.005$⁠, respectively, while the characteristic relaxation time was taken as $λ=2$⁠. The power-law index n was selected as $n=0.2$ for the shear-thinning case and $n=1.4$ for the shear-thickening case, allowing the model to reproduce both types of non-Newtonian behavior. These parameter choices lie within a physically reasonable range and serve as reference values for assessing the qualitative similarities and differences between the proposed internal-variable framework and the Carreau viscosity law.

For the internal-variable-based formulation, the material constants are chosen as follows. For the nonlinear viscosity scaling, the model uses $μ˜=5.03947N·s/m2$⁠, while in the Type 2 model, the strain-energy contribution is scaled by $E˜=0.1027MPa$⁠. These constants were selected empirically within a physically reasonable range so as to produce qualitative viscosity variations similar to those observed in the Carreau model, including both shear-thinning and shear-thickening trends.

Based on the physical properties, the Reynolds number for a Newtonian fluid with $f(γ)=0$ is approximately $Re=82$⁠. This dimensionless quantity measures the ratio between inertial and viscous forces and is expressed as $Re=ρULμ$⁠, with the characteristic velocity fixed to $U=1$ in the present study. In channel flow, the characteristic length L corresponds to the channel height H, giving $Re=ρUHμ$⁠, while in flow past a cylinder, it is taken as the cylinder diameter D, leading to $Re=ρUDμ$⁠. Time discretization is carried out using the implicit Euler scheme with a fixed time step of $Δt=0.02s$⁠. The total simulation time is $T=15s$⁠, resulting in 750 time steps. At the converged state after 750 time steps, the velocity profiles at the cross-section located at $x=0.5$ are shown in Figure 2. This section corresponds to the fully developed region of the channel flow and provides a suitable basis for evaluating qualitative and quantitative differences between the internal-variable-based models (Type 1 and Type 2) and the Carreau model.

Under shear-thinning conditions (Type 2 comparison), the proposed internal-variable model produces a noticeably flatter velocity profile in the channel center, consistent with the plug-like behavior that arises due to the reduction of effective viscosity. The Carreau model also exhibits shear-thinning, yet its viscosity relaxation at high shear rates is comparatively milder, leading to differences in the magnitude of the central velocity gradient and the deceleration near the walls. These discrepancies can be attributed to the fact that the internal-variable formulation incorporates spatial and temporal evolution effects, whereas the Carreau model relies solely on an algebraic shear rate-dependent viscosity law.

In contrast, under shear-thickening conditions (Type 1 comparison), both models produce nearly identical parabolic velocity profiles. The increase in effective viscosity suppresses profile flattening, and the maximum centerline velocity as well as the overall distribution agree closely between the two formulations. This is expected, as the qualitative viscosity increase in shear-thickening flows is similar for both models, and the influence of internal-variable evolution becomes less pronounced in steady-state regimes.

Based on these observations, while the proposed internal-variable model exhibits a more pronounced plug-like behavior than the Carreau model under shear-thinning conditions, the two models show nearly identical behavior for shear-thickening flows. This demonstrates that the internal-variable framework is capable of reproducing a broad spectrum of non-Newtonian responses and highlights its ability to capture history-dependent effects in shear-thinning regimes.

For the nonlinear viscosity scaling, the model uses $μ˜=7.03947N·s/m2$⁠, while for the Type 2 model, the strain energy term is scaled by $E˜=0.0654MPa$⁠. These parameters were chosen empirically within a physically reasonable range to ensure numerical stability and convergence of the nonlinear solver. They act as reference scaling constants rather than fitted material properties, and small variations within the same order of magnitude do not affect the qualitative flow behavior or the conclusions of the study. The boundary conditions used in this channel flow example are summarized in Table 1. With the chosen physical parameters (⁠$U=2$⁠), the Newtonian fluid yields a Reynolds number of approximately $Re=164$⁠. Here, $Re=ρUL/μ$ is evaluated using the channel height as the characteristic length, consistent with standard definitions for internal flows. In the case of the Type 1 model, the effective Reynolds number varies within a wide range from $Re ≈0.1$ to 164 due to strong viscosity changes induced by the internal variable $γ$⁠. By contrast, the Type 2 model defines viscosity evolution governed by the strain energy, which exerts a pronounced influence on the effective Reynolds number. As a result, the Reynolds number varies approximately from $Re≈2$ to 164, with its temporal evolution being driven by deformation. This behavior highlights the fundamental differences in flow characteristics compared to the Newtonian case.

Time integration is performed using an implicit Euler scheme with a constant time step of $Δt=0.02<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠, resulting in a total of 400 steps over the simulation horizon of $T=8<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠.

Figure 3 illustrates the temporal evolution of the velocity profile [Figure 3(a)] and the viscosity function $f(γ)$ at the cross-section $x=0.5m$⁠. In the channel center, the viscosity exhibits almost no variation, whereas near the walls, $f(γ)$ increases over time, leading to localized shear-thickening [Figure 3(e)]. As a result of this viscosity behavior, the flow in the central region maintains a relatively high velocity while in the near-wall regions the increased viscosity significantly reduces the flow speed [Figure 3(b)]. Consequently, a parabolic velocity profile is observed [Figure 3(d)].

By contrast, the Type 2 model exhibits an inverse tendency. Near the walls, $f(γ)$ decreases [Figure 3(f)], resulting in a locally reduced viscosity and hence faster near-wall flow. Meanwhile, in the central region, the variation of $f(γ)$ is relatively small, and thus, a flatter velocity distribution compared to the Type 1 case [Figure 3(c)]. This behavior indicates that the gradient of the displacement field plays a significant role in the evolution of viscosity in the Type 2 model.

Figure 4 illustrates the temporal evolution of the channel flow over the entire domain obtained with the Type 1 model. The figure illustrates the temporal evolution of the velocity profile and the viscosity function $f(γ)$ at successive time steps. In the channel center, the values of $f(γ)$ remain nearly constant, indicating a low-viscosity region in which the flow maintains relatively high velocity. By contrast, near the walls, $f(γ)$ approaches unity, leading to a significant increase in the effective viscosity. This results in localized shear-thickening behavior, which suppresses the growth of velocity in the boundary layer and stabilizes the shear layer. Consequently, the flow develops into a characteristic parabolic velocity profile. These findings demonstrate that the evolution of the internal variable $γ$ directly governs the local viscosity distribution and thereby strongly influences the overall flow characteristics. In contrast to the Type 1 case, the Type 2 model (Figure 5) exhibits progressive shear-thinning near the inlet and along the channel walls. Over time, the values of $f(γ)$ decrease in the near-wall regions, leading to a pronounced reduction in viscosity and consequently faster flow close to the boundaries. This mechanism prevents the formation of a pronounced viscous boundary layer typically observed in Newtonian fluids, as the velocity gradients near the walls are relaxed. As a result, the overall velocity distribution becomes increasingly flattened across the channel. Hence, the Type 2 model demonstrates that strain-energy-driven viscosity evolution leads to accelerated near-wall flow and a more uniform velocity profile throughout the domain.

It is important to emphasize that the observed differences between the Type 1 and Type 2 models do not arise from changes in the algebraic form of a constitutive viscosity law, as is the case in classical generalized Newtonian formulations such as the power-law, Carreau or Bingham models. Instead, the crucial distinction lies in the presence of an internal variable whose evolution governs the effective viscosity in both space and time. Classical generalized Newtonian models prescribe viscosity as an instantaneous function of the local shear rate, without incorporating memory or structural evolution effects. A comparison with traditional thixotropic models further highlights the conceptual differences. Thixotropic models use a structure parameter to mimic the time-dependent breakdown and rebuilding of microstructure; however, the evolution of this parameter is usually defined through empirically chosen rate equations and is not derived from thermodynamic considerations. While such models phenomenologically capture temporal effects, they do not provide a free-energy-based or dissipation-consistent description of viscosity evolution. By contrast, the present formulation introduces a thermodynamically motivated evolution equation for the internal variable $γ$⁠, derived directly from the free energy potential and dissipation functional. This ensures that viscosity changes occur consistently with the second law of thermodynamics, allowing the effective viscosity to adapt dynamically during flow development in a physically grounded manner. Consequently, the flow features observed in the channel configuration reflect not only the current shear-rate distribution but also the accumulated evolution of $γ$⁠, leading to viscosity patterns that differ fundamentally from those obtained with algebraic constitutive models or empirically defined thixotropic formulations.

In this second example, we consider a two-dimensional channel flow with an embedded circular obstacle (Schäfer and Turek, 1996). The computational domain is identical to the previous channel flow case, with a length of 2.5 m and a height of 0.41 m. A fixed cylinder with a radius of 0.05 m is placed at the position (0.2 m, 0.2 m) within the domain. The fluid is assumed to be incompressible and isothermal, with a constant density of $ρ=1.0kg/m³$⁠. The physical parameters governing the evolution of the internal variable $γ$ are set identical to those in the channel flow configuration, namely, $c=4.412×10−8$ kg/$(m·s2)$ and $η=10.0$ kg/$(m·s)$⁠. The nonlinear viscosity scaling function $f(γ)$ forms the core of the proposed model, with the viscosity scaling coefficient set to $μ˜=9.98N·s/m2$⁠, allowing the effective viscosity to vary in both time and space. In this study, two different values of the base viscosity $μ$ are considered: 0.05 and 0.001. Time discretization is carried out using the Crank–Nicolson scheme with a fixed time step of $Δt=0.01s$⁠. The total simulation time is $T=8s$⁠, resulting in 800 time steps. The boundary conditions used in this cylinder flow example are summarized in Table 2.

For the Type 1 model, when $μ=0.02N·s/m2$⁠, the Reynolds number varies within the range of approximately 0.02 to 10.0, maintaining a fully laminar flow throughout the domain. By contrast, when $μ=0.001N·s/m2$⁠, the Reynolds number can reach up to $Re=200$⁠.

For the Type 2 model, viscosity evolution is driven by strain energy, with an additional elastic modulus $E˜=0.1MPa$ assigned. The parameters (⁠$η,μ˜,c$ and $E˜$⁠) for Type 1 and Type 2 were determined empirically, following the same rationale used in the channel flow example. This setting allows the viscosity to respond sensitively to the strain-rate distribution, enabling observation of how the evolving viscosity influences the flow field for Reynolds numbers in the range $2≲Re≤200$⁠. In particular, while in a Newtonian fluid, the well-known Kármán vortex street forms in the cylinder wake when the Reynolds number exceeds 100, the proposed non-Newtonian models (Type 1 and Type 2) may produce different and unsteady flow patterns depending on their viscosity evolution mechanisms. This provides an opportunity to analyze and compare how the temporally and spatially evolving viscosity affects vortex generation, boundary layer development and flow asymmetry. Figure 7 illustrates the temporal evolution of the flow at a location $x=0.5<private-char specific-use="loader-generated">&#xef07;</private-char>m$ (red line in Figure 6) downstream of the cylinder in the low Reynolds number regime (⁠$Remax=10$⁠) for three different fluids. In the Newtonian case [7(a)], the velocity in the central region behind the cylinder decreases, resulting in a concave velocity profile with a reduced centerline speed. By contrast, the Type 1 model exhibits an increase in viscosity {$f(γ)$ increases, [7(e)]} near the walls and a consistence in viscosity at the center over time, leading to a parabolic velocity profile characteristic of laminar channel flow, with high velocities at the center and low velocities near the walls [7(b)].

For the Type 2 model, $f(γ)$ remains nearly unchanged in the initial stage, maintaining a uniform flow. As time progresses, viscosity decreases in both the central and near-wall regions, while in the vicinity of the cylinder it decreases relatively slowly [7(f)]. Consequently, due to shear-thinning effects, the viscosity reduction in the central and near-wall regions leads to velocity increases in these areas [7(c)]. However, the relatively preserved viscosity around the cylinder suppresses acceleration in the wake center, resulting in a concave velocity profile with lower central velocities and relatively higher velocities near the sides, similar to the Newtonian case [7(d)].

Figure 8 presents the temporal evolution of drag and lift coefficients for the Newtonian fluid and the two proposed non-Newtonian models (Type 1 and Type 2). For the Type 1 model, pronounced shear-thickening occurs around the cylinder, leading to a significant local increase in viscosity. As a result, both drag and lift coefficients remain consistently higher than those of the Newtonian fluid and the Type 2 model, and eventually converge to even larger steady-state values. This behavior can be interpreted as the outcome of viscosity enhancement, which amplifies the hydrodynamic forces acting on the cylinder.

By contrast, the Type 2 model initially exhibits relatively high viscosity, resulting in larger drag and lift coefficients compared to the Newtonian case. During the transient phase, these forces further increase; however, as the internal variable $γ$ evolves, viscosity gradually relaxes and a transient response emerges, leading to a subsequent reduction in both drag and lift. This rise-and-fall pattern is not associated with flow instabilities but rather originates from the viscosity evolution mechanism of the model. Since the Reynolds number in this study is relatively low, the flow remains laminar, and the observed decrease in hydrodynamic forces is attributed to the transient viscosity relaxation rather than structural changes in the flow field. Figure 9 [see Subfigures 9(a), 9(b), 9(c), 9(d)] presents the convergence study conducted for the non-Newtonian fluid model. The upper part of the figure [Subfigures 9(a), 9(b)] illustrates the spatial convergence (mesh refinement $r1→r4$⁠), where the corresponding meshes consist of approximately 328, 662, 5,862 and 1,0328 elements, with about 5,972, 11,893, 101,337 and 178,150 degrees of freedom, respectively. The lower part [Subfigures 9(c), 9(d)] shows the temporal convergence obtained by varying the time-step size with a fixed mesh. In both cases, as the mesh is refined or the time step is reduced, the drag and lift coefficients gradually approach stable values. These results demonstrate that the proposed formulation provides grid- and time-independent solutions, thereby confirming its reliability and consistency.

Figure 10 illustrates the temporal evolution of the velocity profile at the location $x=0.5<private-char specific-use="loader-generated">&#xef07;</private-char>m$ in the high Reynolds number regime (⁠$Remax=200$⁠). For the Newtonian fluid, due to the relatively high Reynolds number, periodic vortex shedding is observed in the wake of the cylinder. As a result, oscillatory behavior in the velocity signal is clearly visible in Figure 10(a).

In the Type 1 model, the observed shear consistence in the central region and shear-thickening near the walls each contribute to flow stability in different ways. The velocity profile is shown in Figure 10(b). The shear-thickening near the walls raises viscosity, thickening the boundary layer and consequently suppressing shear-layer instabilities [Figure 10(e)]. This viscosity distribution suppresses the growth of the flow velocity; see Figure 10(c). The Type 2 $f(γ)$ field is shown in Figure 10(f). Consequently, in the wake of the cylinder, a distinct Kármán vortex street develops, similar to that observed in Newtonian fluids under high Reynolds number conditions; see Figure 10(d). Figures 11–13 depict the evolution of the flow past a circular cylinder at $Remax=200$ over time, for three distinct cases, across the entire domain. Figure 11 shows the Newtonian fluid, where wake instabilities progressively intensify over time and lead to periodic vortex shedding. As a result, a characteristic Kármán vortex street is clearly observed in the wake. By contrast, Figure 12 illustrates the Type 1 model, where strong shear-thickening is observed around the cylinder and along the channel walls. The resulting increase in local viscosity effectively suppresses wake instabilities, so that the flow maintains a stable parabolic profile over time without exhibiting periodic vortex shedding. This stabilization mechanism demonstrates how localized viscosity enhancement inhibits the development of shear layers and preserves flow symmetry in the wake.

Figure 13 depicts the Type 2 model, which exhibits the opposite trend. At the initial stage (e.g. time step 2 s), shear-thinning begins around the cylinder and walls. As time progresses, the reduction of viscosity spreads downstream into the wake region, leading to a substantial weakening of viscous damping. The diminished viscosity sharpens the velocity gradients in the shear layers, thereby amplifying Kelvin–Helmholtz-type instabilities. Consequently, the wake becomes increasingly unstable and eventually develops a Kármán vortex street similar to that of the Newtonian fluid. Similar to the channel flow configuration, the cylinder flow results in Figures 11–13 clearly show that the proposed internal-variable-based variational formulation describes viscosity evolution in a manner fundamentally different from classical empirical non-Newtonian models. Traditional generalized Newtonian models define viscosity as an instantaneous algebraic function of the local shear rate; hence, spatial and temporal variations of viscosity are determined solely by the current shear-rate field, without accounting for flow-history effects or any accumulated viscosity evolution. By contrast, the present variational framework determines viscosity through a time-dependent evolution equation for the internal variable $γ$⁠, thereby treating viscosity as a dynamic quantity rather than as a static constitutive function. As a result, the viscosity adapts spatially and temporally to local flow conditions, and this mechanism becomes clearly visible in the cylinder flow simulations. For the Type 1 model, the evolution of $γ$ induces shear-thickening around the cylinder and along the channel walls, increasing the local viscosity and suppressing wake instabilities. In the Type 2 model, accumulated shear-thinning gradually reduces the viscosity in the wake region, enhancing Kelvin–Helmholtz-type instabilities and ultimately re-establishing a Kármán vortex street similar to that observed in Newtonian fluids. Accordingly, the proposed formulation captures a wide range of unsteady flow phenomena – including time-dependent viscosity variations, instability development, transition mechanisms and flow-pattern formation – that cannot be described using classical algebraic viscosity laws, while providing a thermodynamically consistent and physically interpretable framework for modeling non-Newtonian fluids.

In this third numerical test, a three-dimensional configuration is designed. The purpose of this example is to verify that the model and numerical settings validated in the two-dimensional case can be naturally extended to three-dimensional configurations under laminar flow conditions at moderate Reynolds numbers, without introducing additional assumptions or dimension-dependent modifications. The discretization and solution algorithms are of the same type as before.

To this end, the same set of parameters used in the two-dimensional channel flow problem is retained for the three-dimensional computations. Specifically, for the Type 1 model, the viscosity-related scaling parameter is set to $μ˜=7.03947<private-char specific-use="loader-generated">&#xef07;</private-char>N·s/m2$⁠, while for the Type 2 model, the strain-energy term is scaled by $E˜=0.0654<private-char specific-use="loader-generated">&#xef07;</private-char>MPa$⁠, corresponding to a maximum Reynolds number of $Re=164$⁠. This choice allows us to confirm that the proposed model preserves structurally consistent behavior under identical physical parameters, independent of the spatial dimension. All remaining physical and numerical parameters, such as the fluid density, characteristic velocity and time-step size are kept identical to those used in the two-dimensional channel flow case. Specifically, the fluid density is set to $ρ=1.0kg/m3$⁠, the characteristic velocity to $U=2<private-char specific-use="loader-generated">&#xef07;</private-char>m/s$ and the time-step size to $Δt=0.02<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠. The physical constants are chosen as $c=4.412×10−8<private-char specific-use="loader-generated">&#xef07;</private-char>kg/(m·s2)$ and $η=10.0<private-char specific-use="loader-generated">&#xef07;</private-char>kg/(m·s)$⁠, and the kinematic viscosity is $μ=0.005<private-char specific-use="loader-generated">&#xef07;</private-char>m2/s$⁠. In total, 400 time steps are performed, corresponding to a final simulation time of $8<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠.

In the present three-dimensional example, the computational domain is a rectangular channel obtained by uniformly extending the geometric configuration of the two-dimensional channel flow problem in the depth direction (Figure 1). Specifically, the channel has a length of $L=2.5<private-char specific-use="loader-generated">&#xef07;</private-char>m$⁠, a height of $H=0.41<private-char specific-use="loader-generated">&#xef07;</private-char>m$ and a depth of $D=0.41<private-char specific-use="loader-generated">&#xef07;</private-char>m$⁠. The computational mesh is constructed by consistently extending the meshing strategy used in the two-dimensional case, and the entire domain is discretized using three-dimensional hexahedral finite elements. The entire computational domain is discretized using a total of 7,600 finite elements, resulting in an overall number of 250,271 degrees of freedom.

With regard to the boundary conditions, the same configuration as that used in the two-dimensional channel flow problem is adopted. No-slip conditions are imposed on the upper and lower walls as well as on the front and back walls of the channel, while a traction-free condition is prescribed at the outlet (see Table 3). By contrast, the inlet boundary condition is adapted to account for the geometric characteristics of the three-dimensional channel flow, where a three-dimensional velocity profile defined in both the height and depth directions is applied.

Figures 14 and 15 show, for the Type 1 model and the Type 2 model, respectively, the temporal evolution of the streamwise velocity field and the internal variable function $f(γ)$ evaluated on a cross-section taken at the mid-plane of the three-dimensional channel, corresponding to $z=0.205<private-char specific-use="loader-generated">&#xef07;</private-char>m$⁠.

All results exhibit qualitative agreement with the behavior observed in the two-dimensional channel flow case. For the Type 1 model, the spatial distribution of the internal variable leads to shear-thickening behavior in the vicinity of the walls, while a relatively shear-thinning response is observed in the central region of the channel. As a consequence of this distribution, the velocity increases in the channel core and decreases near the walls (Figure 14). By contrast, for the Type 2 model, the internal-variable-based mechanism gives rise to a predominantly shear-thinning behavior near the walls, which results in an increased velocity in the wall-adjacent regions (Figure 15).

Figure 16 presents the cross-sectional distributions of the streamwise velocity field and the internal variable at a fully developed streamwise location, $x=0.5<private-char specific-use="loader-generated">&#xef07;</private-char>m$⁠, in the three-dimensional channel flow. Beyond a mere three-dimensional visualization, the figure is designed to allow for a direct assessment of the spatial correlation between the velocity field and the internal variable at the cross-sectional level. In particular, it serves to verify whether the characteristic distribution of the internal variable and the associated flow structures observed in the two-dimensional channel flow are consistently reproduced in a three-dimensional setting. Through this cross-section-based analysis, it can be confirmed that the mechanism by which the internal variable governs the flow behavior in the proposed model remains structurally preserved despite the increase in spatial dimensionality.

The cross-sectional results of the Type 1 model shown in Figure 16 demonstrate that the distributions of the internal variable and the velocity field are formed symmetrically with respect to the channel geometry. The internal variable function $f(γ)$ attains relatively high values in the vicinity of the channel walls, indicating that shear-thickening behavior is dominant in these regions due to the no-slip boundary condition. By contrast, the central region of the channel exhibits a comparatively shear-thinning response. As a consequence of this spatial distribution of the internal variable, the streamwise velocity increases in the channel core while being reduced near the walls, resulting in a characteristic velocity profile. These results demonstrate that the Type 1 model is capable of clearly distinguishing spatially varying shear states in a three-dimensional channel flow and of simultaneously capturing shear-thickening and shear-thinning behaviors within a unified variational framework.

The cross-sectional results of the Type 2 model shown in Figure 16 reveal a decreasing trend of the internal variable function $f(γ)$ along the four walls of the three-dimensional channel, leading to a predominance of shear-thinning behavior in the wall-adjacent regions. This reduction of the internal variable results in a decrease of the effective viscosity near the walls, which in turn increases the streamwise velocity in the vicinity of the boundaries and promotes a more uniform redistribution of the velocity profile across the channel cross-section. In particular, the fact that the velocity profile is no longer concentrated in the channel core but instead exhibits a smoothly flattened shape extending toward the walls is characteristic of shear-thinning behavior in three-dimensional channel flows. These observations indicate that the wall-centered shear-thinning mechanism identified for the Type 2 model in the two-dimensional channel flow is consistently reproduced in three dimensions while preserving geometric symmetry.

This final section considers a three-dimensional flow around a circular cylinder. The same physical and numerical parameters as those used in the two-dimensional cylinder flow analysis are adopted. The maximum Reynolds number is set to $Re=200$⁠. For the Type 1 model, the viscosity-related scaling parameter is chosen as $μ˜=9.98<private-char specific-use="loader-generated">&#xef07;</private-char>N·s/m2$⁠, while for the Type 2 model, the strain-energy-related scaling parameter is set to $E˜=0.1<private-char specific-use="loader-generated">&#xef07;</private-char>MPa$⁠. All remaining physical and numerical parameters are kept identical to those used in the two-dimensional cylinder flow simulations. Specifically, the fluid density is set to $ρ=1.0kg/m3$⁠, the characteristic velocity to $U=2<private-char specific-use="loader-generated">&#xef07;</private-char>m/s$ and the time-step size to $Δt=0.01<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠. The physical constants are chosen as $c=4.412×10−8<private-char specific-use="loader-generated">&#xef07;</private-char>kg/(m·s2)$ and $η=10.0<private-char specific-use="loader-generated">&#xef07;</private-char>kg/(m·s)$⁠, and the kinematic viscosity is $μ=0.001<private-char specific-use="loader-generated">&#xef07;</private-char>m2/s$⁠. In total, 1,000 time steps are performed, corresponding to a final simulation time of $10<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠.

The three-dimensional cylinder flow configuration is constructed by embedding a circular cylinder into a rectangular channel, obtained through a uniform extrusion of the corresponding two-dimensional cylinder geometry in the depth direction (Figure 6). The channel dimensions are $L=2.5<private-char specific-use="loader-generated">&#xef07;</private-char>m$ in length, $H=0.41<private-char specific-use="loader-generated">&#xef07;</private-char>m$ in height and $D=0.41<private-char specific-use="loader-generated">&#xef07;</private-char>m$ in depth. The cylindrical obstacle is centered at $(x,y)=(0.2<private-char specific-use="loader-generated">&#xef07;</private-char>m,<private-char specific-use="loader-generated">&#xef07;</private-char>0.2<private-char specific-use="loader-generated">&#xef07;</private-char>m)$ in the cross-sectional plane, with its axis aligned along the depth (z) direction and extending uniformly across the entire domain. The computational mesh is generated by consistently extending the two-dimensional cylinder mesh into three dimensions, and the domain is discretized using hexahedral finite elements. The entire computational domain is discretized using a total of 11,580 finite elements, resulting in an overall number of 366,702 degrees of freedom.

With regard to the boundary conditions, the same configuration as that used in the corresponding two-dimensional cylinder-in-channel flow problem is adopted. No-slip conditions are imposed on the upper and lower walls, the front and back walls of the channel, as well as on the surface of the cylindrical obstacle, while a traction-free condition is prescribed at the outlet (see Table 4). The inlet boundary condition is specified using the same velocity profile as that used in the three-dimensional channel flow case.

The three-dimensional results for the Newtonian fluid and the two proposed non-Newtonian models are discussed below, focusing on the temporal evolution of the wake structure and the associated internal variable distribution, as evaluated on a cross-section taken at the mid-plane of the three-dimensional channel (⁠$z=0.205<private-char specific-use="loader-generated">&#xef07;</private-char>m$⁠).

For the Newtonian fluid (Figure 17), shortly after the flow becomes globally developed, weak unsteady flow structures associated with shear-layer instabilities are observed in the cylinder wake at time step 200 (corresponding to $t=2<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠). At this stage, the wake does not yet exhibit periodic behavior, and no spatially organized vortex pattern is formed. As the simulation advances to time step 400 (⁠$t=4<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠), wake instabilities become more pronounced, and vortices are intermittently formed and shed, resulting in an aperiodic and transitional wake regime. From time step 600 (⁠$t≥6<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠) onward, vortical structures arranged in an alternating manner are released from the cylinder wake at relatively regular time intervals, indicating the emergence of a sustained vortex-shedding regime. This flow behavior indicates that vortex-shedding mechanisms characteristic of a three-dimensional Kármán-type wake are persistently established in time, while the wake retains its inherently three-dimensional structure.

In contrast to the Newtonian case, the three-dimensional results for the Type 1 model (Figure 18) exhibit a markedly different wake behavior due to the shear-thickening mechanism. This model is characterized by shear-thickening behavior, for which the viscosity function $f(γ)$⁠, defined through the internal variable, gradually increases over time and approaches $f(γ)→1$ in the vicinity of the channel walls and the cylindrical obstacle. Such an increase in viscosity leads to a local rise in the effective viscosity in regions of concentrated shear, with particularly pronounced effects observed near the cylinder surface and the channel walls. As a consequence, local flow acceleration and the rapid development of shear layers are suppressed, and potential shear-layer instabilities in the cylinder wake are significantly mitigated. An examination of the temporal evolution of the velocity field shows that, after an initial transient phase, the flow maintains a relatively stable structure, without the emergence of time-periodic vortex shedding or pronounced unsteady wake instabilities as observed in the Newtonian case. Instead, the wake region downstream of the cylinder exhibits limited amplification of shear-layer instabilities due to the viscosity increase, resulting in a smooth and stable three-dimensional wake structure over time. These results demonstrate that the characteristic shear-thickening-induced flow stabilization of the Type 1 model is consistently preserved in three-dimensional configurations.

For the Type 2 model (Figure 19), which is characterized by shear-thinning behavior, an opposite trend is observed. This model exhibits a viscosity function $f(γ)$ defined through the internal variable, which gradually decreases during the initial stage and approaches $f(γ)→0$ in the vicinity of the cylindrical surface. At this early stage (time step 200), the viscosity reduction is primarily localized around the cylinder, whereas the downstream wake region has not yet experienced a significant decrease in viscosity. As a consequence, up to time step 400 (corresponding to $t≤4<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠), despite the local viscosity decrease near the cylinder, shear-layer instabilities do not sufficiently amplify in the wake, and a comparatively stable flow structure is preserved. As the simulation advances, from around time step 600 (⁠$t=6<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠) onward, regions where $f(γ)$ transitions toward values close to zero begin to form in the downstream wake, and these low-viscosity regions progressively expand in the streamwise direction. Accordingly, the local reduction in effective viscosity gradually spreads throughout the wake region. This viscosity decrease promotes the development of shear layers and the amplification of instabilities in the cylinder wake, leading to a progressive deformation of the vortical structures. At time step 600, increasing wake instability is observed, while for time steps 800 and beyond (corresponding to $t≥8<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠), a clear vortex-shedding behavior emerges in the downstream region. With further time advancement, by time step 1,000 (⁠$t=10<private-char specific-use="loader-generated">&#xef07;</private-char>s$⁠), the flow reaches a vortex-shedding regime characteristic of a three-dimensional Kármán-type wake.

Overall, the three-dimensional cylinder flow results demonstrate that the key characteristics of the non-Newtonian fluid models observed in the two-dimensional cylinder flow are consistently preserved in the three-dimensional setting. Specifically, the Type 1 model exhibits shear-thickening behavior, leading to an increase in the effective viscosity in the vicinity of the cylinder and in the wake region. As a consequence, the development of shear layers and wake instabilities is effectively suppressed, resulting in a comparatively stable wake structure. By contrast, the Type 2 model is characterized by shear-thinning behavior, for which the effective viscosity decreases in the wake region. This viscosity reduction promotes the development of shear layers and the amplification of instabilities, ultimately giving rise to a pronounced vortex-shedding behavior. These trends are qualitatively consistent with the physical mechanisms identified in the two-dimensional cylinder flow and indicate that the proposed internal-variable-based non-Newtonian fluid model is capable of robustly capturing both shear-thickening and shear-thinning responses even with increasing spatial dimensionality.

In this work, we formulated a non-Newtonian fluid model based on an internal variable within the framework of a space–time variational principle and analyzed the resulting flow behavior through numerical simulations. By deriving the weak formulation from Hamilton’s principle, the proposed framework naturally incorporates the temporal and spatial evolution of viscosity. Unlike classical models based on algebraic stress–strain rate relations, the present approach defines viscosity variation through physically motivated energy potentials, enabling a more consistent and thermodynamically sound description.

Two prototype models (Type 1 and Type 2) were applied to benchmark problems of channel flow and cylinder flow. Both settings were considered in two and three spatial dimensions. The results show that the evolution of the internal variable $γ$ leads to fundamentally different flow behavior depending on the model. In the channel flow case, the Type 1 model exhibited sustained viscosity in the center and thickening near the walls, while the Type 2 model, driven by strain energy, resulted in a flattened velocity profile due to anisotropic viscosity distribution. In the cylinder flow scenario, the Type 1 model suppressed Kármán vortex shedding due to localized viscosity gradients, whereas the Type 2 model exhibited strong shear-thinning after a certain time, leading to low-viscosity states and subsequent vortex shedding similar to the Newtonian case.

Furthermore, the internal-variable-based formulation proposed in this study is conceptually distinct from classical empirical non-Newtonian models such as the power-law, Carreau or Bingham models. Traditional approaches prescribe the apparent viscosity through algebraic relations as an instantaneous function of the local shear rate, thereby capturing nonlinear behavior only in a steady-state sense. By contrast, the present framework models viscosity through a spatio-temporal evolution equation for the internal variable $γ$⁠, enabling a dynamic response that naturally incorporates flow history and time-dependent effects while ensuring thermodynamic consistency. This perspective extends the scope of conventional empirical models and provides a physically grounded basis for describing unsteady and history-dependent rheological behavior.

These findings demonstrate that the proposed model effectively captures both spatial and temporal viscosity variations and provides a reliable prediction of complex flow structures. In particular, the ability to represent shear-thinning and shear-thickening behavior within a unified constitutive framework highlights the flexibility and extensibility of the present variational approach for analyzing non-Newtonian fluids.