Viscoelastic flows studied by smoothed particle dynamics

J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

Viscoelastic flows studied by smoothed particle dynamics Marco Ellero a , Martin Kröger a,b,∗ , Siegfried Hess a a

Institut für Theoretische Physik, Technische Universität Berlin, D-10623 Berlin, Germany b Polymer Physics, Institute of Polymers, ETH Zürich, CH-8092 Zürich, Switzerland Received 8 November 2001; received in revised form 2 April 2002

Abstract A viscoelastic numerical scheme based on smoothed particle dynamics is presented. The concept goes a step beyond smoothed particle hydrodynamics (SPH) which is a grid-free Lagrangian method describing the flow by fluid-pseudo-particles. The relevant properties are interpolated directly on the resulting movable grid. In this work, the effect of viscoelasticity is incorporated into the ordinary conservation laws by a differential constitutive equation supplied for the stress tensor. In order to give confidence in the methodology we explicitly consider the non-stationary simple corotational Maxwell model in a channel geometry. Without further developments the scheme is applicable to ‘realistic’ models relevant for three-dimensional (3D) viscoelastic flows in complex geometries. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Smoothed particle hydrodynamics; Viscoelasticity; Isothermal flow; Fluid-pseudo-particles; Jaumann–Maxwell model

1. Introduction Until recently, numerical simulations of viscoelastic flows in complex geometries have often been based on a purely macroscopic approach where conservation laws together with a specific rheological constitutive equation [1] are taken into account. This approach was in some cases sufficient for managing, controlling and optimizing the manufacturing processes. Although progress in macroscopic viscoelastic flow computations has been very impressive [2,3], the development of numerical schemes for time-dependent and three-dimensional (3D) flows still receives considerable attention. In this article, we extend the grid-free smoothed particle hydrodynamics (SPH) simulation method to smoothed particle dynamics. The novel approach—still abbreviated as SPH in the following—denotes a macroscopic, dissipative, formulation which allows one to describe viscoelastic isothermal flows through the use of a model differential constitutive relation for the stress tensor. Ordinary SPH discretized equations for the conservation laws are then supplemented by this model. Before introducing the extended SPH ∗

Corresponding author. E-mail addresses: [email protected] (M. Ellero), [email protected] (M. Kröger). 0377-0257/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 0 2 5 7 ( 0 2 ) 0 0 0 5 9 - 9

36

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

method we summarize recent developments in closely related approaches to the discretization of complex fluid flows. A ‘micro–macro’ approach to viscoelastic flow simulations emerged in [4] that combines the solution of the conservation laws with the direct use of a kinetic theory model describing the fluid’s rheology (e.g. [5–8]). The approach is computationally expensive and combines standard finite element methods with stochastic simulation. In this case, there is no need for a rheological constitutive equation to describe the fluid. In [9] a numerical technique referred to as the Lagrangian particle method (LPM) had been proposed. LPM is used for solving time-dependent viscoelastic flows on a grid with either the macroscopic or the micro–macro approach. It combines in a decoupled fashion, the Eulerian solution of the conservation laws using Galerkin’s finite element technique with a Lagrangian computation of the extra-stress. Evaluations are done at a number of discrete particles that are convected by the flow. The extra-stress is computed by integrating either the relevant differential constitutive equation (macroscopic approach), or the stochastic differential equation associated with the kinetic theory model (micro–macro approach). In the micro–macro LPM simulations, each Lagrangian particle convected by the flow carries an ensemble of particles with internal degrees of freedom. For the case where these degrees are statistically correlated, the numerical approach is equivalent—in the limit of vanishing discretization error—to the method of Brownian configuration fields [7,10]. Relatively new is the attempt to attack flow problems in complex geometries with SPH, a Lagrangian ‘macro’ method that was developed 25 years ago for astrophysical problems by Lucy [11], and Gingold and Monaghan [12]. SPH is a fully Lagrangian scheme permitting the interpolation of the flow properties directly at a discrete set of pseudo-particles, i.e. without the need to define a spatial mesh. Its Lagrangian nature, associated with the absence of a fixed grid, is its main strength. It removes difficulties concerning the convective terms and allows one to tackle fluid and solid flow problems involving large deformations and free surfaces in a relatively natural way. In the last 10 years, many SPH simulations have been applied to different physical situations including compressible Newtonian flows [13], incompressible free surface flows [14], high strain mechanics, ultrarelativistic shocks [15], impact and sliding friction problems between elasto-plastic materials [16–19], numerical fluid [20] and gas dynamics in astrophysics [21]. The SPH method has also been applied to problems in kinetic theory, such as the dynamics of homogeneous liquid crystals [22] based on a Fokker–Planck approach. Implementations of standard SPH for parallel architectures recently became available [23]. The method has received substantial theoretical support in order to make it consistent from a statistical point of view, i.e. correct treatments of thermal fluctuations and consistent fluctuating hydrodynamics [24]. In [25], a close connection is made between SPH and dissipative particle dynamics (DPD), a fully Lagrangian method in which the force between the particles are modeled on the mesoscale, following the philosophy of lattice Boltzmann or lattice gas automata schemes but with the advantage of the flexibility of being gridless. Recently, application of DPD to viscoelastic flow has been proposed consisting in assigning to every DPD particle one or more internal structural variables [26]. There are also other ways of discretizing the Navier–Stokes equations for which boundary conditions are very simple to apply, whereas for SPH (and DPD) these are problems which need to be addressed. One such possibility is through finite volumes, which is essentially at the same level as SPH, but includes thermal fluctuations depending on the physical size of the fluid particles [27]. In this paper, we extend the grid-free ‘macro’ SPH formulation to a macroscopic, dissipative, formulation which allows one to describe viscoelastic isothermal flows through the use of a model differential constitutive relation for the anisotropic (symmetric traceless) part of the stress tensor. Ordinary SPH

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

