Stress Waves in Polymeric Fluids

Abstract

This paper demonstrates the existence, propagation, transmission, reflection, and interaction of deviatoric stress waves in polymeric fluids for which the mathematical models are derived using conservation and balance laws (CBL) of Classical Continuum Mechanics (CCM) and the constitutive theories are based on the entropy inequality and representation theorem. The physical mechanisms of deformation in polymeric liquids that enable the stress wave physics are identified and are demonstrated to be valid using Maxwell, Oldroyd-B, and Giesekus polymeric fluids, and are illustrated using model problem studies. We assume polymeric fluids to be isotropic and homogeneous at the macro scale so that the CBL of the CCM can be used to derive their mathematical models. For simplicity, we assume the polymeric fluids to be incompressible in the present work.

Share and Cite:

Surana, K. and Kitchen, M. (2022) Stress Waves in Polymeric Fluids. American Journal of Computational Mathematics, 12, 87-118. doi: 10.4236/ajcm.2022.121007.

1. Introduction, Literature Review and Scope of Work

A polymeric fluid is synthesized using a solvent consisting of short-chain molecules and a polymer consisting of long-chain molecules. The resulting fluid is referred to as a polymeric liquid, or fluid, or viscoelastic fluid. When the composition of a polymeric fluid is dominated by the solvent, it is referred to as a dilute polymeric fluid. On the other hand, if the composition of a polymeric fluid is dominated by the polymer, then it is called a dense polymeric fluid or polymer melt. Early experimental studies aimed at understanding the physics of motion in polymeric fluids using microscopic photography of the flow reveal very complex Brownian motion at the molecular level. It is observed that in the unstressed or relaxed state of a polymeric fluid, the long-chain molecules of the polymer are mostly in a coiled state in the solvent, either by themselves or collectively in a colony of long-chain molecules. These colonies of long-chain molecules are observed to be interconnected with the neighboring colonies.

A solvent and/or a polymer have their own separate viscosities and other properties. Likewise, a polymeric fluid has its own viscosity and other properties that obviously depend on the properties of the solvent and the polymer and their volume fraction, but they must be determined experimentally, just like the properties of a solvent and a polymer.

Upon applying disturbances to a polymeric fluid, the short-chain molecules collectively behave in a similar manner as a Newtonian fluid does. The long-chain molecules, on the other hand, begin to uncoil in the direction of the disturbance (or flow). This requires that they overcome the viscous resistance offered by the solvent as well as the viscous resistance offered by the neighboring polymer molecules. Thus, the motion of the polymer molecules in a polymeric fluid is like a one-dimensional spring. Collectively, the polymer molecules act like one-dimensional springs in the direction of the flow (or disturbance). Just like a one-dimensional spring is only active in its axial direction and has no resistance or response normal to this direction, the polymer molecules behave in a similar manner. That is, normal to the direction of the flow, the polymeric fluid response is quite weak but not zero due to the random orientation and arrangement of the polymer molecules in the polymeric fluid.

Even though at scales lower than the macroscale, polymers are not homogeneous and isotropic (the Brownian motion confirms this), at the macroscale, polymeric fluids can be considered isotropic and homogeneous, i.e., we can use the conservation and balance laws (CBL) of classical continuum mechanics (CCM) and the constitutive theories to derive the mathematical models for deformation physics of polymeric fluids. Surana et al. [1] [2] [3] have presented the derivation of the mathematical models for Maxwell, Oldroyd-B, and Giesekus fluids, both compressible and incompressible using conservation and balance laws of classical continuum mechanics and the derivation of the constitutive theories based on entropy inequality and representation theorem [4] - [20], i.e., theory of isotropic tensors.

There is a vast amount of published work related to the boundary value problems (BVPs) and initial value problems (IVPs) for polymeric fluids. Solutions of IVPs generally report a time response for specific physics of interest, whereas BVPs report solutions that are stationary states of IVPs. To our knowledge, wave propagation physics in polymeric fluids has not been studied and reported. Thus, we do not have any published work to report or discuss in this area. This is in fact the incentive for the work presented in this paper. Generally speaking, we can have two types of waves in continua: 1) waves associated with volumetric deformation; these are pressure waves. Existence of such waves is possible only in compressible continua. 2) The second type of waves is distortional; these are waves of deviatoric stresses. Existence and propagation of these waves require stiffness as well as mass. For example, in incompressible solid continuum that has stiffness and has mass, deviatoric stress waves can exist and propagate, but pressure waves cannot exist as this solid continuum has no volumetric deformation. Similarly, in a compressible gas with virtually no stiffness, the deviatoric stress wave cannot exist, but the pressure waves can exist and propagate as the medium can support volumetric deformation.

Motivation and Scope of Present Work

In the present work we investigate the possibility of existence, propagation, reflection, transmission, and interaction of stress waves in polymeric fluids. First, we consider pressure waves. If the polymeric fluid is incompressible (no volumetric deformation), then pressure waves cannot exist. If the polymeric fluid is compressible, then pressure waves can exist and propagate. Existence of deviatoric stress waves requires stiffness and mass. The viscous resistance offered by the solvent and the neighboring polymer molecules in the uncoiling process of each molecule in the direction of the flow or disturbance can be collectively viewed as one dimensional springs acting in the direction of the disturbance. This is the main source of stiffness in polymeric fluids, and these fluids naturally have mass. Thus, polymeric fluids can support existence and propagation of a deviatoric stress wave related to distortional physics. Since this mechanism of stiffness exists regardless of whether the polymeric fluid is compressible or incompressible, the deviatoric stress wave can exist and propagate in compressible as well as incompressible polymeric fluids. A complete study of wave existence and propagation in polymeric fluids requires the following three investigations.

1) Existence and propagation of distortional stress waves in incompressible polymeric fluids:

a) Existence, propagation, transmission and reflection of distortional stress waves;

b) Factors influencing the speed of propagation of distortional stress waves;

c) Determination of the stress wave speed.

2) Existence and propagation of distortional stress waves in compressible polymeric fluids:

a) Existence, propagation, transmission and reflection of distortional stress wave;

b) Factors influencing the speed of propagation of distortional stress waves;

c) Determination of the stress wave speed.

3) Existence and propagation of pressure waves in compressible polymeric fluids.

Out of the above three studies, study (1) focuses on the fundamental physics, science, and the corresponding mathematical models for deviatoric stress wave existence, propagation, reflection, and transmission in polymeric fluids. The pressure wave physics in compressible continua is well understood and studied. While deviatoric stress waves in compressible polymeric fluids, perhaps in the presence of pressure waves, is certainly of much interest as well, we believe the work proposed in (1) is of more basic importance, hence address this work here. We consider well known incompressible Maxwell, Oldroyd-B, and Giesekus polymeric fluids. The mathematical models (CBL and the constitutive theories) that are strictly based on CCM [1] [2] [3] are considered in the present work.

Mathematical models are non-dimensionalized and their solution for 1D deviatoric stress waves are obtained using space-time coupled finite element method based on the space-time residual functional (STRF) for a space-time strip with time marching [21]. The space-time local approximations in higher order spaces in space and time used in the present work permit higher order global differentiability of space-time approximations for a space-time discretization Ω ¯ x t T of a space-time domain Ω ¯ x t of a space-time strip. Since the space-time differential operator is nonlinear, Newton’s linear method with line search is employed in the STRF formulation with adjustments made to ensure that the resulting space-time integral form from the STRF is space-time variationally consistent [21]. The rate of entropy production due to viscous dissipation is inherent in the physics under consideration. We assume that the model problems domains are completely insulated. Thus, the temperature change (rise) in the deforming volume of matter is only due to entropy generation. For the type of model problems considered in this paper, the temperature rise due to the rate of entropy production is almost negligible ( O ( 10 3 ) or O ( 10 4 ) ) compared to the initial or ambient temperatures of O ( 10 2 ) . Thus, non-isothermal physics has relatively little influence on the existence and propagation of deviatoric stress waves. This permits us to consider isothermal physics only. Thus, the first and the second laws of thermodynamics (FLT, SLT) need not be considered as part of the mathematical model. An outline of the work presented in this paper is summarized in the following.

