Ilyana Anwar

Norsarahaida Amin

Ioan Pop
email: pop.ioan@yahoo.co.uk

Mixed convection boundary layer flow of a viscoelastic fluid over a horizontal circular cylinder

Keywords: Viscoelastic fluid, Boundary layer, Mixed convection, Circular cylinder, Numerical method

The steady mixed convection boundary layer flow of a viscoelastic fluid over a horizontal circular cylinder in a stream flowing vertically upwards is numerically studied for both cases of heated and cooled cylinders. The governing partial differential equations are transformed into dimensionless forms using an appropriate transformation and then solved numerically using the Keller-box method. The comparison between the solutions obtained and those for a Newtonian fluid is found to be very good.

Effects of the mixed convection and elasticity parameters on the skin friction and heat transfer coefficients for a fluid having the Prandtl number equal to one are also discussed. It is found that for some values of the viscoelastic parameter and some negative values of the mixed convection parameter (opposing flow) the boundary layer separates from the cylinder. Heating the cylinder delays separation and can, if the cylinder is warm enough, suppress the separation completely. Similar to the case of a Newtonian fluid, cooling the cylinder brings the separation point nearer to the lower stagnation point.

However, for a sufficiently cold cylinder there will not be a boundary layer.

Introduction

There are many fluids whose behaviour cannot be described by the classical Navier-Stokes and boundary layer equations. The homogeneous and incompressible second grade fluid is one of the many models that have been proposed to describe the non-Newtonian behaviour of such fluids (see Kumari et al. [START_REF] Kumari | Nonsimilar mixed convection flow of a non-Newtonian fluid past a vertical wedge[END_REF]). Mechanics of non-linear fluids present a special challenge to researchers due to its many A c c e p t e d m a n u s c r i p t 2 practical applications, for example, in the design of thrust bearings and radial diffusers, drag reduction, transpiration cooling and thermal oil recovery, to mention just a few. One of the simplest ways in which the viscoelastic fluids have been classified is the methodology proposed by Rivlin and Ericksen [START_REF] Rivlin | Stress deformation relations for isotropic materials[END_REF] and Truesdell and Noll [START_REF] Truesdell | The Non-Linear Field Theories of Mechanics[END_REF], who presented constitutive relations for the stress tensor T as a function of the symmetric part of the velocity gradient D , and its higher (total) derivatives (see Massoudi [START_REF] Massoudi | Local non-similarity solutions for the flow of a non-Newtonian fluid over a wedge[END_REF]). The boundary layer theory for second-grade fluids has been developed by Oldroyd [START_REF] Oldroyd | On the formulation of rheological equations of state[END_REF], Beard and Walters [START_REF] Beard | Elastico-viscous boundary layer flows[END_REF] and Rajagopal et al. [START_REF] Rajagopal | Flow of a viscoelastic fluid over a stretching sheet[END_REF]. However, the boundary layer equations are an order higher than those for the Newtonian (viscous) fluid and the adherence boundary conditions are insufficient to determine the solution of these equations completely. Therefore, we need a boundary condition in addition to the usual adherence boundary conditions. Rajagopal [START_REF]On boundary conditions for fluids of the differential type[END_REF][START_REF] Sequeira | On the creeping flow of the second grade fluid[END_REF] and Rajagopal and Kaloni [START_REF] Rajagopal | Some remarks on boundary conditions for fluids of the differential type[END_REF] have discussed this question in detail. Garg and Rajagopal [START_REF] Garg | Stagnation point flow of a non-Newtonian fluid[END_REF][START_REF] Garg | Flow of a non-Newtonian fluid past a wedge[END_REF] studied the flow of a fluid of second grade near the stagnation point of a semi-infinite wall by augmenting the boundary condition at infinity. It has been shown that their results agree well with the results of Rajeswari and Rathna [START_REF] Rajeswari | Flow of a particular class of non-Newtonian visco-elastic and viscoinelastic fluids near a stagnation point[END_REF] based on the series expansion for small values of the viscoelastic parameter  , which multiplies the highest-order spatial derivative in their equation. The advantage of augmenting the boundary conditions over the perturbation approach is that the analysis is valid even for large values of the parameter  and, as shown by Garg

and Rajagopal [START_REF] Garg | Stagnation point flow of a non-Newtonian fluid[END_REF], significant deviations from the Newtonian behaviour are possible for even moderately large values of  . The augmenting of the boundary condition at infinity has also been used by Kumari et al. [START_REF] Kumari | Nonsimilar mixed convection flow of a non-Newtonian fluid past a vertical wedge[END_REF] for the steady non-similar mixed convection boundary layer flow of a viscoelastic fluid over a permeable vertical wedge, by Ariel [START_REF] Ariel | On extra boundary condition in the stagnation point flow of a second grade fluid[END_REF] for the two-dimensional stagnation point flow of a second grade fluid, by Garg [START_REF] Garg | Non-Newtonian flow over a wedge with suction[END_REF] for the flow of an incompressible fluid of second grade past a wedge with suction at the surface and by Cortell [START_REF] Cortell | A note on flow and heat transfer of a viscoelastic fluid over a stretching sheet[END_REF] for the flow and heat transfer of a second grade fluid past a stretching sheet. Very recently, Hayat and Sajid [START_REF] Hayat | Analytic solution for axisymmetric flow and heat transfer of a second grade fluid past a stretching sheet[END_REF] presented an analytic solution using the homotopy analysis method (HAM) for the flow and heat transfer of a second grade fluid over a radially stretching sheet. Also, Mushtaq et al. [START_REF] Mushtaq | Mixed convection flow of second grade fluid along a vertical stretching flat surface with variable surface temperature[END_REF] studied the effects of thermal buoyancy on flow of a second grade fluid along a vertical, continuous stretching sheet of which the velocity and temperature distributions are assumed to vary according to a power-law form. The governing partial differential equations are non-similar and they were solved using (1) the series expansion method together with the Shanks transformation, (2) the local non-similarity method with second level of truncation and (3) the Keller-box method for some values of the mixed convection parameter. However, the authors have not mentioned whether they have used the extra boundary condition at infinity or not.

In this paper, the steady non-similar mixed convection boundary layer flow of a viscoelastic fluid over a horizontal circular cylinder is considered. The coupled nonlinear partial differential equations have been compared with those of Merkin [START_REF] Merkin | Mixed convection from a horizontal circular cylinder[END_REF] and Nazar et al. [START_REF] Nazar | Mixed convection boundary-layer flow from a horizontal circular cylinder in micropolar fluids: case of constant wall temperature[END_REF]. The results have also been compared with those of Hiemenz [START_REF] Hiemenz | Die Grenzschicht an einem in den gleichfrmigen Flssigkeitsstrom eingetauchten geraden Kreiszylinder[END_REF] and Eckert [START_REF] Eckert | Die Berechnung des Wrmeberganges in der laminaren Grenzschicht umstrmter Krper[END_REF] for a viscous fluid and with those of Ariel [START_REF] Ariel | On extra boundary condition in the stagnation point flow of a second grade fluid[END_REF] for a viscoelastic fluid when the buoyancy forces are absent. It is shown that the agreement between all these results is very good. We wish also to mention to this end that to our best knowledge this classical very important problem has not been studied before for a viscoelastic fluid so that the results are new for these fluids.

Basic equations

The problem that we will study in this paper is the steady mixed convection boundary layer flow past an isothermal horizontal circular cylinder of radius a placed in a viscoelastic fluid. . Under these assumptions along with the Boussinesq approximation, the boundary layer equations can be written as follows, see Garg and