37

discretized equations for the conservation laws are then supplemented by this model. For purely illustrative purposes we choose the nonlinear Jaumann–Maxwell (i.e. corotational Maxwell) model, which in turn is obtained by replacing the partial derivative appearing in the linear Maxwell model with a Jaumann derivative which contains a nonlinear advective term plus a term which takes into account solid-body rotations. Numerical results are discussed for unsteady relaxation processes of an initial sinusoidal velocity profile. An analytic reference solution for the linearized regime is readily available in that particular nonequilibrium problem.

2. Basic SPH formalism The SPH has been motivated in the context of interpolation theory. In order to illustrate the underlying idea, we consider an arbitrary well behaved function, f defined over a domain of interest. The function may represent some physical variable or a density. It can be expressed in terms of its values at a discrete set of disordered points (SPH particle positions) by a suitable definition of an interpolation kernel. This is illustrated next. We start out with some basic definitions. Given a function f defined over an entire domain, we introduce an integral estimate f  at point x as follows
f (x) = f (x  )W (x − x  , h) dx  , (1) where W (x − x  , h) is an interpolating kernel. By definition the kernel has the following properties
W (x − x  , h) d x  = 1, (2) and lim W (x − x  , h) = δ(x − x  ).

h→0

(3)

These features assure proper normalization and consistency in the continuum limit. Here, h represents the range of the interpolation function if it is supposed to have a compact support, or, for example in the case of Gaussian kernel, its width at half height. It was shown in [28] that we can expand f (x  ) in a Taylor series around x since W (x − x  , h) is a strongly peaked function at x = x  . If W (x, h) is an even function in its first label, in particular for the case W (x, h) = W (|x|, h), the term of order O(h) vanishes automatically. Then Eq. (1) becomes f (x) = f (x) + C∇ 2 f (x)h2 + O(h3 ),

(4)

where the coefficient C is insensitive to h. Provided that f varies smoothly on a length scale of order h, the relation gives a leading error proportional to h2 and results in an accuracy for the SPH discretization which is second-order in space. An additional approximation is required concerning the estimate of the integral in (1) as a sum over points of the domain. We, therefore, introduce an instantaneous number density n(x), with  n(x) ≡ δ(x − x j ). (5) j

38

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

We can then associate the infinitesimal integration volume dx  with the quantity φj (volume occupied by the j th particle) by using the replacement rule mj dx  → φj ≡ . (6) ρj Here, we introduced the concepts of mass mj and density ρj = mj n(x j ) associated with the j th particle. The final SPH approximation for a function f , at position x, reads  mj f (x) fj W (|x − x j |, h), (7) ρ j j with fj ≡ f (x j ). An exhaustive review on the mathematical foundation of SPH is provided, e.g. by Benz [28] or Monaghan [29,30]. 3. Momentum equation The derivation of the SPH equations for the conservation of mass and momentum is obtained following the approach of Benz [28]. It actually consists in an integration over all the domain of the partial differential equations describing the conservation laws of continuum mechanics. An integration by parts removes a surface term permitting direct differentiation of the analytic interpolation kernel. Particularly, in the case of momentum conservation, the equation which we discretize is dV α 1 (8) = ∇ α σ αβ . dt ρ We use tensorial notation with Greek indices indicating spatial coordinates and the summation convention for repeated indices is used. The velocity field is denoted by V α and σ αβ is the general stress tensor. The resulting SPH momentum equation reads [28,31]  αβ  αβ  σj dViα σi α = mj + 2 ∇iα Wij + aext,i , (9) 2 dt ρ ρ i j j α represents the where roman indices indicate particles, greek indices indicate components, and aext,i acceleration acting on the ith particle due to the existence of boundaries. In Eq. (9), we abbreviated Wij ≡ W (|x i − x j |, h). The property ∇i Wij = −∇j Wij implies the validity of Newton’s third law and also the total momentum of bulk particles is exactly conserved. This is a property of this particular set of SPH equations, though one should keep in mind that there is no unique derivation of the set. An infinite number of different sets all describing the same physics to the second-order accuracy can in principle be formulated. However, Eq. (9) is special in the respect that it conserves momentum. Alternative forms of SPH have been discussed by Monaghan [29].

4. Density evaluation By substituting the mass density f (x) = ρ(x) into (7), we obtain the following expression:  ρ(x) mj W (|x − x j |, h). j

(10)

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

39

This expression offers a direct way to evaluate the density explicitly as a sum over the particles and justifies the original denomination SPH. Particles are not points but their masses are smoothly smeared over a region with radius h. It should be noted that relation (10) does not involve a computational complexity of order N 2 , where N is the total number of particles. As long as the support of the kernel function is compact, a finite number of particles contributes to a quantity evaluated at an arbitrarily chosen space coordinate. In this case many acceleration techniques can be used, such as ‘linked cell’ algorithms [32] being linear in the number of particles N. An advantage of using expression (10) competing with other formulas is that the mass is exactly conserved as long as the number of particles is constant. On the other hand, it poses a number of problems concerning boundary and edge effects. Artificial boundary layers can cause instabilities and/or slow down the code. In order to produce reliable trajectory evaluations near the walls the time step chosen must be very small. We will reconsider this problems in detail in the section on boundary conditions. In this work we follow the approach described by Monaghan [29]. Thus, we solve the continuity equation for the density dρ/dt = −ρ∇ α V α rather than evaluating a sum over particles to extract the current density. As a consequence, the density becomes an additional dynamical variable. In analogy to the procedure already described for the momentum equation, for the density we obtain the following SPH discretization  mj dρi = −ρi (Viα − Vjα )∇iα Wij . (11) dt ρ j j This strategy solves the problems associated with boundaries and also establishes a computational advantage, i.e. all the variable fields can be evaluated in one step without the need of a pre-loop for the density calculation. On the other hand, the relation (11) does not retain the exact mass conservation. 5. Constitutive equation The set of conservation equations for density and momentum is not closed, we therefore, need a closure relationship which expresses the stress tensor as function of known variables. In the case of an isothermal flow we require σ αβ = σ αβ (ρ, V α ). The symmetric stress tensor is identically rewritten as the sum of an isotropic plus a symmetric traceless (anisotropic) tensor σ αβ = −pδ αβ + σ αβ ,