1) Identification and discussion of the mechanism of stiffness.

2) Mathematical models based on CCM [1] [2] [3] for Maxwell, Oldroyd-B, and Giesekus polymeric fluids, their dimensional form, and their reduction to 1 for 1D axial deviatoric stress wave studies.

3) Details of space-time integral form for a space-time strip and the space-time finite element formulation based on the space-time residual functional (STRF).

4) Model problem studies using Maxwell, Oldroyd-B and Giesekus constitutive models.

a) Wave propagation and reflection.

b) Parametric studies using parameters influencing wave speed.

c) Wave propagation, reflection, transmission, and interaction in domain containing bi-material interfaces.

d) Wave speed determination.

A summary of the work presented in this paper and some conclusions drawn from it are given in the last section of this paper.

2. Mechanisms of Elasticity in Polymeric Fluids

We have explained in the introduction section that drag or resistance offered to the motion of a polymer molecule during uncoiling or stretching due to applied external disturbance is similar to stretching of a one-dimensional spring. Collectively, this physics of polymer molecules in the direction of the disturbance is the source of stiffness in the polymeric fluids. Since the viscous drag forces are proportional to viscosity ( η ) of the polymeric fluids, viscosity η is a parameter that controls presence and the extent of stiffness in polymeric fluids. Higher values of polymeric fluid viscosity η naturally results in more pronounced stiffness. The viscosity of the polymeric fluid naturally depends upon solvent viscosity ( η s ) as well as polymeric viscosity ( η p ). η p is generally much higher than η s . Polymeric fluid viscosity ( η ) is not the sum of η s and η p (even though used so in many published works), but must be determined experimentally using the polymeric fluid. Since η is dependent on η s and η p , both η s and η p influence stiffness in the same manner as η does, i.e., increasing values of η s and η p will result in increasing η that would yield higher stiffness is the polymeric fluid.

The second source of changing elasticity in polymeric fluids is the relaxation time. A smaller relaxation time implies less physical time for the polymer molecules to revert to the relaxed state, hence increased viscous resistance during the relaxation process compared to a larger relaxation time that requires longer physical time for the relaxation process. Thus in this case, the relaxation process is spread over a larger time implying decreased stiffness.

Thus, in polymeric fluids, the speed of propagation of deviatoric stress waves is directly proportional to polymeric fluid viscosity but inversely proportional to the relaxation time. This is obviously due to the fact that the dynamic stiffness of a polymeric fluid is proportional to the viscosity of the polymeric fluid but is inversely proportional to the relaxation time.

In the present work we present model problem studies to determine and demonstrate the influence of viscosity η and relaxation time λ on the speed of the deviatoric stress wave as well as determine the wave speed. We remark that this mechanism of stiffness only exists in polymeric fluids when the polymer molecules are in motion, i.e., the stiffness in polymeric fluids is dynamic. Thus, for a stationary column of polymer fluid the stiffness mechanism is totally absent.

3. Mathematical Model

The CBL of CCM constitutes the core mathematical model for describing the deformation physics of polymeric fluids. These are naturally augmented by the constitutive theories so that the resulting set of partial differential equations has closure. In the early development of the polymer science, the constitutive theories were phenomenological [22] [23], i.e., experimental observation were given mathematical empirical forms and then used in conjunction with CBL so that the mathematical model would have closure. Surana et al. [1] [2] [3] discussed that the polymeric fluid must be assumed isotropic and homogeneous at the macroscale otherwise CBL of CCM may not be valid, and if the polymeric fluid is isotropic and homogeneous, then we must be able to derive the constitutive theories for the polymeric fluids using entropy inequality and the representation theorem [4] - [24]. This is what is in fact presented by Surana et al. [1] [2] [3]. In this paper, we do not present the derivations of the constitutive theories for deviatoric Cauchy stress tensor and, instead, borrow these equations from references [1] [2] [3] and use them in the model problem studies.

Conservation and Balance Laws and the Constitutive Theories

As pointed out earlier in Section 2, in the present work it is valid to assume that the deformation process is isothermal (insulated system and negligible heat generation due to rate of entropy production). Thus, the CBL consists of conservation of mass (CM), balance of linear momenta (BLM) and balance of angular momenta (BAM). The symmetric Cauchy stress tensor (contravariant, see ref. [1] [2] ) σ ¯ ( 0 ) is additively decomposed into equilibrium ( σ ¯ e ( 0 ) ) and deviatoric ( σ ¯ d ( 0 ) ) stress tensors. The constitutive theory for equilibrium stress tensor σ ¯ e ( 0 ) describes volumetric deformation and the distortional deformation is described by the constitutive theory for σ ¯ d ( 0 ) .

σ ¯ ( 0 ) = σ ¯ e ( 0 ) + σ ¯ d ( 0 ) (1)

Since we are considering fluent continua, we must consider Eulerian description for CBL as well as constitutive theories. Conservation and balance laws in Eulerian description for incompressible fluent continua are given by [1] [2] (using notations of reference [1] [2] ).

Conservation and balance laws:

ρ ¯ ( ¯ v ¯ ) = 0 ( CM ) (2)

ρ ¯ D v ¯ D t ρ ¯ F ¯ b ¯ σ ¯ ( 0 ) = 0 ( BLM ) (3)

σ ¯ ( 0 ) = ( σ ¯ ( 0 ) ) T ( BAM ) (4)

Constitutive theories:

We need constitutive theories for σ ¯ e ( 0 ) and σ ¯ d ( 0 ) . Since the deformation is isothermal and the polymeric fluid is assumed incompressible, the constitutive theory for σ ¯ e ( 0 ) is much simplified and is given by [1] [2]

σ ¯ e ( 0 ) = p ¯ δ (5)

in which p ¯ is mechanical pressure, which is assumed positive when compressive. Derivation of the constitutive theory begins with the conjugate pair σ ¯ d ( 0 ) : D ¯ in SLT (see references [1] [2] ) in which D ¯ is the symmetric part of the velocity gradient tensor. Thus, for the incompressible and isothermal physics we can write

σ ¯ d ( 0 ) = σ ¯ d ( 0 ) ( D ¯ ) (6)

If γ ( i ) ; i = 1,2, , n are the convected (covariant) time derivatives of Green’s strain tensor ε [ 0 ] (see [1] [2] ) up to order n, then we find that

D ¯ = γ ( 1 ) (7)

Thus, D ¯ in (6) can be replaced by γ ( 1 ) . Furthermore, if we generalize the argument tensor of σ ¯ d ( 0 ) in (6), then we can replace D ¯ by γ ( i ) ; i = 1,2, , n .

σ ¯ d ( 0 ) = σ ¯ d ( 0 ) ( γ ( i ) ) ; i = 1,2, , n (8)

The argument tensors of σ ¯ d ( 0 ) in (8) suggest that deviatoric Cauchy stress tensor depends upon the convected time derivatives of the Green’s strain tensor up to order n. This is strongly supported by inductive reasoning.