Rajagopal [START_REF] Garg | Flow of a non-Newtonian fluid past a wedge[END_REF] and Mushtaq et al. [START_REF] Mushtaq | Mixed convection flow of second grade fluid along a vertical stretching flat surface with variable surface temperature[END_REF],

Continuity equation 0       y v x u (1) Momentum equation ) / ( sin ) ( 2 2 3 3 2 2 0 2 2 a x T T g y v y u y u v y u u x k y u x d u d u x u v x u u e e                                             (2)
Energy equation

2 2 y T y T v x T u          (3)
subject to the boundary conditions 0 , , 0 / ), ( 0 , 0 at , 0

              x y T T y u x u u x y T T v u e w (4)
where x and y are the Cartesian coordinates measured along the surface of the cylinder starting from the lower stagnation point of the cylinder and y is the coordinate measured normal to the surface of the cylinder, u and v are the velocity components along the 

x and  y axes, ) (x u e is the velocity outside the boundary layer, T is the fluid temperature, g is the acceleration due to gravity 0 k is the viscoelasticity and    , , and  are the thermal diffusivity, thermal expansion coefficient, dynamic viscosity and density of the viscoelastic fluid.

We introduce now the following non-dimensional variables [START_REF] Oldroyd | On the formulation of rheological equations of state[END_REF] where Re = U ∞ a/v is the Reynolds number. Substituting (5) into Eqs.

                 T T T T U x u x u U v v U u u a y y a x x w e e  , ) ( ) ( ) ( Re , ), ( Re , 2 / 1 2 / 1
                                       (7) 2 2 Pr 1 y y v x T u           (8)
subject to the boundary conditions (4), which become 0 , as 0 , 0 ), (

0 , 0 at 1 , 0              x y y u x u u x y v u e   ( 9
)
where Pr is the Prandtl number, K is the viscoelastic parameter and  is the mixed convection parameter, which are defined as

2 0 Re Gr λ ,     a U k K (10) with 2 3 / ) (   a T T g Gr w   
being the Grashof number. It should be mentioned that 0   corresponds to assisting flow (heated cylinder), 0   corresponds to opposing flow (cooled cylinder) and 0   corresponds to the forced convection flow, respectively.

Solution

To solve Eqs. ( 6)-( 8) with the boundary conditions ( 9), we follow Merkin [START_REF] Merkin | Mixed convection from a horizontal circular cylinder[END_REF] and assume that x x u e sin ) (  . Thus, we look for a solution of these equations of the following form

) , ( ), , ( y x y x f x      ( 11
)
where  is the stream function defined as

y u    /  and x v     / 
. Substituting [START_REF] Garg | Stagnation point flow of a non-Newtonian fluid[END_REF] into Eqs. ( 7) and ( 8), we obtain

                                                                                               2 2 2
y f x f y x f y f x y x f y f y x f y f y f x f y f y x f x y f y f f y f y f K x x x x x y f y f f y f   (12)
A c c e p t e d m a n u s c r i p t

6                        y x f x y f x y f y     2 2 Pr 1 ( 13
)
subject to the boundary conditions 0 , as 0 , 0 , sin

0 , 0 at 1 , 0 2 2                  x y y f x x y f x y y f f   (14)
It is worth mentioning that when 0  K , Eqs. ( 12) and ( 13) reduce to the equations governing the mixed convection boundary layer flow of a viscous and incompressible (Newtonian)) fluid studied by Merkin [START_REF] Merkin | Mixed convection from a horizontal circular cylinder[END_REF].

At the lower stagnation point of the cylinder, 0  x , Eqs. ( 12) and ( 13) reduce to the following ordinary differential equations

  0 ' ' ' ' ' ' 2 1 ' ' ' ' ' ' 2 2         f f f f f K f f f f iv   (15) 0 ' ' ' Pr 1     f ( 16
)
with the boundary conditions

0 ) ( , 0 ) ( ' ' , 1 ) ( ' 1 ) 0 ( , 0 ) 0 ( ' ) 0 (            f f f f ( 17
)
where primes denote the differentiation with respect to y .

The physical quantities of principal interest in this problem are the skin friction coefficient f C and heat transfer coefficient w Q . We define these coefficients in non-dimensional form as

) ( Re , Re 2 / 1 2 2 / 1       T T k q a Q U C w w w w f   ( 18
)
where w  and w q are the skin friction and heat flux from the surface of the cylinder, which are given by

0 0 2 2 2 0 0 2                                               y w y y w y T k q y u x u y u v y x u u k y u   ( 19
)
k being the thermal conductivity of the viscoelastic fluid. Using ( 5) and [START_REF] Garg | Stagnation point flow of a non-Newtonian fluid[END_REF], we obtain A c c e p t e d m a n u s c r i p t

7 0 0 2 2 ,                          y w y f y Q y f x C  (20)
For small values of the viscoelastic parameter ) 1 ( K , we look for a solution of Eqs. ( 15) and ( 16)

subject to the boundary conditions [START_REF] Hayat | Analytic solution for axisymmetric flow and heat transfer of a second grade fluid past a stretching sheet[END_REF] in series of the form

... ) ( ) ( ) ( ) ( ... ) ( ) ( ) ( ) ( 2 2 1 0 2 2 1 0         y K y K y y y f K y f K y f y f     ( 21
)
where the functions

) ( y f i and ) ( y i  ( 2 , 1 , 0  i
) are given by the following three sets of equations with the corresponding boundary conditions

; 0 ) ( , 1 ) ( ' 1 ) 0 ( , 0 ) 0 ( ' ) 0 ( 0 ' ' ' Pr 1 0 1 ' ' ' ' ' ' 0 0 0 0 0 0 0 0 0 2 0 0 0 0                     f f f f f f f f (22) ; 0 ) ( , 0 ) ( ' 0 ) 0 ( , 0 ) 0 ( ' ) 0 ( ' ' ' ' Pr 1 ' ' ' ' ' ' 2 ' ' ' ' ' ' ' ' 2 ' ' ' ' ' 1 1 1 1 1 0 1 1 0 1 2 0 0 0 0 0 1 1 0 1 0 1 0 1                         f f f f f f f f f f f f f f f f f (23) 0 ) ( , 0 ) ( ' 0 ) 0 ( , 0 ) 0 ( ' ) 0 ( ' ' ' ' ' Pr 1 ' ' ' ' 2 ) ' ' ' ' ' ' ' ' ( 2 ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' 2 ' ' ' ' ' 2 2 2 2 2 0 2 1 1 2 0 2 1 0 0 1 1 0 0 1 1 0 1 1 2 1 2 2 0 2 0 2 0 2                               f f f f f f f f f f f f f f f f f f f f f f f f f f (24)
Thus, the reduced skin friction and heat transfer from the surface of the cylinder are given by .

.. ) 0 ( ' ) 0 ( ' ) 0 ( ' ) 0 ( ' ... ) 0 ( ' ' ) 0 ( ' ' ) 0 ( ' ' ) 0 ( ' ' 2 2 1 0 2 2 1 0          K K K f K f f f     (25) for 1  K .
A c c e p t e d m a n u s c r i p t 8

Results and discussion

The systems of equations [START_REF] Garg | Flow of a non-Newtonian fluid past a wedge[END_REF][START_REF] Rajeswari | Flow of a particular class of non-Newtonian visco-elastic and viscoinelastic fluids near a stagnation point[END_REF][START_REF] Ariel | On extra boundary condition in the stagnation point flow of a second grade fluid[END_REF] and [START_REF] Hiemenz | Die Grenzschicht an einem in den gleichfrmigen Flssigkeitsstrom eingetauchten geraden Kreiszylinder[END_REF][START_REF] Eckert | Die Berechnung des Wrmeberganges in der laminaren Grenzschicht umstrmter Krper[END_REF](24) were solved numerically for some values of the mixed convection parameter  and viscoelastic parameter K using the Keller-box method, which is described in the book by Cebeci and Bradshaw [START_REF] Cebeci | Physical and Computational Aspects of Convective Heat Transfer[END_REF]. All three cases of the assisting ( 0   ) flow, opposing C and w Q are also given in Tables 1 to 4. It is found that both skin friction f C and heat transfer w Q coefficients decrease as K is increased. A similar trend has been observed by Gard and Rajagopal [START_REF] Garg | Flow of a non-Newtonian fluid past a wedge[END_REF] for the problem of forced convection flow of a viscoelastic fluid past a wedge and by Kumari et al. [START_REF] Kumari | Nonsimilar mixed convection flow of a non-Newtonian fluid past a vertical wedge[END_REF] for the problem of mixed convection flow of a viscoelastic fluid past a vertical wedge. This can be attributed to the thickening of momentum and thermal boundary layers as K increases. Kumari et al. [START_REF] Kumari | Nonsimilar mixed convection flow of a non-Newtonian fluid past a vertical wedge[END_REF] have explained that the increase in the boundary layer thickness with K can be attributed to tensile stress in the boundary layer, which cause an axial contraction and hence the thickening of the boundary layer in the transverse direction. We can also see from Figs. 4 to 7 and tables 1 to 4 that, as it is expected, the boundary-layer separates from the cylinder for some negative values of  (cooling cylinder) and also for some positive values of  (heated cylinder). On the other hand, the results show that, as for the case of a Newtonian fluid, increasing  delays separation of the 12), see Merkin [START_REF] Merkin | Mixed convection from a horizontal circular cylinder[END_REF], . The values obtained by direct numerical integration of Eqs. ( 15) and ( 16) using the boundary conditions [START_REF] Hayat | Analytic solution for axisymmetric flow and heat transfer of a second grade fluid past a stretching sheet[END_REF] with the extra boundary condition