(12)

where p is the ordinary hydrostatic pressure responsible of volumetric changes and σ αβ is an anisotropic tensor related to shape deformations preserving volumes. The anisotropic part also covers the normal stress effects. 5.1. Bulk deformation: equation of state A closure relation for p in (12) is given from the equation of state relating hydrostatic pressure to local density. Many possible forms for this can be adopted depending on the particular problem. In the case of an incompressible fluid, a realistic equation of state is too expensive to be solved computationally.

40

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

We are going to choose the following one:  γ  ρ −1 , p(ρ) = p0 ρ0

(13)

where the reference density ρ0 (always taken as ρ0 = 1 in this work) relates to a vanishing hydrostatic reference pressure p(ρ0 ) = 0. Actually, SPH is a method developed mainly for compressible flow problems, but this artificial equation proposed initially by Batchelor [33] with γ = 7 for the description of water can reduce considerably the effects due to compressibility. As discussed by Monaghan [29] and Morris et al. [34], it is always possible to choose a value of the parameter p0 such that the resulting sound speed c with   ∂p(ρ) p0 7 , (14) c(ρ) ≡ ∂ρ ρ0 is suitably larger than a typical flow velocity V characterizing the model. It has been argued by Monaghan that a choice of c 10 times larger than this typical velocity V diminishes density fluctuations under a level comparable with the squared Mach number M 2 (M ≡ V /c), which in this case corresponds to relative density fluctuations of about 1%. This condition of quasi incompressibility introduces errors which are consistent with the global accuracy of the SPH interpolation process. An exact treatment of the incompressibility should take into account a further kinematic constraint (e.g. vorticity-stream function SPH formulation) on the velocity field which ensures that it is divergence-free [35,36]. No such modifications of the basic SPH scheme are used in our work. Compressibility effects are reduced to a small amount by using the state equation (13). 5.2. Shear and extensional flows: Maxwell model The goal in this section is to apply the procedure to a simple model subjected to shear and extensional deformations. Our choice is a nonlinear corotational Maxwell model for the components of the anisotropic (symmetric and traceless) stress tensor σ [37]. This model is chosen for purely illustrative purpose, and has been used in a limited number of applications [18]. The symmetric strain rate tensor is defined as follows:   1 ∂V β ∂V α ˙ αβ = . (15) + 2 ∂x α ∂x β This tensor characterizes the local total deformation (shape and volume) of the material. Changes in volume are controlled, in our scheme, uniquely from the equation of state, while the constitutive Maxwell equation involves only the symmetric traceless (anisotropic) strain rate tensor which is the traceless part of the symmetric strain rate tensor (15), ˙ αβ ≡ ˙ αβ −

1 αβ γ γ 1 ∂V γ αβ δ , δ ˙ = ˙ αβ − d d ∂x γ

(16)

where d is the dimensionality of the problem. The corotational Maxwell model [37] relates the anisotropic stress to anisotropic strain rate d αβ G σ + σ αγ ωγβ + σ γβ ωαγ = 2G ˙ αβ − σ αβ , dt η

(17)

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

41

where G is the (high frequency) shear modulus, η the viscosity of the specimen and the left hand side of (17) is the Jaumann time derivative of σ describing the time evolution of a fixed volume in a Lagrangian way subjected to rotations. The vorticity tensor ωαβ is defined as ωαβ ≡

(∂V β /∂x α ) − (∂V α /∂x β ) . 2

Finally, the SPH model equations for the anisotropic part of the stress tensor are obtained, from (17), in the same way as those for momentum or mass conservation giving   G αβ 2 d αβ αγ γβ αβ . (18) σ = σi (kγβ − kβγ ) + σi (kαγ − kγ α ) − σi +2G kαβ + kβα − kγ γ δ η d dt i The SPH discretization of the transposed velocity gradient appearing through strain rate (15) and (16) and rotation rate (vorticity) tensors, explicitly appearing in (18), reads 1 ∂V µ 1  mj µ µ ∂Wij (V − V ) ν. (19) j i 2 ∂x ν 2 j ρj ∂x Note that the anisotropic part of the stress tensor has normal stress components, and one-third of the trace of the total stress tensor equals the negative hydrostatic pressure, cf. Eq. (12). There are comparable approaches using differential constitutive equations adapted to simulate elasto-plastic materials [16–19]. In these works, the term due to elasticity in (17) remains the same but a further relation, e.g. based on the von Mises criterion is used to recognize the plastic regime where the stress stays constant.