From the published and experimental work in polymer science [1] [2] [22] [23] we know that a constitutive theory for deviatoric Cauchy stress must at least use the first convected time derivative of the deviatoric Cauchy stress as a constitutive variable with deviatoric Cauchy stress tensor as its argument tensor. Otherwise it is not possible to show the existence of the memory modulus, hence a lack of relaxation physics or Rheology. That is, we must at least have (based on (8))

σ ¯ d ( 1 ) = σ ¯ d ( 1 ) ( γ ( i ) , σ ¯ d ( 0 ) ) ; i = 1,2, , n (9)

Let σ ¯ d ( j ) ; j = 1,2, , m be the convected (contravariant) time derivatives of the deviatoric Cauchy stress tensor up to order m, then we can generalize (9) by choosing σ ¯ d ( m ) on a constitutive tensor and σ ¯ d ( j ) ; j = 0,1, , m 1 as its argument tensors in addition to γ ( i ) ; i = 1,2, , n . Thus, we have the following.

σ ¯ d ( m ) = σ ¯ d ( m ) ( γ ( i ) , σ ¯ d ( j ) ) ; i = 1,2, , n ; j = 0,1, , m 1 (10)

Using (10) we could derive a most general constitutive theory for σ ¯ d ( m ) based on integrity and representation theorem. Surana et al. [1] [2] [3] have given the general derivation of Maxwell, Oldroyd-B and Giesekus constitutive models based on integrity. In the present work we consider the constitutive theories for polymeric fluids that are commonly used in published works.

Maxwell Constitutive model:

This is a linear viscoelastic model. In this constitutive theory we consider n = 1 and m = 1 and the resulting linear constitutive theory is given by [1] [2] [3]

σ ¯ d ( 0 ) + λ σ ¯ d ( 1 ) = 2 η γ ( 1 ) (11)

we note that since the fluid is incompressible t r γ ( 1 ) = 0 . λ is relaxation time and η is viscosity of the polymeric fluid. This model is generally used for dilute polymeric fluids.

Oldroyd-B Constitutive model:

Oldroyd-B Constitutive model is referred to as a quasilinear constitutive model. This model can be derived from the general constitutive theory of order n and m by choosing n = 2 and m = 1 and retaining the lowest degree terms [1] [2] [3]

σ ¯ d ( 0 ) + λ 1 σ ¯ d ( 1 ) = 2 η γ ( 1 ) + 2 η λ 2 γ ( 2 ) (12)

in which λ 1 and λ 2 are relaxation and retardation times and η is the viscosity of the polymeric fluid. This model is also generally used for dilute polymeric fluids.

Giesekus Constitutive model:

Giesekus Constitutive model is a nonlinear constitutive model. This constitutive theory can be derived using n = 1 and m = 1 . All terms in this constitutive model are the same as those in Maxwell’s model except that this model contains generator ( σ ¯ d ( 0 ) ) 2 as an additional term [1] [2].

σ ¯ d ( 0 ) + λ σ ¯ d ( 1 ) = 2 η γ ( 1 ) + α λ η ( σ ¯ d ( 0 ) ) 2 (13)

α is called the mobility factor. It accounts for the formation of the polymer colonies connections as well as their separation, i.e., it accounts for mobility of the polymer colonies. This model is generally considered more suitable for dense polymeric fluids or polymer melts.

Unified single constitutive mathematical model for polymeric fluids:

We note all models above are a subset of the general constitutive theory based on n = 2 and m = 1 with appropriate choices of generators and invariants. We could consider

σ ¯ d ( 0 ) + λ σ ¯ d ( 1 ) = 2 η γ ( 1 ) + 2 η λ 2 γ ( 2 ) + α λ η ( σ ¯ d ( 0 ) ) 2 (14)

From (14) we can recover Maxwell, Oldroyd-B and Giesekus constitutive models by appropriate selection of material coefficients and setting the others to zero, i.e.,

Maxwell model:

λ 2 = 0 , α = 0

Oldroyd-B model:

α = 0

Giesekus model:

λ 2 = 0

Complete mathematical model in 3 :

Equations (2)-(4) and (14) constitute the complete mathematical model in BLM (3), constitutive equation σ ¯ e ( 0 ) (1) and σ ¯ d ( 0 ) (6), ten equations in v ¯ (3), p ¯ (1) and σ ¯ d ( 0 ) (6), ten variables. Hence, the mathematical model has closure.

4. Dimensionless Form of the Mathematical Model

We write (2)-(4) and (14) with hat ( ^ ) on each quantity indicating they have their usual units or dimensions.

ρ ¯ ^ ( ¯ ^ v ¯ ^ ) = 0 ρ ¯ ^ D v ¯ ^ D t ^ ρ ¯ ^ F ¯ ^ b + ¯ ^ ( p ¯ ^ δ ) ¯ ^ σ ¯ ^ d ( 0 ) = 0 σ ¯ ^ d ( 0 ) + λ ^ σ ¯ ^ d ( 1 ) = 2 η ^ γ ^ ( 1 ) + 2 η ^ λ ^ 2 γ ^ ( 2 ) + α ^ λ ^ η ^ ( σ ¯ ^ d ( 0 ) ) 2 } (15)

Additive decomposition of σ ¯ ( 0 ) in (1) has been used in BLM. We choose reference quantities (with subscript (0)) and define dimensionless variables. Let L 0 , ρ 0 , v 0 , t 0 , τ 0 , p 0 , η 0 be reference length, density, velocity, time, stress, pressure and viscosity. Then, we define the dimensionless quantity as follows.

x ¯ = x ¯ ^ L 0 , ρ ¯ = ρ ¯ ^ ρ 0 , v ¯ = v ¯ ^ v 0 , = t ^ t 0 , σ ¯ d ( 0 ) = 1 τ 0 ( σ ¯ ^ d ( 0 ) ) , p ¯ = p ¯ ^ p 0 , η = η ^ η 0 wechoose p 0 = τ 0 = ρ 0 v 0 2 : characteristickineticenergy and t 0 = L 0 v 0 } (16)

Using (16) in (15) we can obtain

ρ ¯ ( ¯ v ¯ ) = 0 (17)

ρ ¯ D v ¯ D t ρ ¯ F ¯ b ( p 0 ρ 0 v 0 2 ) ¯ ( p ¯ δ ) ( τ 0 ρ 0 v 0 2 ) ¯ σ ¯ d ( 0 ) = 0 (18)

σ ¯ d ( 0 ) + D e σ ¯ d ( 1 ) = ( η 0 v 0 L 0 τ 0 ) 2 η γ ( 1 ) + ( η 0 v 0 L 0 τ 0 ) 2 η D e 2 γ ( 2 ) + ( τ 0 t 0 η 0 ) α ( D e η ) ( σ ¯ d ( 0 ) ) 2 (19)

using p 0 = τ 0 = ρ 0 v 0 2 and t 0 = L 0 / v 0 , (17)-(19) can be written as

ρ ¯ ( ¯ v ¯ ) = 0 (20)

ρ ¯ D v ¯ D t ρ ¯ F ¯ b ¯ ( p ¯ δ ) ¯ σ ¯ d ( 0 ) = 0 (21)

σ ¯ d ( 0 ) + D e σ ¯ d ( 1 ) = ( 2 η R e ) γ ( 1 ) + 2 η ( D e 2 R e ) γ ( 2 ) + α ( D e R e η ) ( σ ¯ d ( 0 ) ) 2 (22)

in which

D D t = t + v ¯ ¯ (23)

σ ¯ d ( 1 ) = D D t ( σ ¯ d ( 0 ) ) L ¯ σ ¯ d ( 0 ) σ ¯ d ( 0 ) L ¯ T (24)

γ ( 2 ) = D D t ( γ ( 1 ) ) + L ¯ γ ( 1 ) + γ ( 1 ) L ¯ T (25)