0 sin ) cos ( 0 3 3               x x y f x y s  (26) Though 0 ) / ( 0 2 2     y y f at s x x  , the streamwise velocity component 0 /    y f near 0  y and so   0 / 0 3 3     y y f at s x x  . From (
0 ) ( ' '   f
are also included in Table 5. The classical results reported by Hiemenz [START_REF] Hiemenz | Die Grenzschicht an einem in den gleichfrmigen Flssigkeitsstrom eingetauchten geraden Kreiszylinder[END_REF] and Eckert [START_REF] Eckert | Die Berechnung des Wrmeberganges in der laminaren Grenzschicht umstrmter Krper[END_REF] A . In addition, the results reported by Ariel [START_REF] Ariel | On extra boundary condition in the stagnation point flow of a second grade fluid[END_REF] for a viscoelastic fluid ( 0  K

) obtained using the direct numerical solution of Eqs. ( 15) and ( 16) with the boundary conditions ( 17) when 0   (forced convection flow) and the series solutions for small values of K ( ) 1  are included in Table 5. In all cases he results are found to be in very good agreement. Further, it is noticed that the results based on the series expansion agree well with the exact numerical solution of Eqs. ( 15) and [START_REF] Cortell | A note on flow and heat transfer of a viscoelastic fluid over a stretching sheet[END_REF]. Therefore, the advantage of using the extra boundary condition

0 ) ( ' '   f
over the series expansion is that the analysis is valid even for large values of the parameter K . Garg and Rajagopal [START_REF] Garg | Stagnation point flow of a non-Newtonian fluid[END_REF] have shown that significant deviations from the Newtonian behaviour are possible for even moderately large values of K .

Finally, Figs. 9 and 10 illustrate the velocity and temperature profiles at the lower stagnation point of the cylinder against y for several values of the parameter K when 1   (assisting flow) and 1    (opposing flow), respectively, and 1 Pr  . These figures show how the viscoelastic parameter K affects the fluid velocity and temperature profiles. Thus, Figs. 9 show that the velocity profiles decrease when K is increased and that the values of these profiles are lower for a viscoelastic fluid than for a Newtonian fluid ( 0  K

). Therefore, the thickness of the velocity boundary layer for a viscoelastic fluid is higher than for a Newtonian fluid. Further, it is seen from Figs. 7 that the velocity profiles are lower for the case of the assisting flow ( 1  

) than those for an opposing flow (

1   
), respectively. The reverse trend is observed for the temperature profiles, which are shown in Fig. 10.

Conclusion

The 12) and ( 13) subject to the boundary conditions ( 14) is not possible. The results of Merkin [START_REF] Merkin | Mixed convection from a horizontal circular cylinder[END_REF] for viscous fluids can be recovered easily when the viscoelastic parameter 0  K .

Aknowledgements

The authors wish to express their thanks to the reviewers for the valuable comments and suggestions.

. λ = 0 λ = 1 λ = -1 Series Numerical Series Numerical Series Numerical K ) 0 ( ' ' f ) 0 ( '   ) 0 ( ' ' f ) 0 ( '   ) 0 ( ' ' f ) 0 ( '   ) 0 ( ' ' f ) 0 ( '   ) 0 ( ' ' f ) 0 ( '   ) 0 ( ' ' f ) 0 ( '   0 1.

A c c e p t e d m a n u s c r i p t 3 governing

 3 the flow have been solved numerically using a very efficient finite-difference scheme, known as the Keller-box method, which is described in the book by Cebeci and Bradshaw[START_REF] Cebeci | Physical and Computational Aspects of Convective Heat Transfer[END_REF] after augmenting the boundary condition at infinity. The effects of the mixed convection and viscoelastic parameters on the skin friction and heat transfer around the cylinder are studied. The effects of these parameters on the velocity and temperature profiles near the lower stagnation point of the cylinder are presented here only for the Prandtl number equal to one, although they can be obtained for other values of the Prandtl number, as well. The particular cases of the present results for a viscous (Newtonian) fluid

Fig. 1 .

 1 Fig. 1. Physical model and coordinate system.

  (1) -(3), we get the following non-

  . In order to save space we present here results only the case when the Prandtl number Pr is one. The value of 10 for the computation. The iterations were continued until an accuracy of 10 -6 was achieved.The comparison of the present results for the skin friction coefficient f C and the heat transfer from the cylinder w Q with those of Merkin[START_REF] Merkin | Mixed convection from a horizontal circular cylinder[END_REF] and Nazar et al.[START_REF] Nazar | Mixed convection boundary-layer flow from a horizontal circular cylinder in micropolar fluids: case of constant wall temperature[END_REF] for a Newtonian fluid ( 0  K) is shown in Figs.2 and 3. The comparison shows that the numerical solutions obtained by the present authors are in very good agreement with those of Merkin[START_REF] Merkin | Mixed convection from a horizontal circular cylinder[END_REF] and Nazar et al.[START_REF] Nazar | Mixed convection boundary-layer flow from a horizontal circular cylinder in micropolar fluids: case of constant wall temperature[END_REF]. We are, therefore, confident that the present results are very accurate. The trend exhibited by the curves is consistent with the expected effect of favourable pressure gradient. A positive mixed convection or buoyancy forceparameter ( 0  ) induces a favourable pressure gradient that enhances the fluid motion, which in turn increases the skin friction coefficient f C and hence the local heat transfer coefficient w Q also increases and the heat transfer from the cylinder is increased with  (see Fig.3).The variation off C and w Q with x for a viscoelastic fluid ( 0  K ) is shown in Figs. 4 to 7 at different positions x at the surface of the cylinder and some values of  when 1  . Numerical values of f

A c c e p t e d m a n u s c r i p t 9 boundary

 9 layer from the cylinder and that separation can be suppressed completely in the range on K , below which a boundary layer solution is not possible. As it was explained by Merkin[START_REF] Merkin | Mixed convection from a horizontal circular cylinder[END_REF] for the case of a Newtonian fluid ( cooled and the natural convection boundary layer would start at the top stagnation point of the cylinder (   x) and for sufficiently small  there comes a point where the flow of the stream upwards cannot overcome the tendency of the fluid next to the cylinder to move downwards under the action of the buoyancy forces.Using the boundary conditions (14) and the fact that 0

c c e p t e d m a n u s c r i p t 10

 10 flow) are also included in this Table



  steady mixed convection boundary layer flow of an incompressible viscoelastif fluid past an isothermal horizontal circular cylinder has been investigated numerically. The governing boundary layer equations are transformed into a non-dimensional form and the resulting nonlinear system of partial differential equations is solved numerically using the Keller-box method. Both the cases when the mixed increases with the increase of the viscoelastic parameter K . On the other hand, for a cooler cylinder ( 0  ), the buoyancy forces also retard the fluid and therefore the separation point is brought nearer to the lower stagnation point. A value of ) found for which the boundary layer separates at this point. It is found that this value of ) (K s  decreases as the parameter K increases.

A c c e p t e d m a n u s c r i p t 11 For

 11 solution of Eqs. (

Fig. 2 .

 2 Fig. 2. Comparison of the local skin friction f C for K = 0 (Newtonian fluid), Pr = 1 and various values of λ.

Fig. 3 .

 3 Fig. 3. Comparison of the heat transfer coefficient w Q for K = 0 (Newtonian fluid), Pr = 1 and various values of λ.

A c c e p t e d m a n u s c r i p t 14 Fig. 4 .Fig. 5 .

 1445 Fig. 4. Variation of the local skin friction coefficient C f for K = 0.2, Pr = 1 and various values of λ.

Fig. 6 .

 6 Fig. 6. Variation of the local skin friction coefficient C f for K = 1, Pr = 1 and various values of λ.

Fig. 7 .

 7 Fig. 7. Variation of the local heat transfer coefficient w Q for 1  K , Pr = 1 and various values of λ.

A c c e p t e d m a n u s c r i p t 16

 16

A c c e p t e d m a n u s c r i p t 17 Fig. 8 .Fig. 9 .

 1789 Fig. 8. Variation of the boundary layer separation point x s with  for

A c c e p t e d m a n u s c r i p t 19 Fig. 10 .

 1910 Fig. 10. Temperature profiles at the lower stagnation point for Pr = 1 and various value of K: a) 1   (assisting flow) and b) 1    (opposing flow).