6. Dissipation It has often been observed that numerical solutions usually present large unphysical oscillations and are unstable if no dissipative term is added to the equations. This is because shocks are present mostly in the first stages when initial configurations relax. Strong instabilities can occur if the shock fronts are not smeared out well enough on a length scale larger than the discretization step h. For the moment we restrict our discussion to simple Newtonian flows and to the relative momentum equations. The natural way to proceed should be to discretize, in the SPH formalism, the contribution of the viscous stress in the momentum equations and to add it to the hydrostatic pressure term. This approach has been followed by some authors but, in order to evaluate this term, nested sums over the particles, or second-order differentiation of the kernel function are involved in the corresponding algorithm. As pointed out by Monaghan [29], this produces a very high computational effort in the first case or low accuracy in the second one. The interpolation accuracy (using second-order derivative of the kernel) is strongly sensitive to the particle disorder if the number of neighbors is not very large. The most common remedy is the one following Gingold and Monaghan [12]. They introduce a completely artificial viscosity term (we keep the notation ‘viscosity’ for convenience, even if it does not have the dimension of a viscosity) in the momentum equation acting only when shocks are present and smearing out discontinuities, mimicking in some way a physical mechanism. The artificial viscosity is actually build up by two contributions: one linearly proportional to the divergence of velocity and another one substantially analogous to the Von Neumann–Richtmeyer viscosity [38] able to handle high Mach

42

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

number flows. This artificial viscosity, Πij which creates a force on particle i in the presence of particle j as explicitly used in Eq. (24), reads  2  (−α c¯ij µij + βµij ) , Vijα xijα < 0, Πij = (20) ρ ¯ ij  α α 0, Vij xij ≥ 0, where ρ¯ij ≡ 21 (ρi + ρj ),

c¯ij ≡ 21 (ci + cj ),

Vijα ≡ Viα − Vjα ,

xijα ≡ xiα − xjα ,

(21)

and µij ≡

hVαij xijα xijα xijα + 0.01 h2

,

(22)

were used. These expressions introduce dissipative effects able to reproduce quite accurately many viscous flow situations. However, there is no analytical calculation which would prove that this motion corresponds to the one prescribed by the Navier–Stokes equation and consequently, the resulting measured viscosity of the specimen is not exactly known a priori, but has to be measured. In Section 7, we implement this form of an artificial viscosity only near the wall in order to enforce a no-slip boundary condition. We do not need to take it into account for the bulk fluid. Because we consider low Mach number flows we have set β = 0 and we found that a value of α of order 10 or slightly above is needed to obtain a vanishing tangential velocity near the surface. Smaller values of this parameter produce a measurable slip wall velocity. In our simulations we choose α = 10. The term 0.01h2 is introduced in order to prevent singularities when two particles become too close. Although numerical results have been obtained which are in good agreement with experimental observations for complex flows [14], it has also been noted that for simple test cases, where analytical solutions and/or previous exact calculations are available, this term produces inaccuracies in the velocity profiles. These are actually due to an excessive artificial shear viscosity giving a too large vorticity decay and unphysical momentum transfer. Several ways to escape this difficulty have been proposed. One possibility is to take into account the vorticity associated with every particle, thus, assuring the artificial viscosity vanishes in pure shear flows [34]. Alternative approaches are based on the definition of a further viscosity parameter associated with a source-decay equation causing increase and decay of this parameter in an entering and outgoing shock, respectively [39]. Our approach is based on the solution of a further dynamic equation for the stress tensor based on the Maxwell model. Dissipation does not enter into the momentum equation but rather through the model equation via the corresponding time-dependent stress tensor (and the related time-dependent viscosity). To this end let us consider the dimensionless Deborah number, De = τ/T , which gives an estimate of the elastic effects in the simulated flow by comparing the typical relaxation time τ = η/G for the Maxwell model with a characteristic time T = L/V where the inertial effects are dominant. The quantities L and V denote, respectively, a typical length and velocity characterizing the problem. Here, L is the linear box dimensions of the cubic simulation cell. In the limiting case De → 0 the elastic behavior is negligible and the evolution equation (17) corresponds to the usual constitutive stress-strain relation for Newtonian viscous flow σ αβ = 2η ˙ αβ = 2η(kαβ +kβα )−(4/d)ηkγ γ δ αβ , where d is the dimensionality of the system. In the same limit, we have a continuum expression for the stress tensor which corresponds exactly to

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

43

that one used by Takeda et al. [13] in their simulations. However, the stress tensor discretized in the momentum equation by, e.g. Takeda et al. involved a second-order differentiation of the kernel function while, in the present work, two independent first-order derivations are applied, respectively, to the stress tensor in the momentum Eq. (9), and to the velocity field in (16). The viscosity enters implicitly in the momentum equation via the stress tensor. The way we treat dissipation quite naturally suggests avoiding an artificial viscosity. It retains dissipation in a form consistent with the Navier–Stokes equations. Neither nested sums over particles nor second derivatives of the kernel function are involved in the time loop. Compared with classical SPH, there is a weak increase of memory costs by use of further independent variables, i.e. the stress components. 7. Boundary conditions Boundary conditions receive special attention in SPH. A number of strategies have been followed during the last years. One may consider fixed boundary particles at the wall surfaces interacting with the SPH bulk particles via an artificial repulsive force which prevents penetration. Problems arise if the particle density has to be evaluated near the wall. As mentioned in Section 6, when we consider just one layer of equispaced boundary particles the density of particles decreases rapidly from a bulk value to approximatively half the bulk value upon approaching the surface. This is because for a particle just on the surface, only the single hemisphere intersecting the inner domain is taken into account in the smoothing process. Via Eq. (13), the resulting density gradient produces a pressure gradient in the momentum equation. Consequently, particles accelerate towards the wall producing an inhomogeneity in the dynamical grid and artificial layers parallel to the surface. An obviously improved strategy is to fill the exterior, i.e. at least at a range of depth comparable with the support h of the weight function, with boundary particles which reproduce the desired internal averaged density. Then one also considers the exterior particles for the evaluation of averaged quantities. This approach, due to Morris et al. [34], is known as “image boundary particle”. A third and more convincing approach is to evaluate the contribution of a gas of ‘ghost particles’ existing outside the bulk domain. It is convincing since one exactly recovers SPH in the continuum limit. No artificial terms are introduced, and there is no external particles. Rather complicated surface integrals must be solved when complex geometries come into play, cf. [13]. In this work, yet another approach suggested by Monaghan is adopted. We consider a layer of fixed boundary particles at the wall surface, not being involved in the averaging process, see Fig. 1. The boundary particles interact with the SPH bulk particles via a two-body repulsive force vector of Lennard–Jones (LJ) type, [(r0 /r)4 − (r0 /r)2 ]r −2 rα , r ≤ r0 , α fwall (r) ∝ (23) 0, r > r0 , This force (actually we consider the correcponding acceleration but keep the symbol) prevents particles from crossing the border via adjustment of a proportionality factor in (23), whose precise value has little influence on the overall result. The value for r0 should depend on the particle density via the mean distance between particles and should satisfy r0 < h/2. Explicitly, in this work we have chosen h = 0.08 and r0 = 0.03. Just for particles approaching the wall, i.e. within a distance less than the cutoff radius for Lucy’s function, we take into account a viscous force based on the artificial viscosity expression. This procedure yields boundary conditions suitable for most fluid flow problems, where the parameters α, β,