[ L ¯ ] = [ { v ¯ } { x ¯ } ] , [ D ¯ ] = 1 2 ( [ L ¯ ] + [ L ¯ ] T ) (26)

In which R e , D e , and D e 2 are Reynolds number and Deborah numbers defined as R e = ρ 0 v 0 L 0 η 0 , D e = λ ^ t 0 , D e 2 = λ ^ 2 t 0 .

5. Dimensionless Mathematical Model in 1

We consider pure axial deformation in x 1 direction. For this case, the continuity equation is identically satisfied, hence does not yield any equations. For 1D in 1 (o-x1 axis of the x-frame) we have the following

Let, L ¯ = v ¯ 1 x ¯ 1 , ¯ = x ¯ 1 , σ ¯ d ( 0 ) = σ ¯ d 11 ( 0 ) , D D t = t + v ¯ 1 x ¯ 1 , σ ¯ d ( 1 ) = σ ¯ d 11 ( 0 ) t + v ¯ 1 σ ¯ d 11 ( 0 ) x ¯ 1 2 v ¯ 1 x ¯ 1 ( σ ¯ d 11 ( 0 ) ) , γ ( 2 ) = 2 v ¯ 1 t x ¯ 1 + v ¯ 1 2 v ¯ 1 x ¯ 1 2 + 2 ( v ¯ 1 x ¯ 1 ) 2 } (27)

using (27) in BLM (21) we can obtain the following in 1 for BLM

ρ ¯ ( v ¯ 1 t + v ¯ 1 v ¯ 1 x ¯ 1 ) ρ ¯ F ¯ 1 b + p ¯ x ¯ 1 x ¯ 1 ( σ ¯ d 11 ( 0 ) ) = 0 (28)

Again, using (27) in (22) we can obtain the following for the constitutive theory in 1

σ ¯ d 11 ( 0 ) + D e ( t ( σ ¯ d 11 ( 0 ) ) v ¯ 1 x ¯ 1 ( σ ¯ d 11 ( 0 ) ) 2 v ¯ 1 x ¯ 1 ( σ ¯ d 11 ( 0 ) ) ) = 2 η R e v ¯ 1 x ¯ 1 + 2 η D e 2 R e ( 2 v ¯ 1 t x ¯ 1 + v ¯ 1 2 v ¯ 1 x ¯ 1 2 + 2 ( v ¯ 1 x ¯ 1 ) 2 ) + α D e R e η ( σ ¯ d 11 ( 0 ) ) 2 (29)

Equation (28) and (29) constitute the dimensionless form of the mathematical model. In this mathematical model v ¯ 1 and σ ¯ d 11 ( 0 ) are dependent variables. p ¯ x ¯ is a specified nonzero value if the flow is pressure driven; otherwise p ¯ x ¯ = 0 . In the case of Maxwell model: D e 2 = 0 and α = 0 , for Oldroyd-B model: α = 0 , and for Giesekus model: D e 2 = 0 .

6. Solution of the Mathematical Model

The mathematical model consisting of (28) and (29) is used to study deviatoric stress waves in Maxwell, Oldroyd-B and Giesekus polymeric fluids. The mathematical model (initial value problem, IVP) consists of two nonlinear partial differential equations in dependent variables v ¯ 1 and σ ¯ d 11 ( 0 ) and independent variables x ¯ 1 and time t. We consider numerical solutions of (28) and (29) using space-time coupled finite element method in which space-time integral form is derived using space-time residual functional (STRF). Newton’s linear method with line search [21] is used to obtain solutions of the non-linear algebraic equations resulting from STRF method.

When studying evolution described by IVPs, it is beneficial to consider discretization Ω ¯ x t T of space-time domain Ω ¯ x t into space-time strips (Figure 1(a)), i.e., Ω ¯ n x t T = e Ω ¯ x t e in which Ω ¯ i x t is the ith space-time strip. The solution is calculated for a space-time strip and then time marched beginning with the first space-time

Figure 1. Discretization of space-time domain into space-time strips, discretization of nth space-time strip into space-time elements. (a) A space-time domain Ω ¯ x t T and its discretization in space-time strips; (b) Discretization of nth space-time strip Ω ¯ n x t in space-time element.

strip [21]. Figure 1(b) shows discretization Ω ¯ x t T = i Ω ¯ i x t of an nth space-time strip Ω ¯ n x t using nine-node p-version hierarchical space-time elements with higher order global differentiability local approximation in space and time. In the following we present details of the STRF finite-element formulation and the solution procedure for nonlinear algebraic equations based on Newton’s linear method with line search [21]. Consider an nth space-time strip Ω ¯ n x t .

Let ( v ¯ 1 ) h e and ( σ ¯ d 11 ( 0 ) ) h e be local approximations of v ¯ 1 and σ ¯ d 11 ( 0 ) over a space-time element Ω ¯ x t e , and let ( v ¯ 1 ) h and ( σ ¯ d 11 ( 0 ) ) h be approximations of v ¯ 1 and σ ¯ d 11 ( 0 ) over the discretization Ω ¯ n x t T of the nth space-time strip Ω ¯ n x t . Then,

( v ¯ 1 ) h = e ( v ¯ 1 ) h e and ( σ ¯ d 11 ( 0 ) ) h = e ( σ ¯ d 11 ( 0 ) ) h e (30)

Upon substituting ( v ¯ 1 ) h and ( σ ¯ d 11 ( 0 ) ) h in (28) and (29), we obtain residual functions E 1 and E 2 (we assume absence of body forces)

E 1 = ρ ¯ ( ( v ¯ 1 ) h ) t + p ¯ h x ¯ 1 x ¯ 1 ( ( σ ¯ d 11 ( 0 ) ) h ) = 0 E 2 = ( σ ¯ d 11 ( 0 ) ) h + D e ( t ( ( σ ¯ d 11 ( 0 ) ) h ) ( ( v ¯ 1 ) h ) x ¯ 1 ( ( σ ¯ d 11 ( 0 ) ) h ) 2 ( ( v ¯ 1 ) h ) x ¯ 1 ( ( σ ¯ d 11 ( 0 ) ) h ) ) 2 η R e ( ( v ¯ 1 ) h ) x ¯ 1 2 η D e 2 R e ( 2 ( ( v ¯ 1 ) h ) t x ¯ 1 + ( ( v ¯ 1 ) h ) 2 ( ( v ¯ 1 ) h ) x ¯ 1 2 + 2 ( ( ( v ¯ 1 ) h ) x ¯ 1 ) 2 ) α D e R e η ( ( σ ¯ d 11 ( 0 ) ) h ) 2 = 0 } (31)

We define the residual functional I using E 1 and E 2 over Ω ¯ n x t T

I = I 1 + I 2 = ( E 1 , E 1 ) + ( E 2 , E 2 ) (32)

We assume that I is differentiable in its arguments (i.e. ( v ¯ 1 ) h and ( σ ¯ d 11 ( 0 ) ) h ), then δ I is unique and δ I = 0 is a necessary condition for an extremum of the residual functional I.

δ I = 2 ( E 1 , δ E 1 ) + 2 ( E 2 , δ E 2 ) = { g 1 } + { g 2 } = { g } = 0 (33)

and the extremum principle or the sufficient condition is given by δ 2 I [21].

δ 2 I = 2 ( δ E 1 , δ E 1 ) + 2 ( δ E 2 , δ E 2 ) + 2 ( E 1 , δ 2 E 1 ) + 2 ( E 2 , δ 2 E 2 ) (34)