Table 1 .

 1 Values of the local skin friction coefficient C f for K = 0.2, Pr = 1 and various values of λ.

	λ

Table 2 .

 2 Values of the local heat transfer coefficient w Q for K = 0.2, Pr = 1 and various values of λ.

	λ

Table 3 .

 3 Values of the local skin friction coefficient C f for K = 1, Pr = 1 and various values of λ.

	λ

Table 4 .

 4 Values of the local heat transfer coefficient C f for K = 1, Pr = 1 and various values of λ.

	λ


Nabil Kerdid
email: nkerdid@ccis.imamu.edu.sa

Hervé Le Dret
email: ledret@ccr.jussieu.fr

Abdelkader Saïdi
email: saidi@math.u-strasbg.fr

Numerical approximation for a nonlinear membrane problem

Keywords: Nonlinear elasticity, Membranes/thin films, Finite element approximation

We present a numerical study of large deformations of nonlinearly elastic membranes. We consider the nonlinear membrane model obtained by H. Le Dret and A. Raoult using Γ-convergence, in the case of a Saint Venant-Kirchhoff bulk material. We consider conforming P 1 and Q 1 finite element approximations of the membrane problem and use a nonlinear conjugate gradient algorithm to minimize the discrete energy. We present numerical tests including membranes subjected to live pressure loads.

INTRODUCTION

The purpose of this article is to devise numerical approximations of large deformations of a nonlinearly elastic membrane. The nonlinear membrane model used here was obtained in [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF], with refinements in [START_REF] Shafrir | The membrane shell model in nonlinear elasticity: A variational asymptotic derivation[END_REF]. The relevance of this model stems from the fact that it was derived from three-dimensional nonlinear elasticity by means of a rigorous convergence method. Similar nonlinear membrane models had already been obtained previously by Pipkin, directly in the context of standard two-dimensional membrane theory and using relaxation, see [START_REF] Pipkin | The relaxed energy density for isotropic elastic membranes[END_REF].

Our numerical study of the nonlinear membrane model is made possible due to the explicit formula for the nonlinear membrane energy given in [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] in the A c c e p t e d m a n u s c r i p t case of the Saint Venant-Kirchhoff bulk material. For a general bulk material, an explicit computation of the corresponding nonlinear membrane energy entails the determination of the quasiconvex envelope of a function defined on the space of 3 × 2 matrices, a hopeless task as a general rule. The work of [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] was motivated by [START_REF] Acerbi | A variational definition of the strain energy for an elastic string[END_REF], the first article to deal with a genuine dimension reduction in nonlinear elasticity via a mathematical convergence result in the case of nonlinearly elastic strings. It was followed in recent years by many, sometimes highly technical developments, including derivations of inextensional bending models and Von Kármán type models, see [START_REF] Friesecke | A theorem on geometric rigidity and the derivation of nonlinear plate theory from three-dimensional elasticity[END_REF][START_REF] Friesecke | The Föppl-von Kármán plate theory as a low energy Γ-limit of nonlinear elasticity[END_REF][START_REF] Friesecke | A hierarchy of plate models derived from nonlinear elasticity by gamma-convergence[END_REF], among others. We are not however aware of the energies found in [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF][START_REF] Shafrir | The membrane shell model in nonlinear elasticity: A variational asymptotic derivation[END_REF] ever being used in a numerical context, even though there are many numerical works on membranes and thin films, see [START_REF] Atai | Coupled deformations of elastic curves and surfaces[END_REF][START_REF] Haseganu | Analysis of partly wrinkled membranes by the method of dynamic relaxation[END_REF] for numerical studies of Pipkin's model using a differential equation approach, including the case of a pressurized membrane [START_REF] Steigmann | Eliza Haseganu's analysis of wrinkling in pressurized membranes[END_REF], see also more or less ad hoc models designed for simplicity or efficiency, for instance in [START_REF] Bonet | Finite element analysis of air supported membrane structures[END_REF][START_REF] Le Maître | Large displacement analysis for ideally flexible sails[END_REF][START_REF] Le Maître | Unsteady model of sail and flow interaction[END_REF]. The modeling and numerical simulation of nonlinear membranes is also attracting increasing interest for materials with more sophisticated material response than just nonlinear elasticity, see [START_REF] Weinberg | A geometrically exact thin membrane model -investigation of large deformations and wrinkling[END_REF] for finite strain, viscoelastic membranes.

One purpose of the present article is to advocate the use of a rigorously derived membrane energy to perform computations that are grounded on an indisputable three-dimensional model and are still efficient. Naturally, we then have to work with what is given by the asymptotic dimensional reduction procedure, and not with an ad hoc energy. This involves a little bit of mathematics, which needs to be carefully done, see section 3 below.

This article is organized as follows: We first briefly present the results of [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] and [START_REF] Shafrir | The membrane shell model in nonlinear elasticity: A variational asymptotic derivation[END_REF]. We consider a three-dimensional hyperelastic homogeneous cylinder of thickness 2ε > 0 made of a given Saint Venant-Kirchhoff material. The body is subjected to a dead loading body force density and a constant pressure differential on its upper and lower surfaces, and a boundary condition of place on its lateral surface. The three-dimensional nonlinear elasticity equilibrium problem is formulated as a minimization problem for the total energy of the body.