44

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

Fig. 1. SPH bulk particles (white) and fixed boundary particles (black). A total number of 992 particles is used in the SPH simulation. Periodic boundary conditions have been imposed in the x-direction, the velocity gradient lies in the orthogonal direction, hence, normal to the rigid walls.

and r0 determine the slip-velocity, which in turn is affected by the flow parameters. Not only the artificial viscosity, but also the LJ interactions couple the velocities of bulk particles to the ones of the boundary particles, but the latter contribution is negligible. The complete surface contribution (active within a layer close to the walls) to the acceleration appearing in (9), therefore, reads   α α aext,i mj Πij ∇iα Wij + fwall,ij , (24) =− j α α with Πij from (20), Wij ≡ W (|x i − x j |, h), and fwall,ij ≡ fwall (|x i − x j |). The term appearing in the sum represents the effect of boundary particle j to the bulk particle i. We solve the physical conservation equations in the bulk domain. Within a thin layer (determined by r0 ) artificial terms are included in the exact discretized SPH bulk equations to capture wall effects.

8. Kernel The next step is to choose a kernel function satisfying the requirements. In the simulation of viscous compressible flows, Takeda et al. used a simple Gaussian function which offers stability properties and is infinitely differentiable, but does not possess compact support. The missing compact support, however, prevents one from using fast methods linear in the number of particles. We choose a kernel which combines

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

45

high computational efficiency with good smoothing properties. The most frequently used kernel is the cubic spline kernel  3  r 2 3  r 3 r   1 − + , 0 ≤ ≤ 1,   2 h 4 h h      r 3 r W (r, h) = w0 1 2− , 1 < ≤ 2,   4 h h      0, otherwise, where w0 is a normalization factor which in two dimensions (2D), for example, takes the value 10/(7π h2 ). The advantage of this particular kernel is, i.e. has compact support, is spherically symmetric, the leading order term is of order h2 . The first and second derivatives are continuous which implies that the interpolation is not too sensitive to particle disorder [29]. The normalization factor may be slightly adjusted in the simulation since the continuous weight function is evaluated at a finite number of discrete points. For the purpose of clarity of the results, here we do not include such an adjustment. This will produce small deviations between analytic and numerical solution in Section 10. Other kernels have been proposed which interpolate more accurately in space (fourth-order), as for example the quintic spline kernel or the super-Gaussian kernel. These kernels possess disadvantages, such as expensive computational performance, and they eventually take locally negative values. In the present work, we choose for the kernel the Lucy weighting function originally introduced by Lucy and Monaghan in the first SPH papers [12] and more recently used in problems dealing with ideal gases [40,41], i.e. we employ  r 3 r    1+3 1− , r ≤ h, h h W (r, h) = w0 (25)   0, r > h, where w0 = 5/(πh2 ) in 2D. It interpolates at second-order in h, has a compact support and its first and second derivatives are continuous. 9. Time integration scheme We adopt a simple two-step predictor–corrector (P–C) scheme with second-order accuracy in time to calculate the evolution of the independent variables. It should be replaced if higher precision is essential. In the following, the superscript indicates the discrete time step. The P–C scheme consists of an Eulerian explicit first evaluation of all quantities (predictor step): x˜in+1 = xin + Vin  t/2, V˜in+1 = Vin + Fin  t/2, ρ˜in+1 = ρin + Gni  t/2, and σ˜ ijn+1 = σijn + Hijn  t/2, where F , G, H denote the sources, cf. Eqs. (9), (11) and (18), evaluated as follows Gn = G(ρ n , Vin ), Fin = Fi (pn , ρ n , σijn ), Hijn = Hij (ρ n , Vin , σijn ). Once we have obtained the predicted values of all the independent variables (for every particle) we evaluate the source over these states in order to calculate the final corrected variables: xin+1 = xin + (Vin + V˜in+1 )  t/2 and analogous equations for Vin+1 , ρin+1 , σijn+1 . This step involves quantities, such as ˜ n+1 ≡ G(ρ˜ n+1 , V˜in+1 ). G

46

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