while ( δ E 1 , δ E 1 ) > 0 , ( δ E 2 , δ E 2 ) > 0 δ E 1 and δ E 2 , this is not true for 2 ( E 1 , δ 2 E 1 ) and 2 ( E 2 , δ 2 E 2 ) terms. Thus, (34) at this stage does not give unique extremum principle or sufficient condition. Let

( v ¯ 1 ) h e = [ N v ] { δ v e } ( σ ¯ d 11 ( 0 ) ) h e = [ N σ ] { δ σ e } } ( v ¯ 1 ) h e and ( σ ¯ d 11 ( 0 ) ) h e V H k , p ( Ω ¯ x t e ) , k 2 (35)

be the local approximation of v ¯ 1 and σ ¯ d 11 ( 0 ) over Ω ¯ x t e in which { δ v e } and { δ σ e } are nodal degrees of freedom. V is the approximation space containing functions of class C 1 or higher. Let { δ v } and { δ σ } be the degrees of freedom for v ¯ 1 and σ ¯ d 11 ( 0 ) for the discretization Ω ¯ n x t T of the nth space-time strip, then { δ v } = e { δ v e } ; { δ σ } = e { δ σ e } . We can write the following using (31) and (32)

I = e I 1 e + e I 2 e = e ( E 1 e , E 1 e ) Ω ¯ x t e + e ( E 2 e , E 2 e ) Ω ¯ x t e (36)

in which E 1 e and E 2 e for Ω ¯ x t e can be obtained from (32) by replacing ( v ¯ 1 ) h and ( σ ¯ d 11 ( 0 ) ) h by ( v ¯ 1 ) h e and ( σ ¯ d 11 ( 0 ) ) h e . Then,

δ I = e 2 ( E 1 e , δ E 1 e ) + e 2 ( E 2 e , δ E 2 e ) = e { g 1 e } + e { g 2 e } = { g ( { δ } ) } = 0 (37)

in which { δ } T = [ { δ v } T , { δ σ } T ] are the nodal degrees of freedom for Ω ¯ n x t T . Since the mathematical model consists of nonlinear partial differential equations, { g ( { δ } ) } in (37) is a nonlinear function of { δ } . We consider Newton’s linear method with line search to find a solution { δ } that satisfies δ I = { g } = 0 (in (37)). Let { δ } 0 be an assumed (known) solution, then

{ g ( { δ } 0 ) } 0 (38)

Let { Δ δ } be a correction to { δ } 0 such that

{ δ } = { δ } 0 + { Δ δ } (39)

which satisfies (37), i.e.

{ g ( { δ } ) } = { g ( { δ } 0 + { Δ δ } ) } = 0 (40)

Expanding { g ( { δ } ) } in Taylor series about { δ } 0 and retaining only up to linear terms in { Δ δ } .

{ g ( { δ } ) } = { g ( { δ } 0 + { Δ δ } ) } { g ( { δ } 0 ) } + [ { g } { δ } ] { δ } 0 { Δ δ } = 0 (41)

{ Δ δ } = [ { g } { δ } ] { δ } 0 1 { g ( { δ } 0 ) } (42)

We note that

[ { g } { δ } ] { δ } 0 = [ { δ I } { δ } ] { δ } 0 = [ δ 2 I ] { δ } 0 (43)

Surana [21] [24] have shown that δ 2 I in (34) can be approximated by

δ 2 I 2 ( δ E 1 , δ E 1 ) + 2 ( δ E 2 , δ E 2 ) > 0 δ E 1 and δ E 2 (44)

Thus, now we have an extremum principle. Approximation (44) has no influence on the STRF that ends with δ I = 0 . It merely alters the slope of the tangent plane to the hypersurface δ I = 0 at { δ } 0 during Newton’s linear method. Rather than using (39) to obtain the updated solution, we choose

{ δ } = { δ } 0 + α { Δ δ } , 0 < α < 2, suchthat I ( { δ } ) I ( { δ } 0 ) (45)

This is referred to as line search. This helps in accelerating the convergence of Newton’s Linear method. It is well known that in all iterative methods one must choose a starting or initial solution. For Newton’s linear method to work well, this choice must be in a close neighborhood of the correct solution. In all numerical studies presented in this paper, { δ } 0 = 0 (null) has been used for all degrees of freedom that are not part of the boundary conditions and initial conditions. This choice has worked well in all numerical studies.

7. Summary of Main Steps in Newton’s Linear Method with Line Search

1) choose { δ } 0 .

2) calculate { g ( { δ } 0 ) } If | g i | Δ , a preset tolerance of zero, then { δ } 0 is the solution, skip remaining steps, otherwise follow remaining steps.

3) calculate [ δ 2 I ] using (44).

4) calculate { Δ δ } using (42).

5) updated solution { δ } is given by { δ } = { δ } 0 + α { Δ δ } .

6) set { δ } = { δ } 0 and go to step 2.

8. Model Problems

In this section we present studies to demonstrate existence of deviatoric stress waves and study their propagation, reflection, transmission, and interaction using model problems in 1 . We consider incompressible Maxwell, Oldroyd-B, and Giesekus viscoelastic fluids. Assumptions used in deriving CBL and the constitutive theories remain valid in these studies as well.

Figure 2(a) show a polymeric fluid column of uniform cross-section and of length L. Left end ( x 1 = 0 ) is considered impermeable and the right end ( x 1 = L ) is subjected to a wave of σ ¯ d 11 ( 0 ) of duration 2 Δ t (Figure 2(d)). If we assume that all points on a given cross-section move in the x 1 direction by the same amount, then we can idealize the column of Figure 2(a) by a line representation shown in Figure 2(b). Figure 2(c) shows a space-time finite element discretization Ω ¯ 1 x t T of the first space-time strip Ω ¯ 1 x t using nine-node p-version hierarchical space-time finite elements [21] [24]. Boundary conditions (BCs) and initial conditions (ICs) are also shown in Figure 2(c). A compressive σ ¯ d 11 ( 0 ) is applied at x 1 = L (shown in Figure 2(d)).

σ ¯ d 11 ( 0 ) = 0 ; σ ¯ d 11 ( 0 ) t = 0 at t = 0.0 σ ¯ d 11 ( 0 ) = σ ¯ d 11 ( 0 ) ; σ ¯ d 11 ( 0 ) t = 0 at t = Δ t σ ¯ d 11 ( 0 ) = 0 ; σ ¯ d 11 ( 0 ) t = 0 at t 2 Δ t } (46)

Figure 2. 1D fluid domain, idealization of 1D fluid domain, discretization of first space-time strip with space-time finite elements and applied disturbance. (a) Polymeric fluid domain; (b) Mathematical idealization of fluid domain; (c) Disretization of first space-time strip using 10 space-time p-version finite elements; (d) Deviatoric stress ( σ ¯ d 11 ( 0 ) ) pulse of duration 2Δt.

We consider the following material coefficients for Maxwell, Oldroyd-B and Giesekus fluids.

Maxwell model:

ρ ¯ ^ = 868 kg / m 3 , λ ^ 1 = 0.1 s , η ^ = 3.0 Pa s , ρ 0 = 868 kg / m 3 , L 0 = 0.02 m , v 0 = 0.05 m / s , t 0 = L 0 / v 0 = 0.4 s , η 0 = 3.0 Pa s , p 0 = τ 0 = ρ 0 v 0 2 = 2.17 Pa (47)

Oldroyd-B model:

ρ ¯ ^ = 868 kg / m 3 , λ ^ 1 = 0.1 s , η ^ = 3.0 Pa s , λ ^ 2 = 0.001 s , ρ 0 = 868 kg / m 3 , L 0 = 0.02 m , v 0 = 0.05 m / s , t 0 = L 0 / v 0 = 0.4 s , η 0 = 3.0 Pa s , p 0 = τ 0 = ρ 0 v 0 2 = 2.17 Pa (48)