Using Γ-convergence arguments, H. Le Dret and A. Raoult showed that deformations that almost minimize the three-dimensional total energy converge when the thickness ε of the body goes to zero towards deformations that minimize a nonlinear membrane energy, see [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF]. The convergence takes place in a rescaled weak W 1,p sense. The limit problem is two-dimensional, with values in R 3 .

The limit two-dimensional nonlinear membrane energy is computed in two steps: First minimize the bulk stored energy function with respect to the third column vector of the deformation gradient-this step produces a function W 0 on the space of 3 × 2 matrices-then take the quasiconvex envelope of W 0 . In the special case of a Saint Venant-Kirchhoff material, an explicit formula for this

A c c e p t e d m a n u s c r i p t

quasiconvex envelope QW 0 is available. This formula is expressed in terms of the right singular values of the membrane deformation. In [START_REF] Shafrir | The membrane shell model in nonlinear elasticity: A variational asymptotic derivation[END_REF], in addition to the case of curved membranes, the zero-thickness limit of a constant live pressure loading term is also computed.

In section 3, we present a conforming finite element approximation of the membrane problem. We consider P 1 and Q 1 discretizations of the three Cartesian components of the deformation. In the P 1 case, deformations are approximated by piecewise affine, globally continuous functions on a triangulation of the domain. In the Q1 case, which is appropriate for rectangular membranes, a rectangular mesh is used with globally continuous approximations that are piecewise of partial degree less or equal to one on the rectangular elements. We prove the weak-W 1,4 convergence of the approximate solutions toward a solution of the continuous minimization problem.

The choice of available numerical methods to solve our FE problem is rather limited since the problem under study is highly nonlinear and the membrane stored energy function is only of class C 1 . Consequently, a method relying on second derivatives of the total energy such as the Newton-Raphson method cannot be appropriate. On the contrary, the nonlinear conjugate gradient method with the Polak and Ribière variant seems to be well adapted to our problem. The convergence of the algorithm is guaranteed for a convex functional, which is only the case here when the pressure differential is zero. There is however a slight difficulty in computing the gradient of the stored energy function. We adapted Ball's results concerning the differentiability of frame-indifferent, isotropic functions on the space of n × n square matrices, see [START_REF] Ball | Differentiability properties of symmetric and isotropic functions[END_REF].

In section 4, we present various numerical tests. Both P 1 and Q 1 elements are alternatively used. The first test is a circular membrane subjected to an upward pressure differential and clamped on its boundary. It should be noted that in our formulation, there as absolutely no need to track the deformed normal vector in order to take into account the live loading pressure differential. This is exemplified by the bubble-like deformation computed in this test which bulges out of the supporting circle. Note again that we have not seen our formulation of the pressure term used in a numerical context.

Next, we perform a few tests taken from [START_REF] Stanuszek | FE analysis of large deformations of membranes with wrinkling[END_REF]: a rectangular airbag inflated by an inner pressure and a square membrane attached by its four corners and subjected to a vertical point force applied at its center. As opposed to [START_REF] Stanuszek | FE analysis of large deformations of membranes with wrinkling[END_REF], our model cannot capture wrinkles in detail, because wrinkles are filtered out in the Γ-limit process, which in turn leads to a well-posed limit minimization problem. Such is the nature of weak convergence. However, wrinkled regions are captured. They are the membrane areas where the deformation gradient lies a region of 3 × 2 space where relaxation occurs, i.e., the quasiconvex envelope is such that QW 0 < W 0 . This occurs in compression when at least one of the singular values is

A c c e p t e d m a n u s c r i p t

less than 1 and the other one is not too large, see [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] for details.

The last two tests are in the context of the modeling of fabrics. The first test is a square piece of fabric attached at its center and subjected to a vertical dead loading body force and the second is a tablecloth with no displacement allowed on the table surface. The obtained deformations develop folds and conical points where the normal vector is ill-defined, but this is of no consequence for our handling of the live pressure term.

Part of the results of this article concerning P 1 elements and without pressure differential were announced in [START_REF] Kerdid | Approximation numérique d'un problème de membrane non linéaire[END_REF].

THE CONTINUOUS PROBLEM