The source terms appearing in the integration scheme are obtained by summations over the SPH particles. The term Fi is explictly given on the right hand side of Eq. (9). For particles closer to boundaries than a distance r0 , the external contribution (24) is active. It involves the LJ force (23) plus the contribution due to the artificial pressure Π (20) –(22). We accommodate a Courant–Friedrichs–Lewy [42] stability condition to make sure that there is no numerical propagation of signals faster than the speed of sound c. Accordingly, for the time step we choose t ≤ h/c. In addition, we consider the strength of forces acting to limit the time step. In order to ensure a safe integration we ensure that t is considerably smaller than the inverse Einstein frequency of the system.

10. Numerical results for the prototype In order to verify the viscoelastic properties of this SPH fluid, which contains the corotational Maxwell model as a prototype, we consider a case for which a very simple analytical solution exists. The problem will be related to an equation of telegraphers type, a detailed discussion can be found in the textbook of Joseph [43]. We consider a 2D channel and impose periodic boundary conditions in the x-direction. An extension to 3D flows does not provide any complication, all equations given are valid in 3D too. A particularly simple time-dependent solution is represented by the viscoelastic relaxation of an initial

Fig. 2. The x-component of the velocity vs. its y-component at four different times, where the flow parameters are De = 1/4, Re = 1. The initial profile has sinusoidal shape. Nine hundred and ninety-two particles are contained in the simulation (unit length L = 1). The cutoff radius for Lucy’s function is set to h = 0.08. This implies an average of 20 neighbors per particle. An appropriate time step chosen is t = 0.0001. The analytical solution is for the linearized model and serves as a reference. Differences are discussed in Sections 8 and 10, as well as in Fig. 4.

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

47

sinusoidal 2D velocity profile V x (y, t = 0) = sin (πy). A peculiarity of a viscoelastic flow is its ability to show oscillatory behavior during the relaxation process. This behavior is intimately connected with the fluids transient response when subjected to unsteady shear flows. For simple Newtonian flows, the initial profile decays monotonically to the stationary rest state. If elasticity is considered, oscillations around this final profile appear. In order to compare the numerical solution with an analytical one, we consider the linearized version of the Jaumann–Maxell equation (17) for which the corotational derivative is replaced by a partial one. This corresponds to linearizing around the stationary state and neglecting advective and rotational terms. After suitable dimensionalization, the equation for the x-component of the velocity field reads ∂2 x 1 ∂ x 1 ∂2 x V (y, t) + V (y, t) = V (y, t), ∂t 2 De ∂t De Re ∂y 2

(26)

where De and Re are Deborah and Reynolds numbers, respectively, De ≡ τ/T = (η/G)/(L/V ), and Re = VρL/η. In these expressions, L and V denote, respectively, the typical width of the channel and the initial velocity in its center. In our simple test case the gradient of the velocity in the x-direction vanishes and the two normal stress coefficients are considered negligible so we can eliminate the advective term in the momentum equation. A particular solution has the following form V x (y, t) = A(y) exp(−βt) sin (ωt + φ),

(27)

Fig. 3. Relaxation of the x component of the velocity at the center of the channel vs. time. Flow parameters are De = 1/10, Re = 5/2, otherwise same conditions as for Fig. 2. The SPH values for the velocity were sampled from those particles whose y-coordinates were within a range of half the maximum interaction distance at the center of the channel. Time and velocity are both expressed in dimensionless units.

48

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

where 1 , β= 2 De

 ω=

π2 1 , − De Re 4 De2

(28)

and A solves the harmonic equation ∂ 2 A/∂t 2 = −π 2 A [43]. The particular initial condition we assume corresponds to A(y) = sin (πy) and φ = π/2. The solution (28) exhibits different behaviors depending of the combination of the two dimensionless numbers De, Re. In the case De → 0, the frequency becomes imaginary and the relaxation process is characterized by a strong exponential overdamped decay, and therefore, similar to ideal Newtonian case. In the opposite regime, if the ratio De/Re is not so small, elastic effects become comparable with the inertial ones and a characteristic underdamped decay exhibiting oscillations in the velocity profile around the stationary rest state is observed. Next we present numerical results for two typical decay processes. First, we consider a viscoelastic flow with De = 1/4, Re = 1. Values for the SPH particle velocities in the center of the channel are presented in Fig. 2 and also compared with the analytical solution, Eqs. (27) and (28). The frequency in (28) is real valued and quite strong oscillations due to elastic effects are observed. The velocity profile still remains sinusoidal but its amplitude becomes sometimes negative. Irrelevant deviations between analytical and

Fig. 4. Relaxation of the x-component of the velocity at the center of the channel for three cases. The rotational terms of the chosen differential model (17) for the stress tensor are included (SPH 1) and disregarded (SPH 2) in the SPH simulations. The latter one (line) is to be compared with the analytical solution of the linearized model (triangles). Remaining small differences are kept which can be reduced by adjusting the normalization factor for the weight function. The flow corresponds in all cases to De = 1/4, Re = 1, cf. Fig. 2. The inset provides a zoom into the region where the expected differences between the three solutions are significant.

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

49