Giesekus model:

ρ ¯ ^ = 800 kg / m 3 , λ ^ = 0.06 s , η ^ = 1.426 Pa s , α ^ = 0.15 , ρ 0 = 800 kg / m 3 , L 0 = 0.02 m , v 0 = 0.05 m / s , t 0 = L 0 / v 0 = 0.4 s , η 0 = 1.426 Pa s , p 0 = τ 0 = ρ 0 v 0 2 = 2.0 Pa (49)

In all numerical studies we consider Ω ¯ x t T = e Ω ¯ x t e , a uniform discretization of Ω ¯ x t = [ 0,1 ] × [ 0, Δ t ] domain for the first and the subsequent space-time strips of same spatial length using space-time p-version hierarchical elements with higher order global differentiability in space and time. Since the mathematical models are a system of first order PDEs in space and time, the approximation space V H k , p ( Ω ¯ x t ) with k = 2 in space and time (equal order) in minimally conforming approximation space. We choose k = 1 in space and time, i.e., local approximations of class C 1 in space and time. This ensures that all space-time integrals over Ω ¯ n x t T of the nth space-time strip are in Riemann sense. Equal degree approximation are considered, thus p ξ = p η = p . We consider p 3 . In all numerical studies, p-convergence studies are conducted for the first space-time strip to determine p-levels for which the integrated sum of squares of the residual functional I is O ( 10 8 ) or lower to ensure that the computed solution is converged. These p-levels suffice subsequent space-time strips as well in yielding similar values of residual functional I. This ensures time accuracy of the entire evolution for model problem solutions presented in the following.

8.1. Existence, Propagation, and Reflection of Deviatoric Stress Waves

A schematic of the compressive deviatoric stress wave of peak value 0.01 applied over 0 t 2 Δ t at x = L = 1.0 as shown in Figure 1(d). We consider Δ t = 0.01 . Evolutions are computed for Maxwell, Oldroyd-B and Giesekus viscoelastic fluent continua. For all three fluids, a p-level of 11 is found to produce sufficient I values of the order of O ( 10 8 ) or lower for each space-time strip.

8.1.1. Maxwell Fluid

We perform computations of evolution using Maxwell fluid properties show in (47). Figure 3 and Figure 4 show plots of σ ¯ d 11 ( 0 ) versus distance x for 0 t 41 Δ t . At t = Δ t and t = 2 Δ t , a half and full wave are observed in the spatial domain confirming existence of the deviatoric stress wave. As time elapses, the wave propagates to the left. At t = 3 Δ t , the full wave has propagated to the left with amplitude decay and some (minor) base elongation. At t = 19 Δ t , the stress wave reaches the impermeable boundary at x = 0 . We observe measurable amplitude decay for the wave at t = 19 Δ t . Reflection from the impermeable boundary begins

Figure 3. Time evolution of deviatoric stress ( σ ¯ d 11 ( 0 ) ) in Maxwell fluid.

Figure 4. Continued time evolution of deviatoric stress ( σ ¯ d 11 ( 0 ) ) in Maxwell fluid.

at t = 20 Δ t and is completed at t = 21 Δ t . The reflected wave at t = 21 Δ t propagates to the right ( t = 22 Δ t , t = 23 Δ t ) and reaches the free boundary at the right at t = 38 Δ t with substantially reduced amplitude but without measurable base elongation. The reflection process from the free boundary occurs between t = 39 Δ t and t = 40 Δ t . At t = 40 Δ t , we have a reflected tensile wave that propagates to the left ( t = 41 Δ t ).

The study confirms existence, propagation, and reflection of deviatoric stress waves in the Maxwell fluid. Amplitude decay of the stress wave due to viscous dissipation is clearly observed as the evolution proceeds in time. We do not observe measurable base elongation of the deviatoric stress wave during propagation.

8.1.2. Oldroyd-B Fluid

Studies similar to the Maxwell fluid are performed for Oldroyd-B fluid using the fluid properties listed in (48). The same discretization and p-levels are used for solutions of C 1 as is the case of Maxwell fluid. Figure 5 and Figure 6 show

Figure 5. Time evolution of deviatoric stress ( σ ¯ d 11 ( 0 ) ) in Oldroyd-B fluid.

Figure 6. Continued time evolution of deviatoric stress ( σ ¯ d 11 ( 0 ) ) in Oldroyd-B fluid.

plots of σ ¯ d 11 ( 0 ) versus x for various values of t ( 0 t 48 Δ t ). At t = Δ t , we observe a half wave in the spatial domain with peak value same as the peak value of the applied wave. At t = 2 Δ t , we observe a full wave in the spatial domain but with significantly reduced peak due to dissipation. At t = 3 Δ t , the full wave propagates further towards the left. At t = 14 Δ t , the stress wave is at the impermeable boundary at x = 0 . We observe significant amplitude decay and base elongation. At t = 20 Δ t , reflection occurs. The reflected wave with significantly reduced peak and elongated base is shown at t = 21 Δ t . Further propagation of the reflected wave to the right, reflection from the free boundary, and the reflected tensile wave and its propagation to the left for 31 Δ t t 48 Δ t is shown in Figure 6. Due to significant amplitude decay and base elongation for 20 Δ t t 48 Δ t , scales for σ ¯ d 11 ( 0 ) are changed for better resolution.

The study confirms existence, propagation, and reflection of deviatoric stress waves in the Oldroyd-B fluid. The Oldroyd-B fluid constitutive theory contains second convected time derivative of Green’s strain tensor, γ ( 2 ) , (not present in Maxwell model) responsible for excessive base elongation.

8.1.3. Giesekus Fluid

We consider Giesekus fluid with properties given in (49). We use the same discretization as in Maxwell and Oldroyd-B fluids with the same p-levels to compute evolution for 0 t 44 Δ t . We discuss results in the following.

Figure 7 and Figure 8 show plots of σ ¯ d 11 ( 0 ) versus x for different values of t ( 0 t 44 Δ t ). The results are very similar to Maxwell model presented in Section 8.1.1. At t = Δ t and t = 2 Δ t , we observe half and full stress waves in the spatial domain. The wave propagates to the left with amplitude decay but without significant base elongation. At t = 20 Δ t , the wave reaches the impermeable boundary (at x = 0 ). Wave reflection occurs for t = 21 Δ t , 22 Δ t , and 23 Δ t . The reflected wave propagates to the right ( t = 23 Δ t , t = 24 Δ t ) and reaches the free boundary (at x = 1.0 ) at t = 41 Δ t . The reflected tensile wave ( t = 43 Δ t ) propagates to the left.

This study also confirms existence, propagation, and reflection of the deviatoric stress wave in the Giesekus fluid. The wave evolution shows significant amplitude decay but almost non-visible base elongation. This behavior is similar to Maxwell fluid.

8.2. Deviatoric Stress Wave Speeds

As discussed earlier in polymeric fluids the dynamic stiffness is a function of viscosity and relaxation time. Higher viscosity and lower relaxation time result in higher stiffness, hence higher speed of propagation of deviatoric stress waves. The purpose of the studies in this section is to determine and report wave speeds for Maxwell, Oldroyd-B and Giesekus fluids as a function of viscosity and relaxation time.

For a chosen value of viscosity and relaxation time, we determine evolution for an applied deviatoric stress wave (as in Section 8.1). For a given time t 1 , the peak of the wave is located in spatial domain x 1 . Evolution is continued and the peak of the wave is relocated at time t 2 > t 1 , then the wave speed is simply

S d = Δ x 1 t 2 t 1 (50)