Let us briefly outline the results of [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] and [START_REF] Shafrir | The membrane shell model in nonlinear elasticity: A variational asymptotic derivation[END_REF], to which we refer the reader for more details. Let ω be an open, bounded subset of R 2 with Lipschitz boundary. For all ε > 0, we consider a hyperelastic homogeneous body occupying the reference configuration Ω ε = ω × ]-ε, ε[. We assume that the stored energy function of this body is a function W : M 3 → R which is continuous, coercive and satisfies growth conditions for an exponent p ∈ ]3, +∞[, where M 3 is the space of real 3 × 3 matrices. We furthermore assume that the body is subjected to a dead loading body force density f and to a constant pressure differential ε∆p on its upper and lower surfaces, which is a live load, that is to say a spatially constant pressure p + ε on the upper surface and another spatially constant pressure p - ε on the lower surface such that p + εp - ε = ε∆p. The equilibrium problem for this body may be formulated as a minimization problem for the energy

J ε (φ ) = Ω ε W (∇φ ) dx - Ω ε f ε • φ dx -P ε (φ ), (1)
where

P ε (φ ) = Ω ε π ε det ∇φ + 1 3 ∇π ε • (cof ∇φ T φ ) dx, (2)
over a set of admissible deformations φ belonging to an appropriate Sobolev space and satisfying given boundary conditions of place on part of the lateral boundary.

Here, π ε is a C 1 -function on Ωε that takes the values p ± ε on the upper and lower surfaces. The term P ε appearing in the energy accounts for the pressure load, see [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]. Note that this term incorporates the fact that a pressure load is a live load that follows the normal vector to the deformed body, without having to keep track of this normal vector. Dead loading tractions on the upper, lower and lateral surfaces can also easily be added as well as boundary conditions of place on part of the boundary.

A c c e p t e d m a n u s c r i p t

In [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] and [START_REF] Shafrir | The membrane shell model in nonlinear elasticity: A variational asymptotic derivation[END_REF], see the latter for the pressure term, Le Dret and Raoult proved that a rescaled version of the above three-dimensional energy Γ-converges when the thickness 2ε of the membrane goes to zero in the sense of the weak topology of W 1,p (Ω; R 3 ), thereby showing that minimizing deformations converge, in an appropriate sense, toward solutions of a two-dimensional minimization problem. The limit, two-dimensional nonlinear membrane problem is described as follows.

Let M 3,2 be the space of real 3 × 2 matrices. If z α , α = 1, 2, are two vectors in R 3 , we note (z 1 |z 2 ) the matrix of M 3,2 whose columns are the vectors z α . For all F = (z 1 |z 2 ) ∈ M 3,2 and z ∈ R 3 , we note (F|z) the matrix whose first two columns are z α , and third column is z and write (z|F) with a similar convention. We now define a function W 0 : M 3,2 → R by

W 0 (F) = inf z∈R 3 W ((F|z)). (3)
The function W 0 is continuous and coercive. Let QW 0 be its quasiconvex envelope, see [START_REF] Dacorogna | Direct Methods in the Calculus of Variations[END_REF]. We introduce the space of admissible membrane displacements

Φ M = ψ ∈ W 1,p (ω; R 3 ); ψ(x 1 , x 2 ) = (x 1 , x 2 , 0) T on ∂ ω , (4)
(this is for the case of a boundary condition of place on the whole lateral surface ∂ ω × ]-ε, ε[ in the 3D problem, other conditions are enforced accordingly). The limit nonlinear membrane energy is then defined by

J(ψ) = 2 ω QW 0 (∇ψ) dx 1 dx 2 - ω f •ψ dx 1 dx 2 - ∆p 3 ω det(∂ 1 ψ|∂ 2 ψ|ψ) dx 1 dx 2 ,
(5) for all ψ ∈ Φ M , where f is a 2D-body force resultant density obtained from f ε .

The limit minimization problem now reads: Find ϕ ∈ Φ M such that

J(ϕ) = inf ψ∈Φ M J(ψ). (6)
The energy functional J is obtained via a Γ-limit process and existence of a solution to problem ( 6) is guaranteed under reasonable technical assumptions, see [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] and [START_REF] Shafrir | The membrane shell model in nonlinear elasticity: A variational asymptotic derivation[END_REF].

In general, it is not possible to compute explicitly the quasiconvex envelope of a given function W 0 . Hence, generically with respect to the bulk material, it is not possible to perform numerical computations using this model. However, in the special case of a Saint Venant-Kirchhoff material, the quasiconvex envelope QW 0 was computed in [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF]. Let us recall that the stored energy function of a Saint Venant-Kirchhoff material assumes the form

W (F) = µ 4 tr((F T F -I 3 ) 2 ) + λ 8 (tr(F T F -I 3 )) 2 ,

A c c e p t e d m a n u s c r i p t

where µ > 0 and λ ≥ 0 are the Lamé constants of the material (here F is a 3 × 3 matrix and I n is the n × n identity matrix). The Saint Venant-Kirchhoff energy is the simplest possible frame-indifferent and isotropic nonlinearly elastic stored energy function with a natural state at the identity. It is well suited to the description of the large deflections that are observed in thin membranes. The intermediate function W 0 is given by

W 0 (F) = µ 4 tr((F T F -I 2 ) 2 )+ λ µ 4(λ + 2µ) h(F) 2 + 1 8(λ + 2µ) [λ h(F)-(λ +2µ)] + 2 ,
with h(F) = tr(F T F -I 2 ) (here F is a 3 × 2 matrix, thus F T F is 2 × 2) and t + denotes the positive part of t. Finally, the quasiconvex envelope QW 0 is expressed in terms of the right singular values 0 ≤ v 1 (F) ≤ v 2 (F) of F (i.e. the eigenvalues of

√ F T F) as QW 0 (F) = Φ(v 1 (F), v 2 (F)) = E 8 [v 2 (F) 2 -1] + 2 + E 8(1-ν 2 ) [v 1 (F) 2 + νv 2 (F) 2 -(1 + ν)] + 2 + E 8(1-ν 2 )(1-2ν) [ν(v 1 (F) 2 + v 2 (F) 2 ) -(1 + ν)] + 2 , (7)
where

E = µ(3λ +2µ) λ +µ
is the Young modulus and ν = λ 2(µ+λ ) the Poisson ratio of the Saint Venant-Kirchhoff material under consideration, see [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] and [START_REF] Dret | Quasiconvex envelopes of stored energy densities that are convex with respect to the strain tensor[END_REF]. It turns out that the function QW 0 is convex on M 3,2 , therefore the quasiconvex envelope of W 0 is also its convex envelope in this case.

In the case of the Saint Venant-Kirchhoff material, we thus have a completely explicit expression of the membrane energy that is rigorously derived from threedimensional nonlinear elasticity by means of a convergence result. The pressure term P(ψ) = ∆p 3 ω det(∂ 1 ψ|∂ 2 ψ|ψ) dx 1 dx 2 also incorporates the fact that we are dealing with a live load that follows the normal to the deformed surface, see [START_REF] Shafrir | The membrane shell model in nonlinear elasticity: A variational asymptotic derivation[END_REF]. This feature is be very advantageous for numerical simulations since the deformed normal vector does not explicitly enter this formula. There thus is no need to recalculate this vector at each step of the approximation procedure, as opposed to [START_REF] Stanuszek | FE analysis of large deformations of membranes with wrinkling[END_REF] for example.

THE DISCRETE PROBLEM

Finite element discretization

We use P 1 and Q 1 discretizations. Let us describe the P 1 case since the Q 1 case is entirely similar at this stage. Let τ h be a regular affine family of triangulations

A c c e p t e d m a n u s c r i p t

covering the domain ω. We discretize the three Cartesian components of the deformation using continuous P 1 finite elements. The discrete space of admissible displacements is given by

Φ h M = ψ h ∈ C 0 ( ω; R 3 ), ψ h (x 1 , x 2 ) = (x 1 , x 2 , 0) T on ∂ ω, ψ h | K ∈ (P 1 ) 3 ; ∀K ∈ τ h . (8

) We clearly have Φ h

M ⊂ Φ M and the approximation is conforming. We approximate the continuous problem as follows: Find

φ h ∈ Φ h M such that J(φ h ) = inf ψ h ∈Φ h M J(ψ h ). (9)
This is a convergent approximation scheme in the following sense.

Proposition 3.1 The weak limit points of the sequence φ h in W 1,4 are minimizers of problem [START_REF] Bonet | Finite element analysis of air supported membrane structures[END_REF].

Proof. We first show that the functional J is sequentially weakly lower semicontinuous on W 1,4 (ω; R 3 ). The first two terms

J(ψ) -P(ψ) = 2 ω QW 0 (∇ψ) dx 1 dx 2 - ω f • ψ dx 1 dx 2
combine to form a convex, strongly continuous functional, hence weakly lower semicontinuous functional on the space W 1,4 (ω; R 3 ). For the pressure term, we observe that the components of the vector ∂ 1 ψ ∧ ∂ 2 ψ are null Lagrangians. Since we are in dimension 2, it is well-known that the mapping ψ → ∂ 1 ψ ∧ ∂ 2 ψ is sequentially weakly continuous from W 1,4 (ω; R 3 ) into L 2 (ω; R 3 )), see e.g. [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]. Let us consider a weakly convergent sequence ψ k ψ in W 1,4 (ω; R 3 ). By Rellich's theorem, we have that ψ k → ψ in L 2 (ω; R 3 ) and since det((u|v|w)) = (u ∧ v) • w, it clearly follows that

P(ψ k ) = ∆p 3 ω (∂ 1 ψ k ∧ ∂ 2 ψ k ) • ψ k dx 1 dx 2 → P(ψ), when k → +∞.
Since the functional J is coercive over Φ M , it follows that the sequence φ h is bounded. Let h be a subsequence such that φ h converges weakly in W 1,4 (ω; R 3 ) to a limit point φ .

Let us now take ψ ∈ Φ M . Since the triangulation family is regular, there exists a sequence ψ h ∈ Φ h M such that ψ h → ψ strongly in W 1,4 (ω; R 3 ). Therefore

J(φ ) ≤ lim inf h →0 J(φ h ) ≤ lim h →0 J(ψ h ) = J(ψ),
hence the result.

A c c e p t e d m a n u s c r i p t

Remark 1 By the Rellich-Sobolev embeddings, we also have strong convergence in C 0,α ( ω; R 3 ) for all α < 1/2. However, even if ∆p = 0 and the functional J is convex, it is not strictly convex. Therefore, we cannot deduce strong W 1,4 convergence in general, nor any convergence rate. Note also that the above weak convergence holds true in W 1,p , not only for the Saint Venant-Kirchhoff material, but for a general material satisfying the hypotheses of [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] with p > 3, since quasiconvexity is sufficient to ensure the weak lower semicontinuity of the elastic term in the energy. Besides, the weak lower semicontinuity result can also essentially be viewed as a consequence of the Γconvergence result.