numerical solution already at t = 0 are due to the fact that we do not account for a small correction concerning proper normalization of the continuous weight function with regard to the finite number of neighbors. Second, results of a simulation with De = 1/10, Re = 5/2 are presented in Fig. 3. The ratio Re/De is roughly six times larger than in the first case and effects due to elasticity are less evident but still present. The x component of the velocity becomes still slightly negative. This is so, because the exponential decay factor β in Eqs. (27) and (28) is inversely proportional to De, the characteristic time for elastic effects to appear is very small. The frequency is still positive, and therefore, oscillations exist. Given dimensionless numbers whose ratio is further reduced, the frequency becomes imaginary and a form different from (27) must be considered. In that case, one obtains an overdamped decay where the velocity profile relaxes monotonically. In order to quantify the approximation made by linearizing the Jaumann–Maxwell model, and to test the validity of the foregoing analysis we present an additional result for De = 1/4, Re = 1 which enables us to identify the origin of deviations between analytical solution and simulation result. The reference run named SPH 1 corresponds to the general constitutive equation based on the Jaumann–Maxwell model. For another run SPH 2 the rotational terms, i.e. the second and third term on the left hand side of Eq. (17) are omitted. The latter one is expected to mostly reproduce the given analytical solution of the linearized model, while the former one represents the numerical result of the full model. Small deviations between the analytical solution and the corresponding numerical solution SPH 2 can be assigned to the finite ‘diameter’ of fluid particles (therefore, irrelevant) and to the advective term which is still present in the numerical solution and which is inherent to the Lagrangian nature of SPH. The transverse component of the velocity in the center of the channel versus time is plotted in Fig. 4 for all three cases. As expected, the result SPH 2 is closer to the analytical solution than SPH 1. In Figs. 2 and 3, the differences between analytic and numerical solutions are, therefore, expected and give an estimate on how to understand the qualitative behavior of the implemented model.

11. Conclusions In this paper, we presented an extension of the grid-free SPH method suitable for viscoelastic problems. The main ingredient of our approach is to add a differential constitutive (evolution) equation for the components of the anisotropic stress tensor. The viscosity is treated in a natural way by solving the continuum conservation laws in the SPH formalism without taking into account an artificial pressure tensor in all momentum equations. Neither nested sums nor second-order differentiations of the kernel are involved in the time loop. The differential constitutive equation for the stress implemented in this work for the case of a simple isothermal flow should be regarded as a prototype. The extended SPH method is equally applicable to study not only viscoelastic fluids but also plastic flow and solid friction. Non-isothermal flow can be modeled in the same spirit by using the heat flux vector as an additional independent variable, i.e. by considering an energy balance equation together with a model constitutive equation for the heat flux. Results for the isothermal viscoelastic Jaumann–Maxwell model have been presented considering a basic 2D channel flow test. Comparison of relaxation processes of an initially sinusoidal velocity profile with the available analytical solution show agreement (within expected extent) between numerical experiment and theory. Based on these findings from the study of simple flow, which confirms the suitability of reported implementation details, and in particular of a weak modification at the boundaries,

50

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

the geometry-free, and flow-type-free method presented here can be used to study applied problems. Here, models appropriate to describe the material in simple viscoelastic flows should be used. We might have illustrated the scheme by using other models, such as the codeformational Maxwell model. The method can be equally applied to complex geometries, such as the flow around a cylinder or driven cavity flow for which experimental results are available [44,45].

Acknowledgements Financial support provided by the Deutsche Forschungsgemeinschaft (DFG) via the Sonderforschungsbereich (Sfb) 605 ‘Elementarreibereignisse’ is gratefully acknowledged. This research was supported in part by the National Science Foundation under grant no. PHY99-07949. References [1] M.J. Crochet, A.R. Davies, K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, 1984. [2] R. Keunings, Simulation of viscoelastic fluid flow, in: C.L. Tucker III (Ed.), Fundamentals of Computer Modeling for Polymer Processing, Hanser, Munich, 1989, p. 402. [3] F.P.T. Baaijens, Mixed finite element methods for viscoelastic flow analysis: a review, J. Non-Newtonion Fluid Mech. 79 (1998) 361. [4] M. Laso, H.C. Öttinger, Calculation of viscoelastic flow using molecular models: the CONNFFESSIT approach, J. Non-Newtonion Fluid Mech. 47 (1993) 1. [5] K. Feigl, M. Laso, H.C. Öttinger, The CONNFFESSIT approach for solving a two-dimensional viscoelastic fluid problem, Macromolecules 28 (1995) 3261. [6] M. Laso, M. Picasso, H.C. Öttinger, 2D time-dependent viscoelastic flow calculations using CONNFFESSIT, AIChE J. 43 (1997) 877. [7] M.A. Hulsen, A.P.G. van Heel, B.H.A.A. van den Brule, Simulation of vicoelastic flows using Brownian configuration fields, J. Non-Newtonion Fluid Mech. 70 (1997) 79. [8] T.W. Bell, G.H. Nyland, J.J. de Pablo, M.D. Graham, Combined Brownian dynamics and spectral simulation of the recovery of polymeric fluids after shear flow, Macromolecules 30 (1997) 1806. [9] P. Halin, G. Lielens, R. Keunings, V. Legat, The Lagrangian particle method for macroscopic micro–macro viscoelastic flow computations, J. Non-Newtonion Fluid Mech. 79 (1998) 387. [10] M.A. Hulsen, A.P.G. van Heel, B.H.A.A. van den Brule, Simulation of viscoelastic flows using Brownian configuration fields, J. Non-Newtonion Fluid Mech. 70 (1997) 79. [11] L.B. Lucy, A numerical approach to the testing of the fission hypothesis, Astronom. J. 83 (1977) 1013. [12] R.A. Gingold, J.J. Monaghan, Smoothed particle hydrodynamics theory and application to non-spherical stars, Mon. Not. R. Astronom. Soc. 181 (1977) 375. [13] H. Takeda, S.M. Miyama, M. Sekiya, Numerical simulation of viscous flow by smoothed particle hydrodynamics, Prog. Theor. Phys. 92 (5) (1994) 939. [14] J.J. Monaghan, Simulating free surface flows with SPH, J. Comput. Phys. 110 (1994) 399. [15] S. Siegler, H. Riffert, Smoothed particle hydrodynamics simulations of ultrarelativistic shocks with artificial viscosity, Astrophys. J. 531 (2000) 1053. [16] L.D. Libersky, A.G. Petschek, T.C. Carney, J.R. Hipp, F.A. Allahdadi, High strain Lagrangian hydrodynamics, J. Comput. Phys. 109 (1993) 67. [17] A.G. Petschek, L.D. Libersky, Cylindrical smoothed particle hydrodynamics, J. Comput. Phys. 109 (1993) 76. [18] L.D. Libersky, A.G. Petschek, Smooth Particle hydrodynamics with strength of materials, in: H.E. Trease, M.J. Fritts, W.P. Crowley (Eds.), Advances in the Free-Lagrange Method, Moran, WY, June 1990; Lecturer Notes in Physics, Vol. 395, Springer, New York, 1991, p. 248.

