- Email: [email protected]
Numerical simulation of thermal convection with Maxwellian viscoelasticity
Journal
of Non-Newtonian
Elsevier Science Publishers
Fluid Mechanics,
39 (1991)
Numerical simulation of thermal with Maxwellian viscoelasticity Helmut Max-Planck
(Received
67
67-88
B.V., Amsterdam
convection
Harder Institut ftir Chemie,
Saarstrasse
23, Postfach
April 12, 1990; in revised form September
3060,
6500 Mainz
(Germany)
10, 1990)
Abstract This paper investigates the geophysically relevant problem of thermal convection in a Maxwell medium. The case of convection in a 2-dimensional cavity heated from below is solved numerically for infinite Prandtl number and constant material parameters. Within a finite difference context, different approximation schemes for the advective terms in the material equations are explored. Using an upwind approximation, solutions for an upper convected Maxwell model up to Deborah numbers De = 1.0 are possible, where the vigour of the convective flow is reduced to approximately 50% of its Newtonian equivalent. At this De the normal stress pattern seems to become singular at the boundaries. For De < 0.5 it is shown that the results are mainly controlled by the choice of the free parameters, a, of a Johnson-Segalman model. Keywords:
convection;
finite differences;
lithosphere;
mantle;
viscoelasticity
1. Introduction Since the pioneering experiments of Benard, thermal convection has become a familiar problem in many branches of fluid dynamics. Besides engineering applications, such as insulation problems in reactor cooling, convection has attracted great interest among geophysicists, because convective motion undoubtedly governs the dynamics of the Earth’s mantle. Within the Earth’s mantle heat is produced by the decay of radioactive elements, mainly uranium and thorium, thus driving convective motion in the mantle, which is manifest at the surface by tectonic movement of plates. Of course, a necessary precondition is that the Earth behaves as a viscous fluid, at least on long time scales. In fact the relaxation time of the mantle 0377-0257/91/$03.50
0 1991 - Elsevier Science Publishers
B.V.
68 below 100 km depth can be estimated to be of the order of lo3 years. Compared to the characteristic time of the flow, which is roughly lo* years, the assumption of a viscous mantle seems at first glance to be correct. However, since viscosity and relaxation time in the Earth are strongly temperature dependent and the temperature varies within the mantle by at least 3000 K, we can expect that in the cooler, uppermost part of the Earth, the lithosphere (roughly 100 km thick), the Maxwell time is comparable to or exceeds the characteristic time of the flow. Thus the rheology of the Earth should be more correctly described as viscoelastic, especially since all direct observable surface effects of the motion in the deep interior are dominated by the rheology of the cool lithosphere. Therefore the problem of thermal convection in a viscoelastic medium is of geophysical importance. Unfortunately this problem has not attracted much attention in the fluid dynamics community. Apart from studies on the stability problem at the onset of convection [l], previous work on the numerical simulation of thermal convection to my knowledge is restricted to a single study of Ivins et al. [2]. Since they investigated only the case of low Rayleigh numbers and the range of Deborah numbers was so small that Newtonian and viscoelastic flow differ very little, more work on this problem is required. Here only the most fundamental problem of thermal convection will be discussed: the case of a homogeneous Maxwell fluid in a two-dimensional rectangular box heated from below at infinite Prandtl number, i.e. inertial effects are neglected, which is a reasonable assumption for this application. Some generalisations to cover the application to the Earth-such as temperature-dependent material parameters-are found in a companion paper [3]. By choosing a Maxwell rheology it is intended to minimize the number of material parameters in the rheological model. More complicated rheologies like the Jeffreys model and its non-linear generalisations, i.e. Petrov and Giesekus models, may offer some advantages for numerical simulation [4] but since it would be difficult to find appropriate values for the additional material parameters in the Earth the Maxwell model is preferred here. It is well known [5] that numerical simulation of a flow in a Maxwell medium is a difficult task; most known algorithms fail to converge for fairly low Deborah numbers. This may be a consequence of a change of type in the equation of motion as discussed by Joseph et al. [6]. For the case studied here-upper convected Maxwell model and inertial effects neglected- this possibility can be ruled out for physical reasons [6]. However, a change of type may occur artificially due to numerical errors in the calculation of the stress fields. Recent work on the numerical simulation of viscoelastic flows through a contraction [7] has shown that the range of Deborah numbers for which solutions can be obtained is enhanced considerably if special tech-
69 niques are used, which include an asymmetric upwind approximation for the advective terms in the constitutive equations. Therefore progress in this respect seems possible by the introduction of new numerical methods in the treatment of the hyperbolic constitutive equations. The convection problem has certain advantages compared to other more frequently studied flow problems which makes it attractive for viscoelastic flow simulation: first very accurate (and well behaved) Newtonian solutions are known from the benchmark studies of de Vahl Davies and Jones [8] for finite Prandtl number and of Blankenbach et al. [9] for infinite Prandtl number. This makes it easier to estimate the accuracy of a numerical code for Newtonian problems. Second, since the flow is confined to an enclosed cavity, boundary conditions are derived naturally from the constitutive equations and it is not necessary to impose additional assumptions at the inor out-flows as has been done for the most intensively studied problem of a flow through a contraction. Finally, no singularities are known to be present in the Newtonian solution in contrast to the driven-cavity problem which has also been considered as a test case for viscoelastic flow simulation [lo]. Although this is true for a flow with low Deborah number (De < 0.4), it will be demonstrated later that for higher Deborah number certain singularities develop in the Maxwellian stress fields of the convective problem. Thermal convection can therefore be regarded-besides its practical importance-as an excellent test case for simulating viscoelastic flows. I will therefore not only try to find how the convection patterns differ between Newtonian and Maxwellian fluids, but also compare the performance of two established finite difference methods to integrate the constitutive equations, which are similiar to techniques already discussed by Crochet et al. [5]. A finite difference approach is preferred here instead of a finite element method because the geometry of the problem is simple (a rectangular box) and because of the ease of implementing and modifying a finite difference (FD) method. After discussing briefly the set of governing equations in Section 2, the numerical methods are outlined in Section 3. The results obtained on various grids are discussed in Section 4, where both the upper convected and the co-rotational Maxwell models are considered. The results are restricted mainly to steady-state solutions. In a single example, however, the time evolution of the flow is studied. Finally, the performance of the numerical methods at higher Deborah number is discussed. 2. Governing equations Only the most basic problem of thermal convection in a Maxwell medium with homogeneous material parameters will be considered here: the case of a
70 creeping flow in a two-dimensional rectangular box heated from below. The Boussinesq approximation is assumed, i.e. density variations enter only in the buoyancy term within the equation of motion. By using a stream function approach the continuity condition is fulfilled implicitly and the pressure is eliminated from the problem. Under these assumptions the system of equations to be solved is the following:
v”#=
Raz
3T
+
s=r-2qd, T +
AZ
g
=v2T--
(1) (2)
= 2qd,
(uv)T,
where T is the temperature, $ is the stream function, u is the velocity and n is the viscosity. Here eqn. (1) corresponds to the equation of motion, eqn. (2) defines the non-Newtonian contribution s to the stress tensor 7, eqn. (3) is the Maxwell constitutive equation and eqn. (4) gives the energy equation. In eqn. (3) 9/9t indicates an invariant time derivative LBr ar _=at +(uV)~+w-u~7-a(d~+~d), .9t
where - 1 I a I 1 is a free parameter, d = (VU + vu')/2 the deformation tensor and w = (vu - our)/2 the vorticity. With the choice a = 1 in eqn. (5) the upper convected derivative is obtained, with a = 0 the co-rotational one. The system of governing equations (l-5) is given in dimensionless form and the variables are scaled as in Christensen [ll]. The problem is then controlled by two dimensionless numbers: the Rayleigh number Ra = agpATh3/( KQ) and the elasticity parameter X = A’K/h2 introduced by Ivins et al. [2], which is the relaxation time A’ scaled by the thermal diffusion time h2/K. Here K denotes the thermal diffusivity and h the height of the cell. Since the time t is scaled by the thermal diffusion time, from the elasticity parameter X the relative importance of elastic effects cannot be estimated. This would require the scaling of the relaxation time by a characteristic (overturn) time of the flow, which is, however, not known in advance. The Deborah numbers De given in this study are therefore calculated a posteriori. The Deborah number is here defined as De = max X 11d 11,where the maximum is taken over the whole convection cell. Another important result of the calculation is the Nusselt number Nu, which gives the mean vertical
71 heat flux in the convective system relative to the heat flux in a purely conductive case, i.e. Nu = - (iU’/ay)/( AT/h). Numerically Nu is calculated as in Schubert and Anderson [12]. By the introduction of the non-Newtonian stress contribution s the familiar biharmonic operator is restored on the left-hand side of eqn. (1) so only the right-hand side differs from the Newtonian case. This is a standard method to decouple the equation of motion (1) and the constitutive equations (3) in a finite difference context [5]. Other non-decoupling techniques, suitable for finite element (FE) methods, have been discussed by Marchal and Crochet [13] and Luo and Tanner [ 141. As already pointed out, the problem does not require special boundary conditions for the stress tensor, since at the boundaries simplified versions of (3) derive from the conditions for the stream function 4. As the velocity components (u, w) relate to $ by u = a#/ay and w = - a$~/iZlx, the investigated condition of an enclosed cavity and symmetric boundaries requires 4 = 0 and a2G/a xi = 0 at the boundaries-here x, denotes the direction normal to the boundary. From this it follows that the normal velocity component, vorticity, shear deformation and shear stress vanish at all of the boundaries. Physically these conditions apply to a convection cell surrounded horizontally by identical cells, whereas at the bottom and top of the convecting layer free slip is possible due to the low viscosity of the surrounding layers. These conditions are more appropriate for the application studied than, for example, no-slip conditions. The temperature is set to T = 0 at the upper boundary and T = 1 at the lower boundary. The sidewalls are symmetric, i.e. i3T/dx = 0 there. As an initial condition for the temperature field the conductive profile plus a small sinusoidal perturbation is assumed. In the viscoelastic cases the initial conditions for stress and temperature are those calculated for a lower elasticity parameter. 3. Numerical
method
To solve the system eqns. (l-5) a finite difference approach is used. At the start of the simulation the band matrix resulting from the discretised biharmonic operator in eqn. (1) is factorized once by Gaussian elimination and then stored. This requires only a small amount of the total numerical work of the whole simulation. However, the necessity to store the factorized matrix limits the number of grid points to less than approximately 100 x 100. Because at Rayleigh numbers Ra 2 lo5 thin thermal boundary layers develop, and in cases with Deborah numbers De 2 0.2 also huge stress gradients exist at the margins, it is desirable to increase the resolution there beyond what would be possible with an equidistant mesh. Therefore the
72 possibility is included to use a stretched grid which is successively refined from the centre of the cell towards the margins in geometrical fashion. Harder [15] has shown both by consideration of the local truncation error and by numerical calculations on grids with different resolution that second-order accuracy is preserved in cases with non-equidistant mesh. All terms of the forcing function on the right-hand side (RHS) of eqn. (1) are approximated by central differences. The temperature equation (4) is solved by an alternating direct implicit (ADI) method using a Peaceman-Rachford split. Both diffusion and advection terms are approximated by central differences. If the time development of the convective pattern is of interest, the time step At should be limited by :he Courant criterion, i.e. At I min(AXi/Ui, Ayj/wj) at all grid points (i, j). This criterion is identical to the stability limit of a simple explicit scheme. If only the steady-state solution is desired, however, it is possible with the ADI method to use longer time steps and a steady state is reached with far fewer iterations than would be possible with an explicit scheme. Usually At I n min( Axi/ui, Ayj/wj) is a good choice, where n is the number of grid points in one direction. A further enlarged At is not advisable since the method is less efficient with a too large At and even numerical instability can occur. For a Newtonian fluid only the system of the equation of motion (1) and the temperature equation (4) has to be solved, since the non-Newtonian terms on the RHS of eqn. (1) vanish. The results obtained with the outlined method for a variety of Newtonian test cases are documented in the benchmark study by Blankenbach et al. [9]. From these results it can be estimated that in problems with Ra = lo4 a mesh of 40 x 40 equidistant grid points is sufficient to achieve an accuracy of at least 0.3% in both velocity and temperature fields. This accuracy should be representative for the Newtonian results given in this study. In order to investigate convection in a Maxwellian medium the numerical method described is generalised. The material equations (3) are solved as a decoupled system using two alternative finite difference (FD) techniques already discussed by Crochet et al. [5]. (a) The upwind method: here the time derivatives in eqn. (5) are approximated by forward difference quotients and the spatial derivatives by first-order upwind quotients, i.e. one-sided difference approximations dependent on the sign of the local velocity components. Introducing the notation 7 = 7x,, v = rYv- 7,,, p = 7Yv+ 7,, and labelling the time level by a superscript, while subscripts refer to grid points, the following scheme arises for the shear stress 7:
13 where ~1 =
$At( ) u( + U)/( Xi - Xi-l)>
~3= iAt( I I+ I +W>/(.Yj-Yj-,>,
+_=tAt( c4 =
IuI-~u)/(~i+~ -xi),
I w I - w)/(Yj+I -Vj)-
iAt(
For the other stress components similiar schemes apply. Since only one-sided differences are used, this scheme is only a first-order approximation of the material equations (3). It is therefore desirable to use a method of higher accuracy. The second technique will fulfil this demand. (b) The DuFort-Frankel method: this scheme was introduced by Townsend [16] to solve viscoelastic material equations. Time as well as spatial derivatives are approximated by central differences. Therefore the discretisation extends over three time levels and has the following form:
=-
cl(
Ti”+l,j
-
'l?l,j
) - c2(Tiyj+l
Tiyj_l)-
-
(1- A) q’
+ xov:,-’
+ Xd,,ap~l~l + 4qdxy,
~:j+’ + 2Xd,,ap~~’
4hw$+l = -‘I(
‘in+l,j
-
‘El,,
- 2hd,,ap~,~1 -4Xd,,a$+’ =
(74
)
-
- 2Xd xxav”-l I,
-
“;j_l)
-
4xLVi;-1 -
- 8qd,,,
+ 2hd,,av~f1
-cl(P~+l,j-P~-l,j
c2(v;j+l
1
x at
1
n-1
“ij
(W
+
)-~2(/~r,j+l
-/.~:,j-l)+
4Xd,,a~~-’ (74
where c1 = 2hU/(Xi+l -x~_~) and c2 = 2hw/(~~+~ -Y~_~). To calculate the new stress components (7, v, p)“+l a (3 x 3) matrix has to be inverted for each grid point. In case of a co-rotational model, i.e. a = 0, the system reduces to a set of (2 X 2) matrices, since for the stress trace p = 0 holds within the whole cell. This algorithm is to my knowledge the only FD technique applied to the solution of Maxwellian material equations which posseses second-order accuracy. Although it can be expected that this method is superior to the first-order upwind method described before, a direct comparison is not available yet. This will be done for an example in the next section.
74 For both methods no rigorous stability analysis is available. If a fixed velocity field is considered, the upwind method requires a time step At < A, whereas the DuFort-Frankel scheme allows longer time steps. If the coupled system eqns. (l-4) are solved, however, only a much smaller time step is possible, At -+c h, provided the Deborah number is high enough to induce a substantial coupling of the equation of motion (1) and the material equations (3). In the numerical calculations a time step At in the range of X/20 to h/5 is used. The iteration scheme for the solution of the non-linear system eqns. (l-4) is of a straightforward Picard type, where from the new stress components (7, v)“+’ the non-Newtonian terms on the RHS of eqn. (1) are calculated by the stress decomposition eqn. (2). A modification, similiar to that of Townsend [16], which consists of a sub-iteration for stresses and stream function, i.e. the calculation of the new stress components (7, v, p)‘+’ is repeated with updated coefficients in eqns. (2) and (3), did not improve the stability of the iteration in the problem considered here. Therefore in general the simpler and computationally inexpensive Picard scheme is preferred. If a time-dependent solution is to be calculated, the more restrictive timestep limit arising from the temperature eqn. (4) and material equations (3) must be used for the solution of the coupled system. This is here the condition At c A, following from the explicit scheme for the material equations. For a typical problem considered in this study, around lo4 iterations are needed with such a method. If only the steady-state solution is of interest, the number of iterations can be reduced by using different ‘pseudo-timesteps’ in both equations. With the choice At, = n min( Ax;/ u,, A yj/wj) in the temperature equation and At,s h in the material equations a steady-state solution could usually be calculated with less than lo2 iterations-an obvious improvement compared to a ‘ truely’ time-marching method. Since no analytical, and only very few published numerical, solutions on viscoelastic convection problems are available, it is not easy to verify the numerical scheme quantitatively. Although a qualitative comparison with the results of Ivins et al. [2], given as iso-line plots, was possible in a case of weak coupling, additional verification is desirable. Therefore the following relations are useful as internal consistence checks. Parmentier et al. [17] have shown that for a convection cell, regardless of the rheology applied, the identity Raj, wT d f = lF r :d df holds, where the integration is taken over the whole volume F of the cell. With the Maxwellian material equations (4) it is easy to show that assuming a steady-state solution this identity can be extended to the integrals over the invariants p = tr{ 7) and J = - det{ T} of the stress tensor, i.e.:
75 Combined with a convergence test on refined grids, a check of these criteria is not only a safeguard against a systematic error in the numerical code, but also allows a rough estimate of the accuracy of the calculated solution. 4. Results and discussion 4.1 Convergence tests and stability of the numerical scheme In order to examine the accuracy of the numerical code and to compare the different methods of solving the material equations, a set of convergence tests has been completed. As test problem a case with the low Rayleigh number Ra = 9487 was chosen. This value exceeds the critical Rayleigh number Ra, at which convection sets in by only a factor 12. Such a low value is chosen primarily to avoid the additional complication of thin boundary layers in the temperature field developing at high Rayleigh numbers, Ra > lOORa,. The aspect ratio of the cell is one and the upper convected Maxwell rheology (a = 1) is assumed, since this choice of the parameter a in the Oldroyd derivative (5) is most frequently studied in computational fluid dynamics. In Fig. 1 the relative change in the maximum of the stream function #,,, in the viscoelastic flow field compared to its Newtonian counterpart is 16
I
1
8
I
L
14
-
12
-
10
-
8
-
6
-
4
-
2
-
0
’ 0.0
I
1 + : upwind, 20x20 l : upwind,kOx40 o : DuFort-Frankel, x : DuFort-Frankel,
I
0.1
20x20 40x40
I
I
0.3
0.2
I
0.4
0.5
De
Fig. 1. Relative change of the maximum of the stream function Newtonian result for various Deborah numbers De, grid resolution scheme in the material equations.
$,, compared to the and different numerical
X
X
Fig. 2. Stress trace p at the left part of the lower boundary (y = 0 and x < 0.4). For x > 0.4 the stress profile is essentially flat. The results are obtained by the upwind scheme and the DuFort-Frankel scheme (dotted profiles) at different grid resolution. (Ra = 9487, De = 0.34, h = 0.0015, a = 1).
displayed with increasing Deborah number. As the grid resolution is increased, the solutions obtained with different numerical methods obviously converge against each other, as expected. The results of Fig. 1 show further that with the DuFort-Frankel scheme higher accuracy is obtained than with the simple upwind scheme. This, of course, follows from the fact that the first method is a second-order approximation whereas the upwind scheme is only of first order. In the problem considered in Fig. 2 the maximum misfit in the integral relations (8) is, with the DuFort-Frankel scheme, 1.1% on a 20 x 20 grid and decreases to 0.3% on a 40 x 40 grid. The corresponding values with the upwind scheme are 5.6% and 3.0%. These results are consistent with the approximation error of both schemes and confirm the methods used. The upwind method, however, has the inherent deficiency of introducing an artificial numerical diffusivity into the material equations, which is responsible for smoothed stress gradients at the boundaries. This is obvious from Fig. 2, where the stress trace p at the lower boundary in the vicinity of the convergent stagnation point is shown, calculated with both schemes on two identical grids. Nevertheless, the numerical diffusivity of the upwind scheme allows solutions to be obtained at considerably higher Deborah numbers. At De > 0.4 the upwind method gives an acceptable solution whereas the
77 2400
stress trzca i
a)
I
Fig. 3. As Fig. 2, but De = 0.42, X = 0.0020.
DuFort-Frankel scheme fails at values De = 0.4. The breakdown of the scheme is manifest by artificial oscillations with grid periodicity in the stress fields (Fig. 3), which increase at higher De and finally the iteration diverges. This failure is restricted to the solution of the coupled system eqns. (l-4) and accompanied by a loss of positive definiteness of the stress tensor. Increasing the grid resolution does not change the situation much as demonstrated in Fig. 3. Even with an additional grid refinement to 80 x 80 points or stretched grids, which increase the resolution near the boundaries further by a factor 5, the problem changes only insignificantly. Similiar behaviour is known for most numerical methods in a wide class of viscoelastic flow problems, but the special case of Maxwellian rheology seems the most difficult to handle [4,5]. Only recently finite element techniques have been introduced, which incorporate a one-sided approximation of the advective terms in the material equations and which-combined with a subdivision of the finite elements into finer sub-elements for approximating the stresses-allow a solution at substantially higher De. It is therefore remarkable that even the established FD upwind improves the range of De where well behaved solutions can be obtained. The discretisation of the convective terms in the material equations seems therefore to be of crucial importance. The solution of high Deborah number problems will be further discussed in Section 4.5, but first some aspects of the convection problem will be considered in the low Deborah number regime (De c 0.4), where both numerical schemes yield acceptable solutions.
78 4.2 Comparison of upper convected/ co-rotational models In the Oldroyd derivative eqn. (5) of the constitutive eqn. (3) the constant a is considered as a free rheological parameter. This approach is sometimes referred to as a Johnson-Segalman model. It is now easy to show that by a transition from a = 1 (upper convected) to a = - 1 (lower convected) the solution remains the same with the only exception that the stress trace p changes sign. Therefore only the two extreme cases a = 1 and a = 0 (co-rotational) will be discussed. The choice of an upper convected Maxwell model is much more popular than the co-rotational model since the latter case exhibits unrealistic rheometrical behaviour at De > 0.5 in a simple shear flow, i.e. a decrease in stress with increasing shear rate and, associated with that, the stress tensor will no longer be positive definite, see e.g. ref. 4. Nevertheless, the co-rotational case is a useful rheological model for De -C0.5. The effect of a different choice of the parameter a is therefore investigated for an example already considered in Fig. 2 (Ra = 9487). As a reference, Fig. 4 displays the Newtonian solution, whereas Fig. 5 gives the solution in the case of an upper convected Maxwell model and Fig. 6 refers to the corresponding co-rotational case. The elasticity parameter is in both cases identical, h = 0.0015, the Deborah numbers are different, however, since De depends on the flow field.
TemDerature
St ream-Fc
Shear-stress
Stress-dif.
t .
Vorticitv
Fig. 4. Temperature T, stream function #, vorticity w, shear stress T and normal stress difference Y for Newtonian rheology and Ra = 9487. Isolines: T(1/7,. . . ,6/7); $( - 3, 6,. . . , - 15); ~(30, 60,. . . ,150); T( - 90, - 60,. . . ,60); V( - 400,0,400, 800).
79 Temperature
Stream-Fct.
Vorticitv
Shear-stress
Stress-dif.
Stress-trace
k2!
Fig. 5. As Fig. 4, but with upper convected Maxwell rheology Additional isolines for the stress trace: ~(75,275,. . ,1075).
For both models characteristic exist: whereas in the co-rotational than its Newtonian counterpart,
(a = 1, A = 0.0015, De = 0.34).
differences to the Newtonian solution case the convective flow is more vigorous the upper convected model yields the
a a+ Temperature
Stream-Fct.
Vorticity
0 0
Shear-stress
Stress-dif.
-H_
a
0
l
Fig. 6. As Fig. 4, but with co-rotational
Maxwell
rheology
(a = 0, h = 0.0015, De = 0.52).
80 opposite effect. Although the differences are hardly visible in the isoline plots of temperature and stream function, the relative deviation from the Newtonian values of the Nusselt number and stream function maximum is (- 2.7, - 5.3)% for a = 1 and (+ 7.7, + 8.8)% for a = 0.These values, as the isoline plots, are calculated on the fine 40 X 40 grid using the DuFortFrankel scheme; high accuracy can thus be expected. The fact that both Maxwell models considered produce opposite deviations from the Newtonian flow, is plausible. A co-rotational Maxwell model is characterised in a simple shear flow by an effective viscosity decreasing with shear rate, whereas in a pure shear flow the effective viscosity is constant and equals the Newtonian values. The upper convected case is different: in a simple shear flow the effective viscosity is constant and equals the Newtonian value whereas in a pure shear flow the effective viscosity increases with strain rate. Therefore it seems reasonable that in a ‘real’ flow, which consists of both pure and simple shear, the co-rotational model will exhibit features of a decreased effective viscosity and the upper convected model similiarly an increased effective viscosity. Correspondingly, the changes in the vorticity field are completely different for both Maxwell models. With a = 1 the deviations from the Newtonian solution are minor at this De-only the maxima are slightly reduced, which is obviously due to the less vigorous flow. In the co-rotational case on the other hand, zones of high vorticity gradients develop at the convergent margins, which are not present in the Newtonian flow and which become larger with increasing Deborah number. In fact some small wiggles in the vorticity have already developed at De = 0.52. This seems to be a consequence of the fact that the ‘critical’ Deborah number De = 0.5 for a simple shear flow is slightly exceeded. Since this is a deficiency of the particular rheological model, both numerical methods fail at higher Deborah numbers. The differences in the stress fields of both Maxwell models are mainly restricted to the stress trace p, which vanish identically for a = 0, as in the Newtonian case. For a = 1,however, in p pronounced boundary layers form in the vicinity of convergent zones of neighbouring convection cells. The normal stress difference Y and the shear stress 7, on the other hand, do not differ much for both Maxwell models. Of course characteristic differences exist if the stress fields are compared with their Newtonian counterparts. A comparison of the shear stresses in Figs. 5 and 6 with Fig. 4 illustrates the fact that in a Maxwellian fluid stresses are advected by the flow. For the same reason large stress gradients build up in the vicinity of convergent margins, whereas the gradients are smoothed at divergent margins. The differences described are accompanied by corresponding effects on the quantities measurable at the upper surface (y = 1) of the convection cell
81 10 Heatflux
’
0
0.0
a)
I 0.2
0.4
0.6
0.8
1.0
X 80 Surface velocity
bl
70 60 -
0.0
0.2
0.4
0.6
0.8
1.0
X
Fig. 7. (a) Heat flux and (b) horizontal velocity at the upper surface (y =l) for x = 0.0015 (a = 0, De = 0.52 and a =l, De = 0.35). The dotted profiles give the Newtonian results.
-the most important are the surface heat flux and the horizontal surface velocity. As shown in Fig. 7 these quantities are mainly controlled by the vigour of the flow. With a = 0 the heat flux at the surface is enlarged, whereas with a = 1 it is decreased. The same is true for the surface velocity. Summarising the result, we can say that with Deborah numbers De < 0.5 the impact of the viscoelasticity is small and the changes of the Maxwellian flow compared to its Newtonian equivalent are mainly controlled by the choice of the parameter a in the Oldroyd derivative eqn. (5).
82 4.3 Flow at high Rayleigh numbers From various numerical simulations, e.g. McKenzie et al. [18], it is well known that convection at high Rayleigh numbers differs significantly from the flow at slightly over-critical Rayleigh numbers. The main new features are narrow boundary layers in the temperature and high vorticity gradients at the cell boundaries, In order to investigate the influence of these new features for a Maxwell medium, models were calculated for Rayleigh numbers Ra = (12, 40, 120) X Ra,. This range includes the three characteristic domains of slightly overcritical, intermediate and boundary layer flow [19]. The results of Fig. 8 demonstrate that the changes between convection in a Maxwellian and a Newtonian medium are mainly independent of the Rayleigh number, at least at this low Deborah number. Similiar conclusions were reached by Ivins et al. [2]. However, in their models the Rayleigh number was less than Ra < 30Ra, and the Deborah numbers were so low in the models with higher Rayleigh number, that only a change of less than 0.1% in the maximum stream function can be observed. However, in addition, in Fig. 8 a small tendency at De > 0.3 is visible to reduce the variation of $,, with increasing Rayleigh number.
12
10
I
I
I
I
1
-AlCIr,m[%,l
+
-
a -
*
+ : ua487
.
* : Ra=?1623
X
De Fig. 8. Relative number De.
changes
of $,,
dependent
on the Rayleigh
number
Ra and
Deborah
83 4.4 Time-dependent solutions Up to this point only stationary solutions were considered. Since it is possible that the stationarity of the solutions was artificially induced by the numerical technique, for one example the flow has been calculated with a time-marching algorithm. The time step is At = 0.1X and the initial condition for both temperature and stress fields is the Newtonian solution. The time evolution given in Fig. 9 shows that the calculation converges towards the stationary solution independent of the numerical time step. From Fig. 9 it is obvious that for convection in a Maxwellian medium two characteristic times are relevant: first the relaxation time X and second the transit time t,, which is defined as the time needed for a fluid particle to cross a cell with height h at a mean velocity (u), i.e. t, = (u)/h. The width of the initial drop scales with A, whereas the long-term evolution is controlled by the transit time t,. This typical behaviour is unchanged by varying the time step, grid resolution or approximation of the stress advection terms. The superimposed short scaled variations at times t -c 0.01, however, are highly dependent on these features, hence it must be
I
17.4 /-
I
F-Y-
1
ii 16.8
1
16.6 P,
-1 De=0.34
16.4 tL
16.2 t\
15.8
1
.JL
......,..... ..... ,.
/
"A
0.00
' tll
I
I 0.02
-I
I 0.04
t 0.06
I 0.08
I 0.10
0.12
time Fig. 9. Time evolution of #,, in the case Ra = 9487, De = 0.34, X = 0.0015. The initial condition is the Newtonian solution. The steady-state solution is indicated by the horizontal line. The bars give the relaxation time X and the transit time t,, defined in the text.
84 suspected nique.
that this aspect is an artefact
introduced
by the numerical
tech-
4.5 Convection at high Deborah numbers In the foregoing discussion the range of Deborah numbers has been restricted to De < 0.5. For an upper convected Maxwell model it is possible to extend the simulation to substantially higher values of De provided the upwind approximation is used in the material equations. In Fig. 10 the steady-state solution is given for a similiar case already discussed in Fig. 5, but now with De = 1.06 and X = 0.009. As in the case with low Deborah number (De = 0.34, X = 0.0015, Fig. 5) the convective flow is less vigorous than in the Newtonian case, but the effect is much more pronounced here. For example the Nusselt number Nu, the stream function maximum #,, and the mean upper surface velocity (us) are reduced to (69, 59, 39)s of the corresponding Newtonian solution. The flow is thus effectively damped at higher De. This of course reduces the strain rates within the cell, which explains why De is increasing more slowly
Temperature
Stream-Fct.
Vorticity
+ + + 1 [ Lilizl -B_ +
:i I3
Shear-stress
Stress-dif.
(1
(J
Stress-trace
a
u
u
-
0
+
a
Fig. 10. As Fig. 4, but with upper convected Maxwell rheology (A = 0.0090, De = 1.065, a = 1). Isolines: T(1/7,. . . ,6/7); $( - 1.25, - 2.50,. . . , - 8.75); w( - 80, - 40, 0,. . . ,120); T( - 80, - 40, 0,. . . ) 120); V(- 8000, - 4000, - 2000, - 1000, 0,. . . ) + 8000); ~(125,250, 500 ,..., 16000).
85 than the elasticity parameter X. As the stream lines in Fig. 10 show, the damping effect is most pronounced in the vicinity of the cell boundaries. Therefore the mean surface velocity is more effected by the increased elasticity of the flow than the other cited quantities. A further new feature of the solution is the broad zone of negative vorticity around the boundary. In isoviscous Newtonian convection such a zone exists only in cases with no-slip boundary conditions. Similarly there is a major change in the stress fields. At first glance it seems that the effects visible at lower Deborah number are only enhanced: due to stress advection the shear stress pattern is now rotated by almost 90” relative to the Newtonian solution and the maxima of the normal stress patterns are squeezed towards the boundary. However, although a refined mesh of 80 x 80 non-equidistant grid points is used, it is impossible to resolve the normal stress pattern near the boundary, where a singularity seems to exist. This is a fundamentally new feature at higher De. For low De grid refinement yields convergent solutions all over the cell, including the boundaries (Fig. 2). At the higher Deborah number De = 1.065 this is no longer true. If the grid is refined, the normal stresses diverge along the boundaries starting from the stagnation point at the cell corner (Figs. 11(a) and (c)). Since the numerical values at the boundary are approximately inversely proportional to the mesh size, it seems that an artificial numerical diffusivity in the material equationsas introduced by the first-order upwind approximation-is needed to cope with such singularity. This would explain, why the DuFort-Frankel scheme fails at higher De. In spite of the singularity at the boundary the solution converges in the interior of the cell. As an example this is demonstrated in Figs. 11(c) and (d), where the stress component p calculated at various grids is shown just one grid point (on the coarsest grid) apart from the boundary. Here grid refinement yields a convergent solution. Thus convergence in the interior of the cell is not affected by the singularity at the boundary. Since all the models presented in this chapter are calculated by the first-order upwind scheme, only limited accuracy can be expected. In fact the maximal numerical error in the integral identities (8) is (26, 15, 9.6)s using the three stretched grids of Figs. 11(b) and (d). Although an error in the range of 10% is still considerable for an 80 X 80 mesh, which includes grid refinement towards the boundaries, the main qualitative features of the solutions are not affected by the mesh size. The solution displayed in Fig. 10 seems, therefore, to be a reliable approximation to the Maxwellian flow problem considered. A further but limited, increase in the Deborah number is possible, although, the misfit in the integral identities eqn. (8) then grows rapidly. Using the finest grid, the upwind method finally fails, if the elasticity parameter h (which is the input parameter) is further increased by
I
x*x
r,
I ni
l
,
x
%"
)*"r
.x
,
1200
%"I,
stress trace
800
%x
Il.
4OxkO
. : grid
:
1600
: . +.
*. .
..
80x80
,
grid
x :
X.
* : pr1.i20: 20
110x10
Fig. 11. (a) Stress trace JJ at the left boundary trace p close to the left boundary at x = l/20 parameters are identical to Fig. 9.
0.0 0
0.2 -
0.4 -
100
o’nI
l.O(
0.)
Y.
r
0.6
grid
.: .
.
' x
5oa
*
x
x
x
x
10000
15Oa
20000 25000
. : prld 20x20
. : prld 40x40
x : grl6 BOX80
I---* *
.. x . xX
30000
for three equidistant grids; (b) as (a), but grids are refined towards the boundaries; (c) stress (grids are equidistant); (d) as (c), but at x = 0.011 and non-equidistant grids. All other model
2000
c)
0
0.6
0.8
1.0
87 a factor of 2. In this case similiar problems arise as with the DuFort-Frankel scheme at lower De. It is tempting to speculate how the performance of the algorithm could be improved for high-De problems. Obviously the presence of a normal stress singularity starting at a stagnation point and evolving in the upstream direction is of crucial importance. But in cases as considered here, where in one cell only one convection roll exists, no internal stagnation points are present and singularities evolve only at the cell boundaries starting from the corner points. This could be avoided by a discretisation suggested by Gatski and Lumley [20], where the normal stresses are approximated at midpoints of grid cells, i.e. (v, P)i+i/2,j+i/2> and the velocity components at the centres of grid cell boundaries, i.e. ~~+i,~, j, IJ_ i,*,, and w equivalent. The advantage is that there is no grid point for the normal stresses where all velocity components vanish. Therefore no stagnation point would coincide with a grid point and a singular normal stress should be avoided. Although such an approximation may not be favourable in general, it appears very suitable for a convection problem of an enclosed cavity with symmetric boundaries. For such a problem the method of Gatski and Lumley [20] should be reconsidered and compared with the more straightforward approach used here. Nonetheless, convergence for an unlimited range of De cannot be expected with such a modified method. It should be remembered, however, that even with the simplest method used here Maxwellian convection problems can be solved with a flow which differs very significantly from its Newtonian counterpart. 5. Concluding remarks This study has demonstrated the major differences between convection with Maxwellian versus Newtonian rheology and shown that thermal convection is a very suitable test case for numerical methods simulating viscoelastic flow. It has been possible to extend the simulation up to Deborah numbers De = 1.0, which is sufficient to induce significant changes in the flow fields. A main new feature at high De is the presence of a normal stress singularity along the boundary, which is absent in Newtonian or in low-De flow. Since this singularity is starting from the stagnation points at the cell corners, this behaviour is presumably related to the well known singularity of the upper convected Maxwell model in a pure shear flow [4,5]. It is common experience [4,7,13] that rheological models with more realistic response are numerically easier to handle. Examples are the Jeffreys model and its non-linear generalisations, i.e. the Giesekus and Leonov models. Certainly the numerical approach used can be extended to such models, but for geophysical applications these models are up to now only of
88
limited use, since it would be difficult to relate the additional model parameters to the Earth. In further studies [3,15] it is thus preferred to investigate first the geophysically better understood effects of material parameters which depend on temperature and hydrostatic pressure. This approach can also remove the singularities arising from stagnation points and the numerical simulation is more stable than in the cases studied here. Further studied may show that the Jeffreys model and related ones are relevant to the Earth. From a fluid dynamical point of view, thermal convection with such rheological models is certainly of interest. Acknowledgements I wish to thank W.R. Jacoby and U. Christensen for suggesting this problem to me and for many helpful discussions during this work. Financial support by the Deutsche Forschungsgemeinschaft is gratefully acknowledged. References 1 I.A. Eltayeb, Proc. R. Sot. London, Ser., A, 356 (1977) 161-176. 2 E.R. Ivins, T.W.J. Unti and R.J. Phillips, Geophys. Astrophys. Fluid Dyn., 22 (1982) 103-132. 3 H. Harder, U. Christensen and W. Jacoby, in preparation. 4 J. van der Zanden and M. Hulsen, J. Non-Newtonian Fluid Mech., 29 (1988) 93-117. 5 M.J. Crochet, A.R. Davies and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, 1984. 6 D.D. Joseph, M. Renardy and J.C. Saut, Arch. Ration. Mech. Anal., 87 (1985) 213-251. 7 B. Debbaut, J.M. Marchal and M.J. Crochet, J. Non-Newtonian Fluid Mech., 29 (1988) 119-146. 8 G. de Vahl Davies and I.P. Jones, Int. J. Num. Meth. Fluids, 3 (1983) 227-248. 9 B. Blankenbach, F. Busse, U. Christensen, L. Cserepes, D. Gunkel, U. Hansen, H. Harder, G. Jarvis, M. Koch, G. Marquart, D. Moore, P. Olson, H. Schmeling and T. Schnaubelt, Geophys. J. Int., 98 (1989) 23-38. 10 M.A. Mendelson, P.-W. Jeh, R.A. Brown and R.C. Armstrong, J. Non-Newtonian Fluid Mech., 10 (1982) 31-54. 11 U. Christensen, Geophys. J.R. Astron. Sot., 77 (1984) 343-384. 12 G. Schubert and C.A. Anderson, Geophys. J.R. Astron. Sot., 80 (1985) 575-601. 13 J.M. MarchaI and M.J. Crochet, J. Non-Newtonian Fluid Mech., 26 (1987) 77-114. 14 X.-L. Luo and R.I. Tanner, J. Non-Newtonian Fluid Mech., 31 (1989) 143-162. 15 H. Harder, Dissertation, Johannes Gutenberg-Universitat Mainz, Department of Geosciences, 1988 (in German). 16 P. Townsend, J. Non-Newtonian Fluid Mech., 14 (1984) 265-278. 17 E.M. Parmentier, D.C. Turcotte and K.E. Torrance, J. Geophys. Res., 81(1976) 183991846. 18 D.P. McKenzie, J.M. Roberts and N.O. Weiss, J. Fluid Mech., 62 (1974) 465-538. 19 G.T. Jarvis and W.R. Peltier, Geophys. J.R. Astron. Sot., 68 (1982) 389-427. 20 T.B. Gatski and J.L. Lumley, J. Comput. Phys., 27 (1978) 42-70.
of Non-Newtonian
Elsevier Science Publishers
Fluid Mechanics,
39 (1991)
Numerical simulation of thermal with Maxwellian viscoelasticity Helmut Max-Planck
(Received
67
67-88
B.V., Amsterdam
convection
Harder Institut ftir Chemie,
Saarstrasse
23, Postfach
April 12, 1990; in revised form September
3060,
6500 Mainz
(Germany)
10, 1990)
Abstract This paper investigates the geophysically relevant problem of thermal convection in a Maxwell medium. The case of convection in a 2-dimensional cavity heated from below is solved numerically for infinite Prandtl number and constant material parameters. Within a finite difference context, different approximation schemes for the advective terms in the material equations are explored. Using an upwind approximation, solutions for an upper convected Maxwell model up to Deborah numbers De = 1.0 are possible, where the vigour of the convective flow is reduced to approximately 50% of its Newtonian equivalent. At this De the normal stress pattern seems to become singular at the boundaries. For De < 0.5 it is shown that the results are mainly controlled by the choice of the free parameters, a, of a Johnson-Segalman model. Keywords:
convection;
finite differences;
lithosphere;
mantle;
viscoelasticity
1. Introduction Since the pioneering experiments of Benard, thermal convection has become a familiar problem in many branches of fluid dynamics. Besides engineering applications, such as insulation problems in reactor cooling, convection has attracted great interest among geophysicists, because convective motion undoubtedly governs the dynamics of the Earth’s mantle. Within the Earth’s mantle heat is produced by the decay of radioactive elements, mainly uranium and thorium, thus driving convective motion in the mantle, which is manifest at the surface by tectonic movement of plates. Of course, a necessary precondition is that the Earth behaves as a viscous fluid, at least on long time scales. In fact the relaxation time of the mantle 0377-0257/91/$03.50
0 1991 - Elsevier Science Publishers
B.V.
68 below 100 km depth can be estimated to be of the order of lo3 years. Compared to the characteristic time of the flow, which is roughly lo* years, the assumption of a viscous mantle seems at first glance to be correct. However, since viscosity and relaxation time in the Earth are strongly temperature dependent and the temperature varies within the mantle by at least 3000 K, we can expect that in the cooler, uppermost part of the Earth, the lithosphere (roughly 100 km thick), the Maxwell time is comparable to or exceeds the characteristic time of the flow. Thus the rheology of the Earth should be more correctly described as viscoelastic, especially since all direct observable surface effects of the motion in the deep interior are dominated by the rheology of the cool lithosphere. Therefore the problem of thermal convection in a viscoelastic medium is of geophysical importance. Unfortunately this problem has not attracted much attention in the fluid dynamics community. Apart from studies on the stability problem at the onset of convection [l], previous work on the numerical simulation of thermal convection to my knowledge is restricted to a single study of Ivins et al. [2]. Since they investigated only the case of low Rayleigh numbers and the range of Deborah numbers was so small that Newtonian and viscoelastic flow differ very little, more work on this problem is required. Here only the most fundamental problem of thermal convection will be discussed: the case of a homogeneous Maxwell fluid in a two-dimensional rectangular box heated from below at infinite Prandtl number, i.e. inertial effects are neglected, which is a reasonable assumption for this application. Some generalisations to cover the application to the Earth-such as temperature-dependent material parameters-are found in a companion paper [3]. By choosing a Maxwell rheology it is intended to minimize the number of material parameters in the rheological model. More complicated rheologies like the Jeffreys model and its non-linear generalisations, i.e. Petrov and Giesekus models, may offer some advantages for numerical simulation [4] but since it would be difficult to find appropriate values for the additional material parameters in the Earth the Maxwell model is preferred here. It is well known [5] that numerical simulation of a flow in a Maxwell medium is a difficult task; most known algorithms fail to converge for fairly low Deborah numbers. This may be a consequence of a change of type in the equation of motion as discussed by Joseph et al. [6]. For the case studied here-upper convected Maxwell model and inertial effects neglected- this possibility can be ruled out for physical reasons [6]. However, a change of type may occur artificially due to numerical errors in the calculation of the stress fields. Recent work on the numerical simulation of viscoelastic flows through a contraction [7] has shown that the range of Deborah numbers for which solutions can be obtained is enhanced considerably if special tech-
69 niques are used, which include an asymmetric upwind approximation for the advective terms in the constitutive equations. Therefore progress in this respect seems possible by the introduction of new numerical methods in the treatment of the hyperbolic constitutive equations. The convection problem has certain advantages compared to other more frequently studied flow problems which makes it attractive for viscoelastic flow simulation: first very accurate (and well behaved) Newtonian solutions are known from the benchmark studies of de Vahl Davies and Jones [8] for finite Prandtl number and of Blankenbach et al. [9] for infinite Prandtl number. This makes it easier to estimate the accuracy of a numerical code for Newtonian problems. Second, since the flow is confined to an enclosed cavity, boundary conditions are derived naturally from the constitutive equations and it is not necessary to impose additional assumptions at the inor out-flows as has been done for the most intensively studied problem of a flow through a contraction. Finally, no singularities are known to be present in the Newtonian solution in contrast to the driven-cavity problem which has also been considered as a test case for viscoelastic flow simulation [lo]. Although this is true for a flow with low Deborah number (De < 0.4), it will be demonstrated later that for higher Deborah number certain singularities develop in the Maxwellian stress fields of the convective problem. Thermal convection can therefore be regarded-besides its practical importance-as an excellent test case for simulating viscoelastic flows. I will therefore not only try to find how the convection patterns differ between Newtonian and Maxwellian fluids, but also compare the performance of two established finite difference methods to integrate the constitutive equations, which are similiar to techniques already discussed by Crochet et al. [5]. A finite difference approach is preferred here instead of a finite element method because the geometry of the problem is simple (a rectangular box) and because of the ease of implementing and modifying a finite difference (FD) method. After discussing briefly the set of governing equations in Section 2, the numerical methods are outlined in Section 3. The results obtained on various grids are discussed in Section 4, where both the upper convected and the co-rotational Maxwell models are considered. The results are restricted mainly to steady-state solutions. In a single example, however, the time evolution of the flow is studied. Finally, the performance of the numerical methods at higher Deborah number is discussed. 2. Governing equations Only the most basic problem of thermal convection in a Maxwell medium with homogeneous material parameters will be considered here: the case of a
70 creeping flow in a two-dimensional rectangular box heated from below. The Boussinesq approximation is assumed, i.e. density variations enter only in the buoyancy term within the equation of motion. By using a stream function approach the continuity condition is fulfilled implicitly and the pressure is eliminated from the problem. Under these assumptions the system of equations to be solved is the following:
v”#=
Raz
3T
+
s=r-2qd, T +
AZ
g
=v2T--
(1) (2)
= 2qd,
(uv)T,
where T is the temperature, $ is the stream function, u is the velocity and n is the viscosity. Here eqn. (1) corresponds to the equation of motion, eqn. (2) defines the non-Newtonian contribution s to the stress tensor 7, eqn. (3) is the Maxwell constitutive equation and eqn. (4) gives the energy equation. In eqn. (3) 9/9t indicates an invariant time derivative LBr ar _=at +(uV)~+w-u~7-a(d~+~d), .9t
where - 1 I a I 1 is a free parameter, d = (VU + vu')/2 the deformation tensor and w = (vu - our)/2 the vorticity. With the choice a = 1 in eqn. (5) the upper convected derivative is obtained, with a = 0 the co-rotational one. The system of governing equations (l-5) is given in dimensionless form and the variables are scaled as in Christensen [ll]. The problem is then controlled by two dimensionless numbers: the Rayleigh number Ra = agpATh3/( KQ) and the elasticity parameter X = A’K/h2 introduced by Ivins et al. [2], which is the relaxation time A’ scaled by the thermal diffusion time h2/K. Here K denotes the thermal diffusivity and h the height of the cell. Since the time t is scaled by the thermal diffusion time, from the elasticity parameter X the relative importance of elastic effects cannot be estimated. This would require the scaling of the relaxation time by a characteristic (overturn) time of the flow, which is, however, not known in advance. The Deborah numbers De given in this study are therefore calculated a posteriori. The Deborah number is here defined as De = max X 11d 11,where the maximum is taken over the whole convection cell. Another important result of the calculation is the Nusselt number Nu, which gives the mean vertical
71 heat flux in the convective system relative to the heat flux in a purely conductive case, i.e. Nu = - (iU’/ay)/( AT/h). Numerically Nu is calculated as in Schubert and Anderson [12]. By the introduction of the non-Newtonian stress contribution s the familiar biharmonic operator is restored on the left-hand side of eqn. (1) so only the right-hand side differs from the Newtonian case. This is a standard method to decouple the equation of motion (1) and the constitutive equations (3) in a finite difference context [5]. Other non-decoupling techniques, suitable for finite element (FE) methods, have been discussed by Marchal and Crochet [13] and Luo and Tanner [ 141. As already pointed out, the problem does not require special boundary conditions for the stress tensor, since at the boundaries simplified versions of (3) derive from the conditions for the stream function 4. As the velocity components (u, w) relate to $ by u = a#/ay and w = - a$~/iZlx, the investigated condition of an enclosed cavity and symmetric boundaries requires 4 = 0 and a2G/a xi = 0 at the boundaries-here x, denotes the direction normal to the boundary. From this it follows that the normal velocity component, vorticity, shear deformation and shear stress vanish at all of the boundaries. Physically these conditions apply to a convection cell surrounded horizontally by identical cells, whereas at the bottom and top of the convecting layer free slip is possible due to the low viscosity of the surrounding layers. These conditions are more appropriate for the application studied than, for example, no-slip conditions. The temperature is set to T = 0 at the upper boundary and T = 1 at the lower boundary. The sidewalls are symmetric, i.e. i3T/dx = 0 there. As an initial condition for the temperature field the conductive profile plus a small sinusoidal perturbation is assumed. In the viscoelastic cases the initial conditions for stress and temperature are those calculated for a lower elasticity parameter. 3. Numerical
method
To solve the system eqns. (l-5) a finite difference approach is used. At the start of the simulation the band matrix resulting from the discretised biharmonic operator in eqn. (1) is factorized once by Gaussian elimination and then stored. This requires only a small amount of the total numerical work of the whole simulation. However, the necessity to store the factorized matrix limits the number of grid points to less than approximately 100 x 100. Because at Rayleigh numbers Ra 2 lo5 thin thermal boundary layers develop, and in cases with Deborah numbers De 2 0.2 also huge stress gradients exist at the margins, it is desirable to increase the resolution there beyond what would be possible with an equidistant mesh. Therefore the
72 possibility is included to use a stretched grid which is successively refined from the centre of the cell towards the margins in geometrical fashion. Harder [15] has shown both by consideration of the local truncation error and by numerical calculations on grids with different resolution that second-order accuracy is preserved in cases with non-equidistant mesh. All terms of the forcing function on the right-hand side (RHS) of eqn. (1) are approximated by central differences. The temperature equation (4) is solved by an alternating direct implicit (ADI) method using a Peaceman-Rachford split. Both diffusion and advection terms are approximated by central differences. If the time development of the convective pattern is of interest, the time step At should be limited by :he Courant criterion, i.e. At I min(AXi/Ui, Ayj/wj) at all grid points (i, j). This criterion is identical to the stability limit of a simple explicit scheme. If only the steady-state solution is desired, however, it is possible with the ADI method to use longer time steps and a steady state is reached with far fewer iterations than would be possible with an explicit scheme. Usually At I n min( Axi/ui, Ayj/wj) is a good choice, where n is the number of grid points in one direction. A further enlarged At is not advisable since the method is less efficient with a too large At and even numerical instability can occur. For a Newtonian fluid only the system of the equation of motion (1) and the temperature equation (4) has to be solved, since the non-Newtonian terms on the RHS of eqn. (1) vanish. The results obtained with the outlined method for a variety of Newtonian test cases are documented in the benchmark study by Blankenbach et al. [9]. From these results it can be estimated that in problems with Ra = lo4 a mesh of 40 x 40 equidistant grid points is sufficient to achieve an accuracy of at least 0.3% in both velocity and temperature fields. This accuracy should be representative for the Newtonian results given in this study. In order to investigate convection in a Maxwellian medium the numerical method described is generalised. The material equations (3) are solved as a decoupled system using two alternative finite difference (FD) techniques already discussed by Crochet et al. [5]. (a) The upwind method: here the time derivatives in eqn. (5) are approximated by forward difference quotients and the spatial derivatives by first-order upwind quotients, i.e. one-sided difference approximations dependent on the sign of the local velocity components. Introducing the notation 7 = 7x,, v = rYv- 7,,, p = 7Yv+ 7,, and labelling the time level by a superscript, while subscripts refer to grid points, the following scheme arises for the shear stress 7:
13 where ~1 =
$At( ) u( + U)/( Xi - Xi-l)>
~3= iAt( I I+ I +W>/(.Yj-Yj-,>,
+_=tAt( c4 =
IuI-~u)/(~i+~ -xi),
I w I - w)/(Yj+I -Vj)-
iAt(
For the other stress components similiar schemes apply. Since only one-sided differences are used, this scheme is only a first-order approximation of the material equations (3). It is therefore desirable to use a method of higher accuracy. The second technique will fulfil this demand. (b) The DuFort-Frankel method: this scheme was introduced by Townsend [16] to solve viscoelastic material equations. Time as well as spatial derivatives are approximated by central differences. Therefore the discretisation extends over three time levels and has the following form:
=-
cl(
Ti”+l,j
-
'l?l,j
) - c2(Tiyj+l
Tiyj_l)-
-
(1- A) q’
+ xov:,-’
+ Xd,,ap~l~l + 4qdxy,
~:j+’ + 2Xd,,ap~~’
4hw$+l = -‘I(
‘in+l,j
-
‘El,,
- 2hd,,ap~,~1 -4Xd,,a$+’ =
(74
)
-
- 2Xd xxav”-l I,
-
“;j_l)
-
4xLVi;-1 -
- 8qd,,,
+ 2hd,,av~f1
-cl(P~+l,j-P~-l,j
c2(v;j+l
1
x at
1
n-1
“ij
(W
+
)-~2(/~r,j+l
-/.~:,j-l)+
4Xd,,a~~-’ (74
where c1 = 2hU/(Xi+l -x~_~) and c2 = 2hw/(~~+~ -Y~_~). To calculate the new stress components (7, v, p)“+l a (3 x 3) matrix has to be inverted for each grid point. In case of a co-rotational model, i.e. a = 0, the system reduces to a set of (2 X 2) matrices, since for the stress trace p = 0 holds within the whole cell. This algorithm is to my knowledge the only FD technique applied to the solution of Maxwellian material equations which posseses second-order accuracy. Although it can be expected that this method is superior to the first-order upwind method described before, a direct comparison is not available yet. This will be done for an example in the next section.
74 For both methods no rigorous stability analysis is available. If a fixed velocity field is considered, the upwind method requires a time step At < A, whereas the DuFort-Frankel scheme allows longer time steps. If the coupled system eqns. (l-4) are solved, however, only a much smaller time step is possible, At -+c h, provided the Deborah number is high enough to induce a substantial coupling of the equation of motion (1) and the material equations (3). In the numerical calculations a time step At in the range of X/20 to h/5 is used. The iteration scheme for the solution of the non-linear system eqns. (l-4) is of a straightforward Picard type, where from the new stress components (7, v)“+’ the non-Newtonian terms on the RHS of eqn. (1) are calculated by the stress decomposition eqn. (2). A modification, similiar to that of Townsend [16], which consists of a sub-iteration for stresses and stream function, i.e. the calculation of the new stress components (7, v, p)‘+’ is repeated with updated coefficients in eqns. (2) and (3), did not improve the stability of the iteration in the problem considered here. Therefore in general the simpler and computationally inexpensive Picard scheme is preferred. If a time-dependent solution is to be calculated, the more restrictive timestep limit arising from the temperature eqn. (4) and material equations (3) must be used for the solution of the coupled system. This is here the condition At c A, following from the explicit scheme for the material equations. For a typical problem considered in this study, around lo4 iterations are needed with such a method. If only the steady-state solution is of interest, the number of iterations can be reduced by using different ‘pseudo-timesteps’ in both equations. With the choice At, = n min( Ax;/ u,, A yj/wj) in the temperature equation and At,s h in the material equations a steady-state solution could usually be calculated with less than lo2 iterations-an obvious improvement compared to a ‘ truely’ time-marching method. Since no analytical, and only very few published numerical, solutions on viscoelastic convection problems are available, it is not easy to verify the numerical scheme quantitatively. Although a qualitative comparison with the results of Ivins et al. [2], given as iso-line plots, was possible in a case of weak coupling, additional verification is desirable. Therefore the following relations are useful as internal consistence checks. Parmentier et al. [17] have shown that for a convection cell, regardless of the rheology applied, the identity Raj, wT d f = lF r :d df holds, where the integration is taken over the whole volume F of the cell. With the Maxwellian material equations (4) it is easy to show that assuming a steady-state solution this identity can be extended to the integrals over the invariants p = tr{ 7) and J = - det{ T} of the stress tensor, i.e.:
75 Combined with a convergence test on refined grids, a check of these criteria is not only a safeguard against a systematic error in the numerical code, but also allows a rough estimate of the accuracy of the calculated solution. 4. Results and discussion 4.1 Convergence tests and stability of the numerical scheme In order to examine the accuracy of the numerical code and to compare the different methods of solving the material equations, a set of convergence tests has been completed. As test problem a case with the low Rayleigh number Ra = 9487 was chosen. This value exceeds the critical Rayleigh number Ra, at which convection sets in by only a factor 12. Such a low value is chosen primarily to avoid the additional complication of thin boundary layers in the temperature field developing at high Rayleigh numbers, Ra > lOORa,. The aspect ratio of the cell is one and the upper convected Maxwell rheology (a = 1) is assumed, since this choice of the parameter a in the Oldroyd derivative (5) is most frequently studied in computational fluid dynamics. In Fig. 1 the relative change in the maximum of the stream function #,,, in the viscoelastic flow field compared to its Newtonian counterpart is 16
I
1
8
I
L
14
-
12
-
10
-
8
-
6
-
4
-
2
-
0
’ 0.0
I
1 + : upwind, 20x20 l : upwind,kOx40 o : DuFort-Frankel, x : DuFort-Frankel,
I
0.1
20x20 40x40
I
I
0.3
0.2
I
0.4
0.5
De
Fig. 1. Relative change of the maximum of the stream function Newtonian result for various Deborah numbers De, grid resolution scheme in the material equations.
$,, compared to the and different numerical
X
X
Fig. 2. Stress trace p at the left part of the lower boundary (y = 0 and x < 0.4). For x > 0.4 the stress profile is essentially flat. The results are obtained by the upwind scheme and the DuFort-Frankel scheme (dotted profiles) at different grid resolution. (Ra = 9487, De = 0.34, h = 0.0015, a = 1).
displayed with increasing Deborah number. As the grid resolution is increased, the solutions obtained with different numerical methods obviously converge against each other, as expected. The results of Fig. 1 show further that with the DuFort-Frankel scheme higher accuracy is obtained than with the simple upwind scheme. This, of course, follows from the fact that the first method is a second-order approximation whereas the upwind scheme is only of first order. In the problem considered in Fig. 2 the maximum misfit in the integral relations (8) is, with the DuFort-Frankel scheme, 1.1% on a 20 x 20 grid and decreases to 0.3% on a 40 x 40 grid. The corresponding values with the upwind scheme are 5.6% and 3.0%. These results are consistent with the approximation error of both schemes and confirm the methods used. The upwind method, however, has the inherent deficiency of introducing an artificial numerical diffusivity into the material equations, which is responsible for smoothed stress gradients at the boundaries. This is obvious from Fig. 2, where the stress trace p at the lower boundary in the vicinity of the convergent stagnation point is shown, calculated with both schemes on two identical grids. Nevertheless, the numerical diffusivity of the upwind scheme allows solutions to be obtained at considerably higher Deborah numbers. At De > 0.4 the upwind method gives an acceptable solution whereas the
77 2400
stress trzca i
a)
I
Fig. 3. As Fig. 2, but De = 0.42, X = 0.0020.
DuFort-Frankel scheme fails at values De = 0.4. The breakdown of the scheme is manifest by artificial oscillations with grid periodicity in the stress fields (Fig. 3), which increase at higher De and finally the iteration diverges. This failure is restricted to the solution of the coupled system eqns. (l-4) and accompanied by a loss of positive definiteness of the stress tensor. Increasing the grid resolution does not change the situation much as demonstrated in Fig. 3. Even with an additional grid refinement to 80 x 80 points or stretched grids, which increase the resolution near the boundaries further by a factor 5, the problem changes only insignificantly. Similiar behaviour is known for most numerical methods in a wide class of viscoelastic flow problems, but the special case of Maxwellian rheology seems the most difficult to handle [4,5]. Only recently finite element techniques have been introduced, which incorporate a one-sided approximation of the advective terms in the material equations and which-combined with a subdivision of the finite elements into finer sub-elements for approximating the stresses-allow a solution at substantially higher De. It is therefore remarkable that even the established FD upwind improves the range of De where well behaved solutions can be obtained. The discretisation of the convective terms in the material equations seems therefore to be of crucial importance. The solution of high Deborah number problems will be further discussed in Section 4.5, but first some aspects of the convection problem will be considered in the low Deborah number regime (De c 0.4), where both numerical schemes yield acceptable solutions.
78 4.2 Comparison of upper convected/ co-rotational models In the Oldroyd derivative eqn. (5) of the constitutive eqn. (3) the constant a is considered as a free rheological parameter. This approach is sometimes referred to as a Johnson-Segalman model. It is now easy to show that by a transition from a = 1 (upper convected) to a = - 1 (lower convected) the solution remains the same with the only exception that the stress trace p changes sign. Therefore only the two extreme cases a = 1 and a = 0 (co-rotational) will be discussed. The choice of an upper convected Maxwell model is much more popular than the co-rotational model since the latter case exhibits unrealistic rheometrical behaviour at De > 0.5 in a simple shear flow, i.e. a decrease in stress with increasing shear rate and, associated with that, the stress tensor will no longer be positive definite, see e.g. ref. 4. Nevertheless, the co-rotational case is a useful rheological model for De -C0.5. The effect of a different choice of the parameter a is therefore investigated for an example already considered in Fig. 2 (Ra = 9487). As a reference, Fig. 4 displays the Newtonian solution, whereas Fig. 5 gives the solution in the case of an upper convected Maxwell model and Fig. 6 refers to the corresponding co-rotational case. The elasticity parameter is in both cases identical, h = 0.0015, the Deborah numbers are different, however, since De depends on the flow field.
TemDerature
St ream-Fc
Shear-stress
Stress-dif.
t .
Vorticitv
Fig. 4. Temperature T, stream function #, vorticity w, shear stress T and normal stress difference Y for Newtonian rheology and Ra = 9487. Isolines: T(1/7,. . . ,6/7); $( - 3, 6,. . . , - 15); ~(30, 60,. . . ,150); T( - 90, - 60,. . . ,60); V( - 400,0,400, 800).
79 Temperature
Stream-Fct.
Vorticitv
Shear-stress
Stress-dif.
Stress-trace
k2!
Fig. 5. As Fig. 4, but with upper convected Maxwell rheology Additional isolines for the stress trace: ~(75,275,. . ,1075).
For both models characteristic exist: whereas in the co-rotational than its Newtonian counterpart,
(a = 1, A = 0.0015, De = 0.34).
differences to the Newtonian solution case the convective flow is more vigorous the upper convected model yields the
a a+ Temperature
Stream-Fct.
Vorticity
0 0
Shear-stress
Stress-dif.
-H_
a
0
l
Fig. 6. As Fig. 4, but with co-rotational
Maxwell
rheology
(a = 0, h = 0.0015, De = 0.52).
80 opposite effect. Although the differences are hardly visible in the isoline plots of temperature and stream function, the relative deviation from the Newtonian values of the Nusselt number and stream function maximum is (- 2.7, - 5.3)% for a = 1 and (+ 7.7, + 8.8)% for a = 0.These values, as the isoline plots, are calculated on the fine 40 X 40 grid using the DuFortFrankel scheme; high accuracy can thus be expected. The fact that both Maxwell models considered produce opposite deviations from the Newtonian flow, is plausible. A co-rotational Maxwell model is characterised in a simple shear flow by an effective viscosity decreasing with shear rate, whereas in a pure shear flow the effective viscosity is constant and equals the Newtonian values. The upper convected case is different: in a simple shear flow the effective viscosity is constant and equals the Newtonian value whereas in a pure shear flow the effective viscosity increases with strain rate. Therefore it seems reasonable that in a ‘real’ flow, which consists of both pure and simple shear, the co-rotational model will exhibit features of a decreased effective viscosity and the upper convected model similiarly an increased effective viscosity. Correspondingly, the changes in the vorticity field are completely different for both Maxwell models. With a = 1 the deviations from the Newtonian solution are minor at this De-only the maxima are slightly reduced, which is obviously due to the less vigorous flow. In the co-rotational case on the other hand, zones of high vorticity gradients develop at the convergent margins, which are not present in the Newtonian flow and which become larger with increasing Deborah number. In fact some small wiggles in the vorticity have already developed at De = 0.52. This seems to be a consequence of the fact that the ‘critical’ Deborah number De = 0.5 for a simple shear flow is slightly exceeded. Since this is a deficiency of the particular rheological model, both numerical methods fail at higher Deborah numbers. The differences in the stress fields of both Maxwell models are mainly restricted to the stress trace p, which vanish identically for a = 0, as in the Newtonian case. For a = 1,however, in p pronounced boundary layers form in the vicinity of convergent zones of neighbouring convection cells. The normal stress difference Y and the shear stress 7, on the other hand, do not differ much for both Maxwell models. Of course characteristic differences exist if the stress fields are compared with their Newtonian counterparts. A comparison of the shear stresses in Figs. 5 and 6 with Fig. 4 illustrates the fact that in a Maxwellian fluid stresses are advected by the flow. For the same reason large stress gradients build up in the vicinity of convergent margins, whereas the gradients are smoothed at divergent margins. The differences described are accompanied by corresponding effects on the quantities measurable at the upper surface (y = 1) of the convection cell
81 10 Heatflux
’
0
0.0
a)
I 0.2
0.4
0.6
0.8
1.0
X 80 Surface velocity
bl
70 60 -
0.0
0.2
0.4
0.6
0.8
1.0
X
Fig. 7. (a) Heat flux and (b) horizontal velocity at the upper surface (y =l) for x = 0.0015 (a = 0, De = 0.52 and a =l, De = 0.35). The dotted profiles give the Newtonian results.
-the most important are the surface heat flux and the horizontal surface velocity. As shown in Fig. 7 these quantities are mainly controlled by the vigour of the flow. With a = 0 the heat flux at the surface is enlarged, whereas with a = 1 it is decreased. The same is true for the surface velocity. Summarising the result, we can say that with Deborah numbers De < 0.5 the impact of the viscoelasticity is small and the changes of the Maxwellian flow compared to its Newtonian equivalent are mainly controlled by the choice of the parameter a in the Oldroyd derivative eqn. (5).
82 4.3 Flow at high Rayleigh numbers From various numerical simulations, e.g. McKenzie et al. [18], it is well known that convection at high Rayleigh numbers differs significantly from the flow at slightly over-critical Rayleigh numbers. The main new features are narrow boundary layers in the temperature and high vorticity gradients at the cell boundaries, In order to investigate the influence of these new features for a Maxwell medium, models were calculated for Rayleigh numbers Ra = (12, 40, 120) X Ra,. This range includes the three characteristic domains of slightly overcritical, intermediate and boundary layer flow [19]. The results of Fig. 8 demonstrate that the changes between convection in a Maxwellian and a Newtonian medium are mainly independent of the Rayleigh number, at least at this low Deborah number. Similiar conclusions were reached by Ivins et al. [2]. However, in their models the Rayleigh number was less than Ra < 30Ra, and the Deborah numbers were so low in the models with higher Rayleigh number, that only a change of less than 0.1% in the maximum stream function can be observed. However, in addition, in Fig. 8 a small tendency at De > 0.3 is visible to reduce the variation of $,, with increasing Rayleigh number.
12
10
I
I
I
I
1
-AlCIr,m[%,l
+
-
a -
*
+ : ua487
.
* : Ra=?1623
X
De Fig. 8. Relative number De.
changes
of $,,
dependent
on the Rayleigh
number
Ra and
Deborah
83 4.4 Time-dependent solutions Up to this point only stationary solutions were considered. Since it is possible that the stationarity of the solutions was artificially induced by the numerical technique, for one example the flow has been calculated with a time-marching algorithm. The time step is At = 0.1X and the initial condition for both temperature and stress fields is the Newtonian solution. The time evolution given in Fig. 9 shows that the calculation converges towards the stationary solution independent of the numerical time step. From Fig. 9 it is obvious that for convection in a Maxwellian medium two characteristic times are relevant: first the relaxation time X and second the transit time t,, which is defined as the time needed for a fluid particle to cross a cell with height h at a mean velocity (u), i.e. t, = (u)/h. The width of the initial drop scales with A, whereas the long-term evolution is controlled by the transit time t,. This typical behaviour is unchanged by varying the time step, grid resolution or approximation of the stress advection terms. The superimposed short scaled variations at times t -c 0.01, however, are highly dependent on these features, hence it must be
I
17.4 /-
I
F-Y-
1
ii 16.8
1
16.6 P,
-1 De=0.34
16.4 tL
16.2 t\
15.8
1
.JL
......,..... ..... ,.
/
"A
0.00
' tll
I
I 0.02
-I
I 0.04
t 0.06
I 0.08
I 0.10
0.12
time Fig. 9. Time evolution of #,, in the case Ra = 9487, De = 0.34, X = 0.0015. The initial condition is the Newtonian solution. The steady-state solution is indicated by the horizontal line. The bars give the relaxation time X and the transit time t,, defined in the text.
84 suspected nique.
that this aspect is an artefact
introduced
by the numerical
tech-
4.5 Convection at high Deborah numbers In the foregoing discussion the range of Deborah numbers has been restricted to De < 0.5. For an upper convected Maxwell model it is possible to extend the simulation to substantially higher values of De provided the upwind approximation is used in the material equations. In Fig. 10 the steady-state solution is given for a similiar case already discussed in Fig. 5, but now with De = 1.06 and X = 0.009. As in the case with low Deborah number (De = 0.34, X = 0.0015, Fig. 5) the convective flow is less vigorous than in the Newtonian case, but the effect is much more pronounced here. For example the Nusselt number Nu, the stream function maximum #,, and the mean upper surface velocity (us) are reduced to (69, 59, 39)s of the corresponding Newtonian solution. The flow is thus effectively damped at higher De. This of course reduces the strain rates within the cell, which explains why De is increasing more slowly
Temperature
Stream-Fct.
Vorticity
+ + + 1 [ Lilizl -B_ +
:i I3
Shear-stress
Stress-dif.
(1
(J
Stress-trace
a
u
u
-
0
+
a
Fig. 10. As Fig. 4, but with upper convected Maxwell rheology (A = 0.0090, De = 1.065, a = 1). Isolines: T(1/7,. . . ,6/7); $( - 1.25, - 2.50,. . . , - 8.75); w( - 80, - 40, 0,. . . ,120); T( - 80, - 40, 0,. . . ) 120); V(- 8000, - 4000, - 2000, - 1000, 0,. . . ) + 8000); ~(125,250, 500 ,..., 16000).
85 than the elasticity parameter X. As the stream lines in Fig. 10 show, the damping effect is most pronounced in the vicinity of the cell boundaries. Therefore the mean surface velocity is more effected by the increased elasticity of the flow than the other cited quantities. A further new feature of the solution is the broad zone of negative vorticity around the boundary. In isoviscous Newtonian convection such a zone exists only in cases with no-slip boundary conditions. Similarly there is a major change in the stress fields. At first glance it seems that the effects visible at lower Deborah number are only enhanced: due to stress advection the shear stress pattern is now rotated by almost 90” relative to the Newtonian solution and the maxima of the normal stress patterns are squeezed towards the boundary. However, although a refined mesh of 80 x 80 non-equidistant grid points is used, it is impossible to resolve the normal stress pattern near the boundary, where a singularity seems to exist. This is a fundamentally new feature at higher De. For low De grid refinement yields convergent solutions all over the cell, including the boundaries (Fig. 2). At the higher Deborah number De = 1.065 this is no longer true. If the grid is refined, the normal stresses diverge along the boundaries starting from the stagnation point at the cell corner (Figs. 11(a) and (c)). Since the numerical values at the boundary are approximately inversely proportional to the mesh size, it seems that an artificial numerical diffusivity in the material equationsas introduced by the first-order upwind approximation-is needed to cope with such singularity. This would explain, why the DuFort-Frankel scheme fails at higher De. In spite of the singularity at the boundary the solution converges in the interior of the cell. As an example this is demonstrated in Figs. 11(c) and (d), where the stress component p calculated at various grids is shown just one grid point (on the coarsest grid) apart from the boundary. Here grid refinement yields a convergent solution. Thus convergence in the interior of the cell is not affected by the singularity at the boundary. Since all the models presented in this chapter are calculated by the first-order upwind scheme, only limited accuracy can be expected. In fact the maximal numerical error in the integral identities (8) is (26, 15, 9.6)s using the three stretched grids of Figs. 11(b) and (d). Although an error in the range of 10% is still considerable for an 80 X 80 mesh, which includes grid refinement towards the boundaries, the main qualitative features of the solutions are not affected by the mesh size. The solution displayed in Fig. 10 seems, therefore, to be a reliable approximation to the Maxwellian flow problem considered. A further but limited, increase in the Deborah number is possible, although, the misfit in the integral identities eqn. (8) then grows rapidly. Using the finest grid, the upwind method finally fails, if the elasticity parameter h (which is the input parameter) is further increased by
I
x*x
r,
I ni
l
,
x
%"
)*"r
.x
,
1200
%"I,
stress trace
800
%x
Il.
4OxkO
. : grid
:
1600
: . +.
*. .
..
80x80
,
grid
x :
X.
* : pr1.i20: 20
110x10
Fig. 11. (a) Stress trace JJ at the left boundary trace p close to the left boundary at x = l/20 parameters are identical to Fig. 9.
0.0 0
0.2 -
0.4 -
100
o’nI
l.O(
0.)
Y.
r
0.6
grid
.: .
.
' x
5oa
*
x
x
x
x
10000
15Oa
20000 25000
. : prld 20x20
. : prld 40x40
x : grl6 BOX80
I---* *
.. x . xX
30000
for three equidistant grids; (b) as (a), but grids are refined towards the boundaries; (c) stress (grids are equidistant); (d) as (c), but at x = 0.011 and non-equidistant grids. All other model
2000
c)
0
0.6
0.8
1.0
87 a factor of 2. In this case similiar problems arise as with the DuFort-Frankel scheme at lower De. It is tempting to speculate how the performance of the algorithm could be improved for high-De problems. Obviously the presence of a normal stress singularity starting at a stagnation point and evolving in the upstream direction is of crucial importance. But in cases as considered here, where in one cell only one convection roll exists, no internal stagnation points are present and singularities evolve only at the cell boundaries starting from the corner points. This could be avoided by a discretisation suggested by Gatski and Lumley [20], where the normal stresses are approximated at midpoints of grid cells, i.e. (v, P)i+i/2,j+i/2> and the velocity components at the centres of grid cell boundaries, i.e. ~~+i,~, j, IJ_ i,*,, and w equivalent. The advantage is that there is no grid point for the normal stresses where all velocity components vanish. Therefore no stagnation point would coincide with a grid point and a singular normal stress should be avoided. Although such an approximation may not be favourable in general, it appears very suitable for a convection problem of an enclosed cavity with symmetric boundaries. For such a problem the method of Gatski and Lumley [20] should be reconsidered and compared with the more straightforward approach used here. Nonetheless, convergence for an unlimited range of De cannot be expected with such a modified method. It should be remembered, however, that even with the simplest method used here Maxwellian convection problems can be solved with a flow which differs very significantly from its Newtonian counterpart. 5. Concluding remarks This study has demonstrated the major differences between convection with Maxwellian versus Newtonian rheology and shown that thermal convection is a very suitable test case for numerical methods simulating viscoelastic flow. It has been possible to extend the simulation up to Deborah numbers De = 1.0, which is sufficient to induce significant changes in the flow fields. A main new feature at high De is the presence of a normal stress singularity along the boundary, which is absent in Newtonian or in low-De flow. Since this singularity is starting from the stagnation points at the cell corners, this behaviour is presumably related to the well known singularity of the upper convected Maxwell model in a pure shear flow [4,5]. It is common experience [4,7,13] that rheological models with more realistic response are numerically easier to handle. Examples are the Jeffreys model and its non-linear generalisations, i.e. the Giesekus and Leonov models. Certainly the numerical approach used can be extended to such models, but for geophysical applications these models are up to now only of
88
limited use, since it would be difficult to relate the additional model parameters to the Earth. In further studies [3,15] it is thus preferred to investigate first the geophysically better understood effects of material parameters which depend on temperature and hydrostatic pressure. This approach can also remove the singularities arising from stagnation points and the numerical simulation is more stable than in the cases studied here. Further studied may show that the Jeffreys model and related ones are relevant to the Earth. From a fluid dynamical point of view, thermal convection with such rheological models is certainly of interest. Acknowledgements I wish to thank W.R. Jacoby and U. Christensen for suggesting this problem to me and for many helpful discussions during this work. Financial support by the Deutsche Forschungsgemeinschaft is gratefully acknowledged. References 1 I.A. Eltayeb, Proc. R. Sot. London, Ser., A, 356 (1977) 161-176. 2 E.R. Ivins, T.W.J. Unti and R.J. Phillips, Geophys. Astrophys. Fluid Dyn., 22 (1982) 103-132. 3 H. Harder, U. Christensen and W. Jacoby, in preparation. 4 J. van der Zanden and M. Hulsen, J. Non-Newtonian Fluid Mech., 29 (1988) 93-117. 5 M.J. Crochet, A.R. Davies and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, 1984. 6 D.D. Joseph, M. Renardy and J.C. Saut, Arch. Ration. Mech. Anal., 87 (1985) 213-251. 7 B. Debbaut, J.M. Marchal and M.J. Crochet, J. Non-Newtonian Fluid Mech., 29 (1988) 119-146. 8 G. de Vahl Davies and I.P. Jones, Int. J. Num. Meth. Fluids, 3 (1983) 227-248. 9 B. Blankenbach, F. Busse, U. Christensen, L. Cserepes, D. Gunkel, U. Hansen, H. Harder, G. Jarvis, M. Koch, G. Marquart, D. Moore, P. Olson, H. Schmeling and T. Schnaubelt, Geophys. J. Int., 98 (1989) 23-38. 10 M.A. Mendelson, P.-W. Jeh, R.A. Brown and R.C. Armstrong, J. Non-Newtonian Fluid Mech., 10 (1982) 31-54. 11 U. Christensen, Geophys. J.R. Astron. Sot., 77 (1984) 343-384. 12 G. Schubert and C.A. Anderson, Geophys. J.R. Astron. Sot., 80 (1985) 575-601. 13 J.M. MarchaI and M.J. Crochet, J. Non-Newtonian Fluid Mech., 26 (1987) 77-114. 14 X.-L. Luo and R.I. Tanner, J. Non-Newtonian Fluid Mech., 31 (1989) 143-162. 15 H. Harder, Dissertation, Johannes Gutenberg-Universitat Mainz, Department of Geosciences, 1988 (in German). 16 P. Townsend, J. Non-Newtonian Fluid Mech., 14 (1984) 265-278. 17 E.M. Parmentier, D.C. Turcotte and K.E. Torrance, J. Geophys. Res., 81(1976) 183991846. 18 D.P. McKenzie, J.M. Roberts and N.O. Weiss, J. Fluid Mech., 62 (1974) 465-538. 19 G.T. Jarvis and W.R. Peltier, Geophys. J.R. Astron. Sot., 68 (1982) 389-427. 20 T.B. Gatski and J.L. Lumley, J. Comput. Phys., 27 (1978) 42-70.