The Polak-Ribière nonlinear conjugate gradient algorithm

The boundary value problem underlying problem ( 6) is highly nonlinear and so is its finite element counterpart. In order to perform numerical computations, we thus need a nonlinear method. The Newton-Raphson method requires the Hessian of J on Φ h M , but the function QW 0 is not of class C 2 . Therefore, this method, or any method relying on second derivatives of J, is not appropriate. Furthermore, the problem is naturally set as an energy minimization problem. We thus use the nonlinear conjugate gradient method with the Polak and Ribière variant, see [START_REF] Saïdi | Analyse mathématique et numérique de modèles de structures intelligentes et de leur contrôle[END_REF], which gives the best results. Let us recall that the conjugate gradient algorithm at the kth iteration is given by:

ψ k+1 h = ψ k h -ρ k d k , d k = ∇J(ψ k h ) -β k-1 d k-1 , β k = - [∇J(ψ k+1 h )-∇J(ψ k h )]•∇J(ψ k+1 h ) ∇J(ψ k h ) , ρ k = d k •∇J(ψ k h ) d k 2 .
(10)

Note that the limit membrane energy J is convex for ∆p = 0 in which case the conjugate gradient iterations converge to a global minimum. It is not necessarily convex for ∆p = 0, and the nonlinear conjugate gradient may conceivably get trapped at a local minimum.

Even though the method itself is easy to use, the computation of ∇J(ψ h ) entails some difficulties. Indeed, since QW 0 is a frame-indifferent, isotropic function defined on M 3,2 via the right singular values of its argument, viz. formula [START_REF] Friesecke | A theorem on geometric rigidity and the derivation of nonlinear plate theory from three-dimensional elasticity[END_REF], the computation of ∂ (QW 0 ) ∂ F , which is a 3 × 2 matrix, requires an adaptation of known results on the differentiation of frame-indifferent, isotropic functions on the space of square matrices M n , see [START_REF] Ball | Differentiability properties of symmetric and isotropic functions[END_REF].

A c c e p t e d m a n u s c r i p t 3.3 Computation of the gradient of J

Let (λ j ) j=1,...,N h be the components of ψ h in the shape function basis (Θ j ) j=1,...,N h of our P 1 or Q 1 finite element space. We clearly have

∂ J ∂ λ j (ψ h ) = ω ∂ (QW 0 ) ∂ F (∇ψ h ) : ∇Θ j dx 1 dx 2 - ω f • Θ j dx 1 dx 2 - ∆p 3 ω λ k λ l det((∂ 1 Θ j |∂ 2 Θ k |Θ l )) + λ i λ l det((∂ 1 Θ i |∂ 2 Θ j |Θ l )) + λ i λ k det((∂ 1 Θ i |∂ 2 Θ k |Θ j )) dx 1 dx 2 , ( 11
)
for j = 1, . . . , N h , where A : B = tr(A T B) denotes the standard dot product on M 3,2 .

We compute the 3 × 2 matrix ∂ (QW 0 ) ∂ F (∇ψ h ) via the following lemma.

Lemma 3.2 The matrix ∂ (QW 0 ) ∂ F (F) has the following explicit form:

∂ (QW 0 ) ∂ F (F) = R   0 0 ∂ 1 Φ(v 1 (F), v 2 (F)) 0 0 ∂ 2 Φ(v 1 (F), v 2 (F))   S, (12)
where (v 1 (F), v 2 (F)) are the right singular values of the matrix F, R ∈ O(3) and S ∈ O(2) are orthogonal matrices such that R T diagonalizes FF T ∈ M 3 and S diagonalizes F T F ∈ M 2 both with eigenvalues in increasing order, and Φ is given by formula [START_REF] Friesecke | A theorem on geometric rigidity and the derivation of nonlinear plate theory from three-dimensional elasticity[END_REF].

Proof. Let us first recall the differentiation result of [START_REF] Ball | Differentiability properties of symmetric and isotropic functions[END_REF]. Let Z : M n → R be a function such that Z(RFS) = Z(F) for all R, S ∈ O(n). It is known that there exists a symmetric function Ψ : R n + → R such that Z(F) = Ψ(v 1 (F), . . . , v n (F)) where (v i (F)) i=1,...,n are the singular values of F. In [START_REF] Ball | Differentiability properties of symmetric and isotropic functions[END_REF], Ball showed that Z is of class C 1 if and only if Ψ is of class C 1 . It is an easy consequence of Ball's results that in this case,

∂ Z ∂ F (F) = R diag ∂ 1 Ψ(v 1 (F), . . . , v n (F)), . . . , ∂ n Ψ(v 1 (F), . . . , v n (F)) S,
where R, S ∈ O(n) are orthogonal matrices occurring in the singular value decomposition of

F = R diag(v 1 (F), . . . , v n (F))S with 0 ≤ v 1 (F) ≤ • • • ≤ v n (F).
Indeed, for any couple of orthogonal matrices R and S, if we let Z R,S (F) = Z(RFS) and differentiate, we obtain that

∂ Z R,S ∂ F (F) = R ∂ Z ∂ F (RFS)S = ∂ Z ∂ F (F),

A c c e p t e d m a n u s c r i p t

hence the result by choosing the right orthogonal matrices. We now recall a few simple facts about the singular value decomposition. By the polar factorization lemma, any F ∈ M n can be written as

F = QU = V Q, where U = √ F T F and V = √ FF T are symmetric, positive and Q ∈ O(n).
The symmetric matrices U and V are unique, the orthogonal matrix Q is unique if det F = 0, it is nonunique otherwise. Let S ∈ O(n) be such that U = S T diag(v i (F))S. Then we have F = (QS T ) diag(v i (F))S, hence a characterization of the orthogonal matrix S to the right of the singular value decomposition of F as a matrix that diagonalizes F T F. Similarly, the orthogonal matrix R to the left of the singular value decomposition of F is a matrix whose transpose diagonalizes FF T , and we have the relation R = QS T . Furthermore, we can always choose these matrices in such a way that the singular values are arranged in increasing order.

Let us go back to the case at hand. Let F be a 3 × 2 matrix. It was established in [START_REF] Dret | The quasiconvex envelope of the Saint Venant-Kirchhoff stored energy function[END_REF] that QW 0 (F) = QW ((F|0)), where QW is the quasiconvex envelope of the Saint Venant-Kirchhoff energy density itself. Since this energy is left and right invariant by O(3), we also have

QW 0 (F) = QW ((0|F)), (13)
Therefore, it follows that the column vectors satisfy

∂ (QW 0 ) ∂ F (F) α = ∂ (QW ) ∂ G ((0|F)) (1+α) , (14)
for α = 1, 2. By [START_REF] Dret | The quasiconvex envelope of the Saint Venant-Kirchhoff stored energy function[END_REF], we have an explicit representation of the quasiconvex envelope of the Saint Venant-Kirchhoff density of the form QW (z|F) = Ψ(v(z|F)) with Ψ : R 3 + → R. Now, by the general result of [START_REF] Ball | Differentiability properties of symmetric and isotropic functions[END_REF], it follows that

∂ (QW ) ∂ G ((0|F)) = R diag ∂ 1 Ψ(0, v 1 (F), v 2 (F)), . . . , ∂ 3 Ψ(0, v 1 (F), v 2 (F)) S, (15)
where R T , S ∈ O(3) respectively diagonalize (0|F)(0|F) T and (0|F) T (0|F), with the requisite order for the singular values. Now we clearly have