M. Ellero et al. / J. Non-Newtonian Fluid Mech. 105 (2002) 35–51

51

[19] C. Maveyraud, W. Benz, A. Sornette, D. Sornette, Solid friction at high sliding velocities: an explicit 3D dynamical smoothed particle hydrodynamics approach, J. Geophys. Res. 104 (1999) 28769. [20] R. Speith, H. Riffert, H. Ruder, Numerical fluid dynamics in astrophysics with smoothed particle hydrodynamics, in: H.-J. Bungartz, F. Durst, C. Zenger (Eds.), High Performance Scientific and Engineering Computing, Lecture Notes in Computational Science and Engineering, 1999, p. 417. [21] O. Flebbe, S. Münzel, H. Herold, H. Riffert, H. Rüder, Smoothed particle hydrodynamics: physical viscosity and the simulation of accretion disk, Astrophys. J. 431 (1994) 754. [22] C.V. Chaubal, A. Srinivasan, Ö. Egecioglu, L.G. Leal, Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems. Part 1. Method, J. Non-Newtonion Fluid Mech. 70 (1997) 125. [23] T. Bubeck, M. Hipp, S. Hüttemann, S. Kunze, M. Ritt, W. Rosenstiel, H. Ruder, R. Speith, Parallel SPH on Cray T3E and NEC SX-4 using DTS, in: E. Krause, W. Jäger (Eds.), Proceedings of the Conference on High Performance Computing in Science and Engineering, 1998, 1999, p. 396. [24] P. Espanol, M. Serrano, H.C. Öttinger, Thermodynamically admissible form for discrete hydrodynamics, Phys. Rev. Lett. 83 (1999) 4542. [25] P. Espanol, Fluid particle dynamics and smoothed particle dynamics, Europhys. Lett. 39 (1997) 606; P. Espanol, Fluid particle model, Phys. Rev. E 57 (1998) 2930. [26] B.I.M. ten Bosch, On a extension of dissipative particle dynamics for viscoelastic flow modelling, J. Non-Newtonion Fluid Mech. 83 (1999) 231. [27] M. Serrano, P. Espanol, Thermodynamically consistent mesoscopic fluid particle model, Phys. Rev. E 64 (2001) 046115. [28] W. Benz, Smoothed particle hydrodynamics: a review, in: J.R. Buchler (Ed.), The Numerical Modelling of Nonlinear Stellar Pulsation: Problems and Prospects, Les Arcs, France, March 1989, NATO ASI Series, Vol. 302, Kluwer Academic Publishers, Boston, MA, 1990, p. 269. [29] J. Monaghan, Smoothed particle hydrodynamics, Annu. Rev. Astronom. Astrophys. 30 (1992) 543. [30] J.J. Monaghan, An introduction to SPH, Comp. Phys. Commun. 48 (1988) 89. [31] O. Kum, W.G. Hoover, Viscous conducting flows with smooth-particle applied mechanics, Phys. Rev. E 52 (1995) 4889. [32] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, McGraw-Hill, New York, 1981. [33] G.K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, 1967. [34] J.P. Morris, P.J. Fox, Y. Zhu, Modeling low Reynolds number incompressible flows using SPH, J. Comput. Phys. 136 (1997) 214. [35] R.L. Panton, Incompressible Flow, Wiley, New York, 1984. [36] J.H. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, Springer, Berlin, 1997. [37] R.B. Bird, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, Vol. 1, Wiley, New York, 1976. [38] J. Von Neumann, R.D. Richtmeyer, A method for the numerical calculation of hydrodynamic shocks, J. Appl. Phys. 21 (1950) 232. [39] J.P. Morris, J.J. Monaghan, A switch to reduce SPH viscosity, J. Comput. Phys. 136 (1997) 41. [40] W.G. Hoover, S. Hess, Equilibrium and nonequilibrium thermomechanics for an effective pair potential used in smooth particle applied mechanics, Physica A 231 (1996) 425. [41] H.A. Posch, W.G. Hoover, O. Kum, Steady-state shear flows via nonequilibrium molecular dynamics and smooth-particle applied mechanics, Phys. Rev. E 52 (1995) 1711. [42] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, Cambridges, 1992, p. 847. [43] D.D. Joseph, Fluid Dynamics of Viscoelastic Liquids, Series in Applied Mathematics, Vol. 84, Springer, Berlin, 1990 [44] P. Pakdel, S.H. Spiegelberg, G.H. McKinley, Cavity flows of elastic liquids: two-dimensional flows, Phys. Fluids 9 (11) (1997) 3123. [45] P. Pakdel, G.H. McKinley, Cavity flows of elastic liquids: purely elastic instabilities, Phys. Fluids 10 (5) (1998) 1058.