Figure 7. Time evolution of deviatoric stress ( σ ¯ d 11 ( 0 ) ) in Giesekus fluid.

where Δ x 1 is the distance traveled by the wave during the time t 2 t 1 . Since for a fixed η and λ the wave speed is constant, this procedure of calculating wave speed works quite accurate.

Figures 9(a)-(c) show graphs of wave speed S d as a function of viscosity η

Figure 8. Continued time evolution of deviatoric stress ( σ ¯ d 11 ( 0 ) ) in Giesekus fluid.

for three different values of relaxation time λ for Maxwell, Oldroyd-B and Giesekus fluids. Since Maxwell and Oldroyd-B constitutive models are for dilute polymeric fluids, the same values of λ are used in Figure 9(a) and Figure 9(b) so that wave speeds can be compared between the two.

Figure 9. Wave speed of σ ¯ d 11 ( 0 ) for Maxwell, Oldroyd-B and Giesekus fluids as a function of viscosity and relaxation time. (a) Deviatoric stress wave speed for Maxwell fluid; (b) Deviatoric stress wave speed for Oldroyd-B fluid; (c) Deviatoric stress wave speed for Giesekus fluid.

From the results presented in Figure 9 we clearly observe increasing wave speed for increasing viscosity and for decreasing relaxation time in all three fluids. Wave speeds in case of Oldroyd-B fluids are slightly higher than Maxwell fluid for the same values of η and λ . In Giesekus fluid the wave speeds are much higher than in the case of Maxwell or Oldroyd-B fluids, this of course is due to much higher viscosity of the Giesekus fluid compared to Maxwell and Oldroyd-B fluids.

8.3. Deviatoric Stress Wave Propagation, Reflection, Transmission, and Interaction in the Presence of Bimaterial Interface

In this study we consider a dimensionless spatial domain of two units with bimaterial interface located at the middle of the domain. We consider Maxwell fluid to present two studies. We consider a twenty space-time element uniform discretization in both studies.

Case (a): In the first study, η 1 > η 2 , where η 1 and η 2 are viscosities of the fluids to the right (material M 1 ) and the fluid to the left (material M 2 ) of the bimaterial interface. We choose η 1 = 9 and η 2 = 3 .

Case (b): In the second study we consider η 1 = 3 for M 1 and η 2 = 9 for M 2 .

In both studies, a compressive deviatoric stress wave of duration 2 Δ t with peak value of 0.01 is applied at the free boundary ( x 1 = 2.0 ). The left boundary at x 1 = 0.0 is assumed impermeable. Evolution in both studies is calculated using Δ t = 0.01 with p ξ = p η = 11 . In both cases for each time strip, I of O ( 10 8 ) or lower is achieved indicating good accuracy of the computed evolution.

8.3.1. Discussion of Results: Case (a)

We discuss results for case (a) shown in Figure 10 and Figure 11 first. At t = 2 Δ t , the deviatoric stress wave is completely in the spatial domain. At t = 12 Δ t , the wave is at the bimaterial interface located at x = 1.0 . At t = 13 Δ t , a part of the compressive stress wave transmits and a part of it reflects as a tensile wave. Both are moving at different speeds to the left and to the right of the bimaterial interface. At t = 22 Δ t , the tensile wave to the right of the interface reaches the free boundary first as it is moving faster than the transmitted compressive wave in material M 2 . At t = 30 Δ t , the compressive wave in material M 2 reaches the impermeable boundary at x 1 = 0.0 and reflects as a compressive wave that propagates to the right. The reflected compressive wave from the free boundary continues to move towards the bimaterial interface. At t = 33 Δ t , the wave in material M 1 reaches the interface and goes through transmission and reflection ( t = 35 Δ t ). The reflected tensile wave in M 1 propagates towards the free boundary, the transmitted compressive wave interacts with the wave in material M 2 ( t = 41 Δ t ). At t = 42 Δ t , the two interacting wave are coincident, but assume their own identity at t = 43 Δ t and continue propagating in the same direction

Figure 10. Wave interaction between two Maxwell fluids of viscosities η1 and η2 where η1 > η2.

( t = 44 Δ t ) as they were before they became coincident. In the meantime, the wave in M 1 continues to propagate towards the free end.

Figure 11. Continued wave interaction between two Maxwell fluids of viscosities η1 and η2 where η1 > η2.

8.3.2. Discussion of Results: Case (b)

Details of the evolution from this study are shown in Figure 12 and Figure 13. Since η 1 < η 2 , the wave propagates slower in medium M 1 compared to medium

Figure 12. Wave interaction between two Maxwell fluids of viscosities η1 and η2 where η1 < η2.

M 2 . This completely changes the propagation, transmission, reflection and interaction details of the waves during evolution. At t = 2 Δ t , the wave enters the spatial domain and propagates to the left toward the bimaterial interface located at x 1 = 1.0 . Transmission and reflection of the incident wave at x 1 = 1.0 is seen at t = 21 Δ t . At t = 22 Δ t , we see the transmitted wave propagating toward the

Figure 13. Wave interaction between two Maxwell fluids of viscosities η1 and η2 where η1 < η2.

impermeable boundary located at x 1 = 0.0 and the reflected wave propagates towards the free boundary located at x 1 = 2.0 . At t = 30 Δ t , the transmitted wave reaches the impermeable boundary at x 1 = 0.0 and the reflected wave continues to propagate towards the free boundary. At t = 37 Δ t , the reflected wave from the boundary at x 1 = 0.0 is propagating toward the interface and the wave in M 1 approaches the free boundary located at x 1 = 2.0 . At t = 43 Δ t , a part of the wave in M 2 is transmitted in M 1 as a compressive wave and the other part is reflected back in M 2 as a tensile wave at the interface, while the tensile wave in M 1 continues to propagate toward the interface. At t = 49 Δ t , the tensile and compressive waves are about to coincide resulting in a compressive wave with smaller peak value than the original compressive wave ( t = 50 Δ t ). At t = 51 Δ t , the two coincident waves (at t = 50 Δ t ) separate into the original waves and propagate as shown for t = 52 Δ t .

Both Studies in case (a) and case (b) show existence, propagation, transmission, reflection and interaction of the deviatoric stress waves in the polymeric incompressible fluids. Viscosities η 1 and η 2 are suitably chosen to demonstrate the influence of two different wave speed in the wave physics interaction during evolution.

9. Summary and Conclusions

In this paper we have demonstrated and presented existence, propagation, reflection, transmission, and interaction of deviatoric stress waves in incompressible polymeric fluids using Maxwell, Oldroyd-B, and Giesekus polymeric fluids. The constitutive models [1] [2] [3] for these fluids are based on CCM and have been derived using entropy inequality and representation theorem [4] - [20]. In the following we present a summary of the work and draw some conclusions.

1) Polymeric fluid viscosity and relaxation time are identified as the main sources of “dynamic stiffness” in polymeric fluids. In a stationary column of polymeric fluid, the stiffness is absent.

2) Resistance offered due to viscous drag to the stretching or uncoiling of long chain polymer molecules in the direction of the flow (or disturbance) is similar to the extension of a one dimensional spring and is one of the sources of dynamic stiffness in polymeric fluids. Since this process only exists when the polymer molecules are in motion, the stiffness is dynamic. Polymeric fluids with higher viscosity have higher dynamic stiffness. That is, the dynamic stiffness of a polymeric fluid is proportional to polymeric fluid viscosity.