(0|F)(0|F) T = FF T and (0|F) T (0|F) = 0 0 0 F T F .
Therefore, if S ∈ O(2) diagonalizes F T F with the eigenvalues in increasing order, it follows that S = 1 0 0 S , since 0 is an eigenvalue of the 3 × 3 matrix (0|F) T (0|F).

A c c e p t e d m a n u s c r i p t

Now, we just have to identify the last two column vectors in formula [START_REF] Dret | Quasiconvex envelopes of stored energy densities that are convex with respect to the strain tensor[END_REF]. This is easy since

∂ (QW ) ∂ G = R   ∂ 1 Ψ 0 0 0 ∂ 2 Ψ 0 0 0 ∂ 3 Ψ   1 0 0 S = R    ∂ 1 Ψ 0 0 ∂ 2 Ψ 0 0 ∂ 3 Ψ S    = R   ∂ 1 Ψ 0 0   0 0 ∂ 2 Ψ 0 0 ∂ 3 Ψ   S   .
The conclusion follows from the fact that [START_REF] Dret | The quasiconvex envelope of the Saint Venant-Kirchhoff stored energy function[END_REF].

Ψ(0, v 1 , v 2 ) = Φ(v 1 , v 2 ) for all 0 ≤ v 1 ≤ v 2 , see
Remark 2 The expressions of ∂ 1 Φ(v 1 , v 2 ) and ∂ 2 Φ(v 1 , v 2 ) are given by ∂ 1 Φ(v 1 , v 2 ) = E 2(1 -ν 2 ) v 1 [v 2 1 + νv 2 2 -(1 + ν)] + + Eν 2(1 -ν 2 )(1 -2ν) v 1 [ν(v 2 1 + v 2 2 ) -(1 + ν)] + and
∂ 2 Φ(v 1 , v 2 ) = E 2 v 2 [v 2 2 -1] + + Eν 2(1 -ν 2 ) v 2 [v 2 1 + νv 2 2 -(1 + ν)] + + Eν 2(1 -ν 2 )(1 -2ν) v 2 [ν(v 2 1 + v 2 2 ) -(1 + ν)] + for all 0 ≤ v 1 ≤ v 2 ,
in view of formula [START_REF] Friesecke | A theorem on geometric rigidity and the derivation of nonlinear plate theory from three-dimensional elasticity[END_REF]. Consequently, to compute the part of the gradient pertaining to QW 0 , we only need to compute the two singular values, which entails solving a second degree equation, and a couple of corresponding orthogonal matrices.

A c c e p t e d m a n u s c r i p t

However, both dead loading and pressure terms require numerical integration in the P 1 as well as Q 1 cases. In particular, in the Q 1 case, the singular values and orthogonal matrices must be evaluated at Gauss points in each element. In order to understand the limit model ( 5) in terms of actual membrane parameters, it should be emphasized that ( 5) is a rescaled, zero-thickness limit model. In particular, it does not involve the actual nonzero thickness of a given membrane of interest. For instance, the pressure differential ∆p stands for the limit of ε -1 ∆p ε when ε goes to zero, where ∆p ε is the actual pressure differential and ε the half-thickness of the membrane. Likewise, body force resultants are rescaled limits.

Now there are infinitely many different ways of embedding an actual nonzero thickness membrane into a family of membranes with vanishing thickness. We can however agree that such a sequence should be made of the same bulk material, hence have material coefficients that are independent of ε, as we have done in the asymptotic analysis. In this context, the computed results for a limit pressure differential ∆p of 1.0e+9 as in the first test below, correspond to an actual pressure differential of 50 kPa for a 0.1 mm thick membrane, and so forth. We will only give the rescaled values below, unless otherwise specified.

Let us first present P 1 tests. The first test (Figure 1) is a circular membrane subjected to an upward pressure differential and clamped on its boundary. A bubble-like deformation is thus created and it is apparent that the live loading character of the pressure is well taken into account, without having to explicitly track the deformed normal vectors. Our second P 1 test is an airbag: a pillow-like structure subjected to an outward pressure differential that inflates it. To compute it, we use two equal square reference domains corresponding to the upper and lower parts of the airbag, with the condition that horizontal displacements agree and vertical displacements vanish on their common boundary. It is instructive to compare this test with the corresponding one in [START_REF] Stanuszek | FE analysis of large deformations of membranes with wrinkling[END_REF]. Because of the relaxed nature of the membrane energy we use, there is no need for a special treatment of the areas located near the middle of the airbag sides, where wrinkling occurs. Such wrinkling is smoothed out by the Γ-limit process, as a result of weak convergence. On the other hand, we do not capture wrinkle details. However, our model can predict wrinkled areas which correspond to areas where the deformation gradient lies in the subset of M 3,2 in which the energy is relaxed. See [START_REF] Dret | The nonlinear membrane model as variational limit of three-dimensional nonlinear elasticity[END_REF] for a description of this subset in terms of the singular values of the deformation gradient.

!!"# !! !$"# $ $"# ! !"# !!"# !! !$"# $ $"# ! !"# $ $"# ! !"#
! !"# !"$ !"% !"& ' '"# '"$ '"% '"& # ! !"( ' '"( # !!"$ !!") !!"# !!"' ! !"' !"# !") !"$
Let us now show the results of a few Q 1 tests. First is the same airbag as before.

A c c e p t e d m a n u s c r i p t

We next compute again the same airbag, with an added point force f 1 = f 2 = 0 and f 3 = -1.0e+9, slightly off-center on the top surface. Back to P 1 -tests, a square membrane attached by its four corners and subjected to a vertical point force applied at the center. Our next tests are in the context of the modeling of fabrics. We take values for the Young modulus E = 2,500 Pa and Poisson ratio ν = 0.01 that are characteristic of cotton fabric, see [START_REF] Kang | Drape simulation of woven fabric using finiteelement method[END_REF]. These values are rescaled with respect to the thickness since the bulk Young modulus of cotton is of the order of 8 GPa.

A c c e p t e d m a n u s c r i p t

In Figure 6, we consider a square piece of fabric attached at its center and subjected to a vertical dead loading body force f 3 = -1,000, slightly counterbalanced by an upward pressure differential ∆p = 100. The lateral sides are hanging free. Our last test is a tablecloth. The material constants and applied forces are the same as above, and the displacement is set to zero on the center square that represents the table, i.e., the fabric is not allowed to slide across the table like an actual tablecloth.

A c c e p t e d m a n u s c r i p tFigure 1 .Figure 2 .

 12 Figure 1. Circular membrane with Young modulus E = 2.7e+2 MPa and Poisson ratio ν = 0.4, corresponding to an elastomer, (rescaled) pressure differential ∆p = 1.0 GPa, zero dead loading body force, 3,243 degrees of freedom.The convergence of the algorithm is as follows

Figure 3 .

 3 Figure 3. P 1 airbag, same mechanical data except pressure differential ∆p = 1.0e+7, 5,046 degrees of freedom.

Figure 4 .

 4 Figure 4. Gradient norm/initial gradient norm ratio plotted against iteration number.

Figure 5 .

 5 Figure 5. Q 1 airbag.

Figure 6 .

 6 Figure 6. Q 1 airbag with point force.

Figure 7 .

 7 Figure 7. Square membrane, with E = 2.1e+8, ν = 0.25 and central point force f 1 = f 2 = 0 and f 3 = -1.0e+10.

Figure 8 .

 8 Figure 8. Q 1 square piece of fabric.

  Figure 9. Q 1 tablecloth.

Note that in the P 1 case, deformation gradients are constant element-wise. Hence the contributions of QW 0 to the energy and its gradient are computed exactly.

A c c e p t e d m a n u s c r i p t 5 Acknowledgement

The authors wish to thank the anonymous referee for pointing out relevant literature to them and improving the manuscript.