3) The second source of dynamic stiffness in polymeric fluids is the relaxtion time. A polymeric fluid with higher relaxation time takes a longer time to come back to the relaxed (unstressed) state upon cessation of disturbance, implying reduced resistive forces during the relaxation process, hence reducing dynamic stiffness. On the other hand, a polymeric fluid with smaller relaxation time resumes the unstressed or relaxed state quicker upon cessation of disturbance. This process naturally corresponds to higher resistive forces during relaxation, hence increased dynamic stiffness. Thus, dynamic stiffness of a polymeric fluid is inversely proportional to the relaxation time.

4) Since the polymer fluid mass is kept fixed or constant, the increased dynamic stiffness corresponds to faster speed of propagation of deviatoric stress waves. Likewise, reduced dynamic stiffness results in slower speed of propagation of deviatoric stress waves in polymeric fluids.

5) Based on (2)-(4), we conclude the higher viscosity of the polymeric fluid and lower relaxation time result in higher speeds of propagation of deviatoric stress waves. Likewise, lower viscosity of polymeric fluid and higher relaxation time results in lower speeds of deviatoric stress waves in polymeric fluids.

6) Model problem studies are presented using Maxwell, Oldroyd-B and Giesekus polymeric fluids to show existence, propagation, reflection, transmission and interaction of deviatoric stress waves in a single polymeric fluid as well as a medium containing two different polymeric fluids with bimaterial interface. In all cases, converged solutions reported in the paper are free of any oscillations and the residual functional (I) values of the order of O (108) or lower are achieved during the entire evolution for each space-time strip confirming good accuracy of the reported solutions.

7) Results presented for the deviatoric stress waves speed as a function of polymeric fluid viscosity for different values of the relaxation time confirm that wave speed is proportional to polymeric fluid viscosity but inversely proportional to the relaxation time.

8) The work presented in this paper is fundamental science, and, to our knowledge, constitute a new discovery related to the dynamics of the polymeric fluid physics.

9) Accuracy of the deviatoric stress wave solution, i.e., the base elongation and the amplitude decay during the propagation and the accuracy of the speed of the wave propagation are dependent on the accuracy of the material properties of the polymeric fluid. The finite element results reported in this paper are virtually dispersion free, and the residual functional (I) values are O (108) or lower ensured accurate solutions of the partial differential equations constituting the mathematical model. The choice of the minimally conforming space (k = 2) ensures that the space-time integrals over the space-time discretization of a space-time strips are Riemann, thus confirming accuracy of computation of the residual functional I, the overall measure of solution accuracy.

Acknowledgements

The first author is grateful for his endowed professorships and the Department of Mechanical Engineering of the University of Kansas for providing financial support to the second author. The computational facilities provided by the Computational Mechanics Laboratory of the Mechanical Engineering Department are also acknowledged.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1] Surana, K.S. (2015) Advanced Mechanics of Continua. CRC/Taylor and Francis, Boca Raton, FL.
[2] Surana, K.S. (2021) Classical Continuum Mechanics, Second Edition. CRC/Taylor and Francis, Boca Raton, FL.
[3] Surana, K.S., Nunez, D., Reddy, J.N. and Romkes, A. (2014) Rate Constitutive Theory for Ordered Thermoviscoelastic Fluids: Polymers. Continuum Mechanics and Thermodynamics, 26, 143-181.
https://doi.org/10.1007/s00161-013-0295-8
[4] Boehler, J.P. (1977) On Irreducible Representations for Isotropic Scalar Functions. Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 57, 323-327.
https://doi.org/10.1002/zamm.19770570608
[5] Prager, W. (1945) Strain Hardening under Combined Stresses. Journal of Applied Physics, 16, 837-840.
[6] Reiner, M. (1945) A Mathematical Theory of Dilatancy. American Journal of Mathematics, 67, 350-362.
https://doi.org/10.2307/2371950
[7] Rivlin, R.S. (1955) Further Remarks on the Stress-Deformation Relations for Isotropic Materials. Journal of Rational Mechanics and Analysis, 4, 681-702.
https://doi.org/10.1512/iumj.1955.4.54025
[8] Rivlin, R.S. and Ericksen, J.L. (1955) Stress-Deformation Relations for Isotropic Materials. Journal of Rational Mechanics and Analysis, 4, 323-425.
https://doi.org/10.1512/iumj.1955.4.54011
[9] Smith, G.F. (1970) On a Fundamental Error in Two Papers of C.C. Wang, “On Representations for Isotropic Functions, Part I and Part II”. Archive for Rational Mechanics and Analysis, 36, 161-165.
https://doi.org/10.1007/BF00272240
[10] Smith, G.F. (1971) On Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 9, 899-916.
https://doi.org/10.1016/0020-7225(71)90023-1
[11] Spencer, A.J.M. (1971) Part III. Theory of Invariants. In: Eringen, A.C., Ed., Mathematics, Academic Press, Cambridge, 239-353.
https://doi.org/10.1016/B978-0-12-240801-4.50008-X
[12] Spencer, A.J.M. and Rivlin, R.S. (1959) The Theory of Matrix Polynomials and Its Application to the Mechanics of Isotropic Continua. Archive for Rational Mechanics and Analysis, 2, 309-336.
https://doi.org/10.1007/BF00277933
[13] Spencer, A.J.M. and Rivlin, R.S. (1960) Further Results in the Theory of Matrix Polynomials. Archive for Rational Mechanics and Analysis, 4, 214-230.
https://doi.org/10.1007/BF00281388
[14] Todd, J.A. (1948) Ternary Quadratic Types. Philosophical Transactions of the Royal Society of London. Series A: Mathematical and Physical Sciences, 241, 399-456.
https://doi.org/10.1098/rsta.1948.0025
[15] Wang, C.C. (1969) On Representations for Isotropic Functions, Part I. Archive for Rational Mechanics and Analysis, 33, 249-267.
https://doi.org/10.1007/BF00281278
[16] Wang, C.C. (1969) On Representations for Isotropic Functions, Part II. Archive for Rational Mechanics and Analysis, 33, 268-287.
https://doi.org/10.1007/BF00281279
[17] Wang, C.C. (1970) A New Representation Theorem for Isotropic Functions, Part I and Part II. Archive for Rational Mechanics and Analysis, 36, 166-223.
[18] Wang, C.C. (1971) Corrigendum to ‘Representations for Isotropic Functions’. Archive for Rational Mechanics and Analysis, 43, 392-395.
https://doi.org/10.1007/BF00252004
[19] Zheng, Q.S. (1993) On the Representations for Isotropic Vector-Valued, Symmetric Tensor-Valued and Skew-Symmetric Tensor-Valued Functions. International Journal of Engineering Science, 31, 1013-1024.
https://doi.org/10.1016/0020-7225(93)90109-8
[20] Zheng, Q.S. (1993) On Transversely Isotropic, Orthotropic and Relatively Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors, and Vectors. International Journal of Engineering Science, 31, 1399-1409.
https://doi.org/10.1016/0020-7225(93)90005-F
[21] Surana, K.S. and Reddy, J.N. (2017) The Finite Element Method for Initial Value Problems: Mathematics and Computations. CRC/Taylor and Francis, Boca Raton.
https://doi.org/10.1201/b22512
[22] Bird, R.B., Armstrong, R.C. and Hassager, O. (1987) Dynamics of Polymeric Liquids, Volume 1, Fluid Mechanics. 2nd Edition, John Wiley and Sons, Hoboken.
[23] Bird, R.B., Armstrong, R.C. and Hassager, O. (1987) Dynamics of Polymeric Liquids, Volume 2, Kinetic Theory. 2nd Edition, John Wiley and Sons, Hoboken.
[24] Surana, K.S. and Reddy, J.N. (2016) The Finite Element Method for Boundary Value Problems: Mathematics and Computations. CRC/Taylor and Francis, Boca Raton.
https://doi.org/10.1201/9781315365718

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.