Carmen Calderer

The dynamical behaviour of nonlinear elastic spherical shells

Introduction

We consider the problem of the radial motion of a thick spherical shell under a constant pressure difference between the inner and the outer surfaces, and under zero body forces. We assume that the material of the shell is homogeneous, isotropic, hyperelastic and incompressible.

The problem of inflation under pressure of an arbitrary (i.e. nonspherical) compressible hyperelastic body has been studied by Ball [14]. He applied the concavity method of Knops, Levine and Payne [ 15] to show that weak solutions to the equations of motion with appropriate initial data do not exist for all positive time. The method of proof suggests, but does not establish, that the mechanism for non-existence is blow-up of the displacement in finite time. However, it seems that this method does not apply if the material is almost incompressible. This is related to the fact that the concavity hypothesis made on the strainenergy function W(F) in [14] cannot hold if the material is almost incompressible. Also the set of initial data corresponding to non-existence is not easy to characterize. In this paper we show that for the radial motion of an incompressible spherical or cylindrical shell the special geometry enables us to determine completely the behaviour of all solutions, and in particular to describe exactly for which initial data the solution exists for all time, and for which the solution becomes unbounded in finite time. This description is possible because the equation of motion reduces to a conservative autonomous ordinary differential equation in the plane.

The constitutive properties of the material enter the differential equation through a function g(x 3 , c5) of the tangential strain x and a parameter c5 giving the thickness of the shell in reference configuration. The strain-energy function W(x 3 ) enters into the expression for g(x 3 , 8) in a complicated way.

The phase-plane analysis is carried out initially under assumptions on the function g(x 3 , 8), rather than W(x 3 ). We make hypotheses on the growth of this function for large and small values of x, on its sign, on the number of its critical points, and on the number of roots

x > 1 of the equation g(x 3 , 5) = P for a given P > 0.

The growth conditions on g are related to the growth properties of W( F) for large F. The importance of such growth conditions in nonlinear elasticity was first emphasized by Antman [16] . In [14] Ball interprets the concavity inequality in terms of the slow growth of W( F) as some norm of the deformation gradient F tends to infmity, this condition being appropriate for a "weak material". This slow growth of W(F) is responsible for the global non-existence of solutions. Following this idea, this paper shows that for "strong materials", those for which g(x 3 , 5) grows faster than certain powers of x for small and large x, all solutions exist and remain bounded for all time, while "weak" materials have unbounded solutions that blow up in finite time, for some boundary pressures and some initial data. (As an example of these growth conditions, let us consider the isotropic strain energy function a(A 1 , X 2 , X 3 ) = xr + X~ + X~, where X 1 , X:z. X 3 describe the principal stretches, and a> 0.

For the radial deformation under study, X 1 = X, X 2 = X -112 = X 3 and x = X -vz, where X is the radial stretch. Then the material is weak for x large, if a< 3, and strong if a> 3, and it is strong for x small for all a> 0).

Since the dynamical description is given in terms of g(x 3 , 5), in Section 3 we consider the problem of determining for 5 > 0 a strain energy function corresponding to a given g(x 3 , 5).

Using finite difference equation methods, we find a unique solution given by an infinite series.

The hypothesis on the sign of g(x 3 , 5) are interpreted in terms of constitutive restrictions on the strain-energy function.

Critical situations arise for some particular values of the external pressure P crit at which the number of equilibrium solutions decreases by two asP passes through P crit• Under some hypotheses on g(x 3 , 5), we, in Section 4, use invariant manifolds to determine when P = P crit and the behaviour of the solutions in the neighbourhood of the relevant equilibrium point.

In Section 5 we consider the dynamical equations of a membrane, for which it is possible to carry out the same analysis. We show that the solutions of the equation of motion are the uniform limit on finite time intervals of the solutions for the thick shell as 5 tends to zero. In this case we construct some examples to show that some of the phase-plane results obtained in Section 2 can occur under constitutive hypotheses on the strain-energy function such as convexity of W(F). (For incompressible materials this is not an unreasonable assumption).

In Section 6 we summarise the corresponding results for a cylindrical shell.

Section 1

We consider a Cartesian system and a set of spherical coordinates (R, e, ~)in IR 3 . Let the reference configuration of the shell be the domain n c IR 3 , n = {(R, e, ~), R1 <R <Rz},

(1.1)
where R 1 , R 2 denote the radii of the inner and the outer surfaces of the shell, respectively.

Let the position (r(t), 8(t), <fJ(t)) at timet of the particle labelled by (R, e, cl>) in the reference configuration be given by

r 3 (t) = R 3 + ,f(t)-RL O(t) = e, <P(t) = cl>.
(1.2) (We refer to [2] for the following elements in continuum mechanics). The physical components of the deformation gradient matrix are

[ R 2 r '] F = diag 7, R.' R , J = det F = I.
In spherical components the right and left Cauchy-Green tensor C, B are given by The principal stretches A 1 , A 2 , A 3 of the deformation (1.2) are

The principal invariants of (1.4) are 2

I = A2 +-• t A' (1.3) (1.4) (1.5) (1.6)
We assume that the solid is hyperelastic and isotropic, i.e. the strain energy function depends on F through the invariants of B, or equivalently, it can be expressed as a symmetric function of the principal stretches Aj, Rl)

T = diag -p+l/lt4+1/l-t4,-p+l/Jt2+1/l-tl"'-p+l/JtR2 +1/1-tl" • r R R r r
(1.10)

The spherical components of the acceleration a for the motion (1.2) are

(1.11)

The balance laws of linear and angular momentum, with zero body force, are given in pointwise form by div T(x, t) = pa(x, t),

T = TT.
The stress vector t measured per unit area of the deformed configuration is

t = Tfi,
(1.12)

(1. [START_REF] Hildebrand | Finite difference equations and simulations[END_REF] where fi denotes the unit outward normal to the boundary of n. We follow the work of Knowles and J acu b [START_REF] Knowles | Finite dynamic deformations of an incompressible medium containing a spherical cavity[END_REF] to deduce the differential equation of the motion of the spherical shell. We write the field equations (1.12) in spherical coordinates From the two last equations we can assert that pis independent of 8 and 1/>, i.e., p = p(r, t).

Differentiating with respect to tin (1.2), we deduce that

.. a § r = -- ar'
where § = ,-1 (2r 1 r~ + r~i\)-!r-1 r1r~.

is the acceleration potential.

(1.14)

(1.15)

(1.16)

(1.17)

We assume that the outer surface is free of traction and the inner one is subjected to a constant pressure p > 0 measured per unit area in the present configuration, T(x,,t)fi(x, t) = -pfi(x, t), = 0, where Hence,

x(r,O,cp) = x(r 1 ,0,cp)

x(r,O,cp) = x(r 2 ,0,cp),

(1.18)

(1.19)

We integrate (1.14) with respect tor using (1.16) and (1.17). Then, we can write (1.19) as

-p( § 1 - §2) + 2 ( ' (~2 2-~:) (\V1 -~22 \V-1) ~ = -p, where §1 =r1r1 +~;L § 1 ( "2 2•• ) 1 -4 4 "2 2 = -2r1r1 + '1'1 -?J2 '1'1• '2
Next we set

x(t) = '~~) > 0, R~ 8 = R~ -1 >O. Then R3 HR, t) = 1 + R! (x
8 )-1/3] [ ( 8 )-1/3 2xx 1 -1 + x 3 + x 2 3 -3 1 + x 3 2 P = -p R 2 ' p 1 x>O, (1.20) (1.21) (1.22) (1.23) (1.24) (1.25) (
We consider the density p as a constant.

or It is important to notice that since Pis constant (I .24) is autonomous. We introduce the one-dimensional strain energy function W(~), a(A.) = W(HA.)).

Differentiating (1.26) and using (1.9), we can write (1.25) as

3 2 •X 3 W'(~) g(x , o) = R2 J6+x' ~-1 d~. p 1 1+6 ..
We seek the solutions of equation (1.24) satisfying the initial conditions

x(O) = x 0 , x(O) = x 0 , for(x 0 ,x 0 )E-~=(O,oo)x R (1.26) (1.27) (1.28)
We consider equation (1.24) for a given g(x 3 , o), o > 0 fixed. We regard g(x 3 , o) as a function of X only, 0 being a parameter. We make the following hypotheses on g(x 3 , o) (1.29):

(A) g E C 2 (0, oo). (B) g(x 3 ,o)>O x>1; g(x 3 ,o)<O x<I.
(C) There exist D ~ C> 0 and 0: > 0, for X large, such that cxa-J .s;;g(x 3 , o) .s;;nxa-3 •

We consider the three cases

(1) o:>3,

(2) 0: < 3, (3) 0: = 3.

(D) There exist M~N>O, and {3>0, for x small, such that-Mx-p .s;;g(x 3 ,8).s;;-Nx-P. .

We consider the two cases

(1) {3~3,

(2) {3 < 3.

Remark. In some examples the exponents o: and {3 may be related. See the examples in Section 5. Moreover, we will make further assumptions later on (a) The number of non degenerate critical points with respect to ~ of g(~, o ).

(b) The number of roots x > 1 of the equation g(x 3 , o) = P, for P > 0 a given constant.

Since we want to analyse equation (1.24) in the phase-plane, it is useful to consider the equivalent first order system of two equations

(1.30) where X = (::)'

(1.31)

(1.32) and (1.28),

x(O) = C:)• (1.33)
It is easy to check that hypothesis (A) on g(x 3 , o) and the smoothness of the functions involved in (1.32) imply that h( •, •) E (C 2 (9)) 2 , and, therefore his Lipschitz in 9 [4], [5].

Proposition:

Let (x 0 ,x 0 ) E9.

The initial boundary value problem (1.30)-(1.33) has a unique maximally defined C 2 solution x(t) on [0, tmax), 0 < tmax o;;;;;oo, satisfyingx 1 (0) =x 0 ,x 2 (0) =x 0 .

Moreover if tmax < 00 , then lx(t)l + lx(t)l ~ oo as t ~ tmax•

Proof:

Local existence and uniqueness of the solution through an initial data point (x 0 ,

x 0 ), follow from standard results in ordinary differential equations.

The last statement will be proved in Section 2, using the energy equation. In particular, we will show that if

lim x(t) = 0, then lim x(t) = -oo. t-.tmax t-+tmax
Therefore, the solution cannot reach a point (0, b) of the boundary of the domain in a finite time.

Energy Equation

Next we want to find the energy equation. We multiply (1.24) by x 2 x and integrate respect

to t obtaining [1 -(1 + : 3 ) - 113 ] x 3 x 2 + G(x 3 , 8)-i Px 3 = E,
where

1 f x3 G(x 3 , 8) = - g(A, 8)dA. 3 1
E is the total energy, which depends on (x 0 , x 0 ). We defme the potential energy

V(x,8) = G(x 3 ,8)-!Px 3 . (1.34) (1.35) (1.36)
Next we want to characterise the equilibrium solutions of the initial boundary value problem ( 1.30)-(1.33). We consider a boundary pressure P > 0. We state the following.

Proposition ( 1.1)

Let P > 0, R 1 > 0, p > 0, 8 > 0, be given constants. We write

I d g(t8) = d~g(t8).
Then:

(i) The equlibrium solutionsx 1 (t) =x,x 2 (t) = 0, tE [0, oo) of(1.30) are the roots of g(x 3 , 8) = P,
(1.37) and they satisfy x> 1. (ii) If g'(x 3 , 8) > 0 then (.X, 0) is a centre and V(x, 8) is a relative minimum. (iii) If g'(x 3 , 8) < 0 then (.X, 0) is a saddle point and V(x, 8) is a relative maximum.

Proof.

(i) This follows from hypothesis (B). (ii), (iii) We consider the eigenvalues of the linearised equation about (.X, 0) of (1.31),

A= -J..L = -3xg'(x3, 8) lvz 2 r1-(1 + :3r 113 J (1.38) Standard results on 0 DE [4],
[5] imply that the equilibrium solution (.X, 0) is a centre of the eigenvalues of the linearised equation are pure imaginary, equal, with opposite sign. This result applies in case (ii).

If g' < 0 then the eigenvalues are real with the same absolute value and opposite sign. By the same results (.X, 0) is a saddle point. This applies in case (iii).

The results on V(x, 8) follow using (1.36) and its derivatives. o

Remark: If g'(x 3 , 8) = 0, then A= 0 = J.1 and the linearisation method fails. In Section 4

we consider this case under further hypotheses on g(.X 3 , 8).

In order to study continuation of the solution of (1.24) for given initial data, we will use the energy equation (1.34). The behaviour of the potential energy V(x, o) for large and small values of x is a basic ingredient in obtaining such results. Therefore, we are going to study the asymptotic growth of V(x, o) for large and small values ofx, for all the cases included in (C) and (D).

Proposition (1.2)

The asymptotic behaviour of V(x, o), for X large, is as follows:

(l) If(C1)holds,then,forxlarge,

.£xQ ~ V(x,o) ~(D + 1 ) xQ.
2a a

(2) If (C3) holds, then, for x large, and a given P > 0,

!(C-P)x 3 ~ V(x,o)<(D-P)x 3 • (3) If (C2) holds, then lim V(x, o) = -1.
x-+oo

P 3 -x 3 Proof:
(l) For a sufficiently large ftxed a> 0 and for sufficiently large x, we can write Using (1.35), (1.36), we have

(D + 1) P V(x o)~ ---xQ --x 3 ' a 3 • Therefore, for a> 3, .E...xQ ~ V(x,o)~(D+ 1 ) xQ.
2a a

(2) and (3) follow similarly.

Proposition (1.3)

The asymptotic behaviour of V(x, o), for x > 0 small, is as follows:

(1.39)

(1.40)

(1.41)

(1) If (Dl) holds, then,

M+ 1 1 N 1 for {3 > 3, {3-3 xfl_ 3 ;;a.V(x,8);;a. 2 ({3-3 )xfl-3 >0, for {3 = 3, lim V(x, 8 ) = 1. x-+0+ -logx (1.42)
(2) If (D2) holds, then

V 0 = lim V(x, 8) x-+0+
(1.43) exists, is finite and positive.

Proof.

We consider a sufficiently small fixed a> 0. For a sufficiently small x > 0, we can apply (D)

to the following inequality

x3 a3 x3 a3 f g(t 8)d~ = f g(~, 8)d~ + f g(t 8)d~ ;;a. f g(~, 8)d~ 1 1 a 3 1 (1.44)
Using (1.35), (1.36), we see that (1.42) follows for {3 =I= 3 and x small. It follows from (1.44) and (B) that co> J:-g(~,8)d~>O.

Then (2) follows immediately. (1), for {3 = 3, follows similarly.

Section 2

Now we analyse the initial-boundary-value problem in the phase-plane under the hypotheses on g(x 3 , 8) previously given. Moreover, we assume

(1) The number of critical points of g(x 3 , 8) is finite.

(

) If g(x 3 , 8) = c, for x E (a, b) 2
, where b >a> I, c > 0, are some given constants, then the boundary pressure P=F c. (The case P = c, can be treated similarly.)

We need the following notation.

Let (x 0 ,x 0 )E,~, and E = E(x 0 ,x 0 ).

We let

Q = {yER:V(y,8) = £}.
It is obvious from Proposition (1.1) under the previous assumptions on g(x 3 , 8), that n is a finite set. We let n = {y 1 , .•. ,Yn}, with, Yi-1 .r;;;.yi, for j.;;;;,n

We state the following proposition.

Proposition 2.1

Let (x 0 ,x 0 )E.9J, E = E(xo,Xo). The intergers r, s are determined by V(x, o), (x 0 , x 0 ) and E(x 0 , x 0 ).

Proof.

1 (a) and 2(a) follow from Propositions (1.2) and (1. Related to the the above proposition we introduce the following notation:

(i) v. = min V(x, o). xE [Ys-1•Ya] (ti) 2 1 Vs = Y-t---, [-1---( 1-+-Y-o.-r-1'-l [ E-V,] .
Defmition.

(i) The positive orbit of the solution x( Then, any initial data (x 0 , x 0 ) E 9 satisfies one of the following possibilities.

(a) If (x 0 , x 0 ) is such that E =I= V(x 28 ), j = I, ... , k, the solution is periodic, with period T given by { [ ( 8)-1/3] }112

x 3 1-I+-

y28 x3 T = 2 dx, JY 28 • 1 E-V(x,8)
for some s ;;;a. I.

(2.3) (b) For (x 0 , x 0 ) such that E = V(x 28 ) for somes E {1, 2, ... , k}, x 28 is a saddle point, and, _ Also lim x(t) = x 28 and lim x(t) = 0.

t-+oo t-+oo

(3) If x 0 does not belong to any of the above intervals then the solution is periodic.

Proof.

(a) We consider Y2s-1 ,Y2 8 such that x 0 E [y 2 .-1 , y 28 ]. Then, the solution (x(t),x(t))

starting at (x 0 , x 0 ) at t = 0 is such that

V(x, o)
1.

.

. X

3.

1. g(x 3 , 6) satisfying hypotheses (Cl), (Dl). 2. V(x, 6) for a given P > 0. 3. Phase-plane diagram.

x(t) E [Y2 8 -1 .Y2sJ for all t E (0, tmax), for otherwise the energy equation would be contradicted. Also 0,..;; x 2 (t) < v1, t E (0, tmax).

Hence, the orbit c+ through (x 0 , x 0 ) is bounded. Moreover, (y 2 . -1 , 0), (y 2 ,, 0) are not equilibrium points by the choice of E. Therefore, the orbit is periodic.

(b) This follows by similar arguments and the fact that (x 2 i> 0) is a saddle point.

Propostion (2.3)

We supppse that (i) (Cl) and (D2) hold, (ii) and (iii) are as in Proposition (2.2).

Then, any initial data (x 0 , x 0 ) E 9 satisfies one of the following possibilities:

(a) If E < V 0 all solutions are bounded. (b) If E ";;> V 0 , then solutions with initial data x 0 > y 1 remain bounded, and solutions with x 0 <y 1 blow up in finite time. In this case lim (x(t), x(t)) = (0,-oo). Proof.

(2.4)

(2.5) (a) By the choice of E this case is included in Proposition {2.2). The same occurs in {b) for

Xo > Yt• {b) If x 0 <y 1 , then x(t) <;y 1 , t E [0, tmax), by defmition of y 1 • This implies that E-V(x(t), c5)'> 0 for t E [0, tmax)•
Therefore, if we suppose that Xo < 0, then x(t) < 0, t E [0, tmax), i.e., x(t) <; Xo is monotonically decreasing. Since

X~~+ x 3 [1-(1 + : 3 r/ 3 ] then lim [E-V(x,c5)] lim x 2 =
x-+0+

x-+ O+ [ ( <' ) )-113]

x 3 1-1 +-x3
By energy equation {1.34), it can not occur that lim (x(t), x(t)) = (0, C), for some constant C.

} 112 = -J: dt. •'(t-(I+,~r•) (E-V(s.6)] {2.8)
{The minus sign is because x(t) < 0).

V(x, o)

1.

2 .

. X X

3.

I. g(x', cS) satisfying hypotheses (Cl), (D2). 2. V(x, c5) for a givenP > 0. 3. Phase-plane diagram.

By the previous statements, it follows that the left integral remains bounded and the right one is unbounded, as t-+ 00 • This contradication shows that tmax < 00 • To obtain estimates for tmax• we consider x 0 sufficiently small and apply (1.44) to the left integral (2.8). We thus get upper and lower bounds for this integral. Retaining the higher order terms in x 0 , using (2.6) and integrating the remaining expression, we see that (2.4) and (2.5) follow immediately.

We have assumed that x 0 < 0. The same conclusions hold in the case x 0 > 0, since this case reduces to the previous one once x becomes negative for x = y 1 • o

Proposition (2.4)

We suppose that (i) (C2) and (D1) hold, (ii) g(x 3 , 6) has 2n + 1 critical points, n;;;;, u, (iii) the roots of g(x 3 , 6) =Pare {.Xi>1, 1 E;;;iE;;;2k+2, kE;;;n}.

Then, any initial data (x 0 ,x 0 ) E 9satisfies one of the following possibilities: (a) For (x 0 ,x 0 ) such that E> V=Max {V(xi), i= 1, ... , 2(k+ 1)}, the solution blows up in finite time.

(b) For (x 0 ,x 0 ) such that E < V(xzk+z) the solution blows up in finite time if x 0 > x2k+Z, and remains bounded if x 0 < x zk + 2 • (c) For (x 0 , x 0 ) such that V>E> V(xzk+z), then if x 0 >x the solution blows up in fmite time and if x 0 < x it remains bounded for all time, where x = max {x E IR+, V(x, 8) = £}. and (1) Ifxo >x 2 1 and.X 0 < 0, or x 0 <x 2 J, then lim x(t) = x 21 and lim x(t) = 0, f...,...oo t-+-oo

If V= V(xzk+z), case (c) is included in (b). (d) ForEE {V(x 21 ),j = 1, ... k + 1}, ifx ~x21 then V(x 21 ) > V(x, 8)},
(2) If x 0 > x 21 ,x 0 > 0 the solution blows up in fmite time. Whenever the solution blows up in finite time, we have the following estimate for !max, provided x 0 is sufficiently large:

!max = If+ ~(xo) {2.14)
where

3C {!__1_ 1 E;;~(xo)E;;3(D+ 1),)!_1_ 1 .
2P ..j P 7 -2o: X o ~ -0:

P P 7 -2o: X o ~ -0:
(2.15)

Proof.

(a) We consider.X 0 >O. Also n= {yt}, and by Proposition (2.1)y 1 E;;x 0 • Then£-V(x, 8) > 0 and x(t) > 0, t E [0, !max).

Hence, x(t) increases monotonically. Also

lim x 3 1 -1 + - = -. [ ( 8 )- 113 ] 8 x~oo x 3 3 (2.16)
We consider (2.8) with positive sign on the right integral. By the previous statements, it follows that the left integral remains bounded and the right one is unbounded, as t _. 00 • This contradiction shows that !max < 00 • To estimate tmax we choose x 0 sufficiently large in the previous integral to apply Proposition (1.2). Using (2.16), taking the higher order terms in x 0 in the remaining expression and integrating, we obtain (2.14) and (2.15).

For any x 0 >x 1 , we can proceed similarly.

If Xo < 0 the arne arguments can be applied after x(t) changes sign in X= Yt. • X

3.

I. q (x', 6), satisfying hypotheses (C2), (Dl). 2. V(x, 6) for a givenP > 0. 3. Phase-plane diagram.

For V 0 < V, (a) If E < V 0 , the conclusions of Proposition (2.4) hold.

(b) If V>E> V 0 ,either (1) x 0 <y 1 , then, the conclusions of Proposition (2.3) hold, or (2)x 0 > y 1 , then, the conclusions of Proposition (2.4) hold. (c) E > V, the solutions are unbounded and blow up in fmite time. If V 0 > V, then the conclusions of Proposition 2.4 hold forE< V 0 and (c) forE> V 0 •

Remark (2.4)

We suppose that (C3) and (Dl) hold. Then (i) If P E;; C, the conclusions of Proposition (2.2) hold, (ii) If P ;;;a. D, the conclusions of Proposition (2.4) hold, where C and Dare given in (1.29).

Section 3

The constitutive properties of the material enter the differential equation through the function g(x 3 , S).

Motivated by this fact we want to find which strain-energy function W(x 3 ) corresponds to a given g(x 3 , S) for S > 0 and W(1) fixed, i.e., we look for a solution of equation (1.27) under our previous hypotheses.

It is interesting to notice, however, that in the derivation of the differential equation

(1.24) and its energy integral (1.34), it is not essential to require the material to be hyperelastic [ 6] , (7] . We construct a solution to (1.27) via a finite difference equation [I3] as follows.

We introduce new variables,

x 3 = I +eY, yE(-oo,oo), w(y) = W'(I + eY). (3.3) (In this section we suppress the dependence of g on 8.) A solution of (3.3) is obtained as follows:

(1) We find a solution WN = wN(y), y E R, for N> I with the following boundary conditions wN(Y) = f.l., yE(-N,-N+.'?). Similarly, for (3.6) using z = y-r 2", we obtain {3.8) and {3.9) give the solution (3.3) satisfying (3.4) for z E R. {2) We take limits in (3.4) as N-+ oo:

We can write from equation (3.3) that

Now we consider (3.5). We let z

= y + r2.'Then r WN(z) = wN(z -r ) + ' Y L ez+(i-r>tg'(l + ez+(J-r>t).

If we take

IJ. = w(-oo); .. w(z) = 1J. + ' Y L ez-KJ:g' {1 + ez=-K.t),, z E (-oo, oo).

K=O

Letting N-+ oo in (3 .9), we again get (3 .1 0).

(3.9) {3.10)

(3.11)
It is easy to check that w given by (3.11) is the unique solution of {3.3). With the original variables, {3.11) becomes {3.12) Now we want to prove uniform covergence of this series in order to integrate it term by term, following classical results in analysis.

Let us consider N ";;!: 1 and ME R large enough, and a flxed x <:. M. Using the mean value theorem we fmd that, for a flxed j

IW'( 3) •W'( 3) ~r(x 3 -1)[ g'(1) + (x 3 -1)K ] X - n X j..,.
5 1)

(1 + 5 )N-1 ( 5 + 2 )( 1 + 5)2N-2 (3.17) ~ r(M 3 -1) r g'(
K(M 3 -1) ] .... 5 (1 + 5)N-1 + (5 + 2)(1 + 5)2N-2 = e(N),
i.e., given e > 0 there exists anN such that the above difference is less thane for n ~ N. o Hence integrating (3.12) for x <.M, we get (3.18)

We have thus proved the following:

Proposition (3.1)
For a given g(x 3 , 5) satisfying hypotheses (A), (B) and for ftxed constants 5 > 0, IJ. E R, the integral equation

2 rx• W'(~) g(x 3 , 5) = R 2 J~3 +o 1: _ 1 • d~, W(1) = 1J. p 1 1 +li ,. (3.19)
has a unique C 2 solution wm, given by (3.18).

Next we want to find what asymptotic growth of W(x 3 ) is implied by hypotheses (C) and (D) ong(x 3 , 5).

Proposition (3.2)

The asymptotic growth of W(x 3 ) for x large, is as follows:

(I) Ifhypotheses(C1)holds, then

(I + 6)<a-3 > 13 ( 3C) I (I + 6)<a- 3 > 13 (

3D)

xa (I + 6 )<a-3)13 _ I D--;; ;> -:y W(x 3 ) ;>xa (I + 6 )<a-3>13 _I C--;-. {3.20)

(2) If {C2) holds, then

I ( 3D) {1 +liP -a.)(N+l)/ 3 -1 J.L -:y W(x3);>xa C--; (I +6)<3-a.)/3 -I +-:yx3 {3.2I)
where

(x3 -I)2 e = e(x3, N) = kl (x3' N) 6{1 + 6)N-l and k 1 (x 3 ,N) = sup ig'(~)l ~E [I, I+ (; 3 +~;N]• (3) If {C3) holds, then J.l. I J.l. -x 3 +(D-C)x 3 ;>-W(x 3 );>(C-D)x 3 +-x 3 r r r (3.22) Proof.
{1), (3) follow from (3.I8) using (CI), (C3), respectively. (2) We approximate the terms with j ;> N, in (3.I8) for aN sufficiently large by e(x 3 , N), and we use (C2). It is important to notice that k 1 tends tog'(I) asNtends to oo. o

Similarly, we can state, for x small, the following.

Proposition (3.3)

We suppose that hypothes (D) hold. Then the asymptotic growth of W(x 3 ), for x small is ~ x-(J .;;;;~ W(x 3 ).;;;;(I +M)x-fJ, where M, N, (j are given in hypotheses (D).

(3.23)

Now we want to discuss the meaning of hypotheses (B) in terms of constitutive restrictions on the strain-energy function. Such restrictions have been widely studied

[6], [7], [8], [9], [10]' [17].
Strong ellipticity or Hadamard's condition [9] for the stability of equilibrium solutions is in some cases physically reasonable. Necessary and sufficient conditions on the isotropic strain-energy function for strong ellipticity to hold are known [ 12] in two dimensions but not in three, and this makes difficult the application of such a condition in specific problems.

It is well known [6], [8], [17], that a convexity of W(F) is not physically reasonable in compressible elasticity because it implies the uniqueness of equilibrium solutions to mixed boundary-value problems. It is, however, a good condition in incompressible elasticity; the strain energy function for a Neo-Neokean material, W( F) = tr( FFT), is convex.

It is easy to check that strong ellipticity applied to our problem implies that a'(A.) > 0 for A.> 1 and a' (A.)< 0 for A.< 1. From this, (B) follows.

We also can state the following.

Proposition (3.4) We suppose that W( F) is convex for FE M 3 x 3 such that det F = 1. If a(A. 1 , A. 2 , A. 3 ) is the corresponding isotropic strain-energy function, then a(A.) = a(A., A. -112 , A. -112 ) is convex on IR+.

Proof.

It follows from [8] that W( F) convex on M 3 x 3 implies that the isotropic strain-energy function of the principal stretches a(A. 1 , A. 2 , A.

~o:a(A., x-1/2' x-1/2) + (1-o:)a(s.t,s.t-1/2 ,s.t-1'2) = o:a(A.) + (1 -o:) a(s.t). So, a(A.) is convex.
o Later, we are going to produce specific examples to show that some of the phase-plane results obtained in Section 2 can occur under convexity of a(A.).

Section 4

The purpose of this section is to study the stability properties of the equilibrium solutions in critical cases. We consider ax> 1 and P > 0 such that g(x 3 ' o) = p and g'(x 3 ' o) = 0.

The eigenvalues of the linearised equation are A= 0 = -IJ., and therefore, the arguments used in Proposition (1.1) do not hold. Following [18], [19] we will compute local invariant manifolds through the equilibrium point.

These values P crtt of P are called "critical pressures" or "bifurcating pressures". Under certain conditions the number of equilibrium solutions decreases by two asP passes through P crit. This problem has been previously studied mainly from a static point of view, cf. [ 11] .

We want to show the existence of an invariant manifold through (x, 0) and also that this manifold can be computed to any degree of accuracy given hypotheses of smoothness on g(x3, c5).

Deimition:

A set S C IR 2 is an invariant manifold for (1.31) if for any (x 0 , x 0 ) E S the solution (x(t), x(t)) of (1.31) starting at (x 0 ,x 0 ) is such that (x(t),x(t))ES, tE [0, T] for some T > 0 sufficiently small.

Proposition ( 4.1)

For a givenP>O, let x> 1 be a solution ofg(x 3 , c5) =Psuch thatg'(x 3 , c5) = 0 andg"(x 3 , c5) * 0, g smooth. Then, there exists an invariant manifold, x 2 = h(xd through (x, 0) defmed for x 1 E [x, x +a), if ..2"(x) ;;a. 0 for x E [x, x +a), for a given a> 0, and for x 1 E (x-b,X], if Sf(x)..;; 0 for x E (x -b,X], for a given b >O.

In this case,

(x) = _[ ____ c5 __ -1,-3-] x3 1 + x3 +3 1 + x3 -3 ' x 1-( 1 + x3) (4.2) 1 3 .
2' (x) =+-(I+ :,)'"(-g(x ,6)].

Proof.

Let x 2 = h(xd be an invariant manifold for equation (1.30). Then Using the variation of constants formulae, we can write the solution of the above equation as

X (X 1 3 Y(xt) = efx'r(x)dx '_ efxr(u)du.f£'(s)ds. • X (4.6)
The conclusion follows by considering the sign of y( x t).

Proposition ( 4.2)

Let the above hypotheses hold. Let e > 0 be a given constant. The invariant manifold through (x, 0) is given locally as follows:

{1) If g"(x 3 , cS) < 0, then h±(x 1 ) is defmed for x 1 such that 0 E;;;x 1 -x <e.

{2) lfg"(x 3 ,cS)>Othenh±(x 1 )isdefinedforx 1 suchthatOE;;;x-x 1 <e.

In both cases

h±(xd = C(x)lx1 -xl 3 ' 2 + O(e 3 ' 2 ),
where lg"(x3, <5)1112 c( x) = v'f x 3 ' 2 -....:..::...~-'----'-'----

[ ( c5 )1/3] 1/2 .
Proof.

I -I +-=3 X (4.7) (4.8)
We consider ( 4.6) for x 1 such that lx 1 -x 0 I <e. We compute all the integrals involved in ( 4.6) developing (s) in a power series about x 1 = x. The result follows by considering both possible signs for g"(?, cS) and using Proposition (4.1). o Proposition ( 4.3)

Let the previous hypotheses hold. Then (x, 0) is unstable.

Proof.

This follows by energy arguments, considering a solution corresponding to E = V(x, l>) starting near (x, 0), taking x 0 > 0 and x 0 < 0 successively, and considering both possible signs of g"(x 3 , c5). o

To get a more specific idea of the asymptotic behaviour of the solutions in the invariant manifold, we consider the following initial value problems:

(a) x 1 = h+(x 1 ), (b) x 1 = h-(x.), x 1 (0) = h+(uo), x 1 (0) = h-(vo)
for given u 0 , v 0 such that lu 0 -xi< e, and lvo -xi< e.

The solution of ( 4.9) for lx 1 -xi < e is as follows.

(I) Ifg"(x 3 ,o)<O,then:

In case (a),

+ -_2_ [ 1 - 1 ] + 0( 1/2) t -C(x) (uo -x)1t2 (x1 -x)1'2 e
and lim X1 =x+.

t-+-oo

In case (b),

---- - + 1/2 2 [ 1 1 ] t -C(x) (vo -x)112 (x. -x)1t2
O(e ) and lim x. =x+.

t-+oo (2) If g"(x 3 , o) > 0, then: In case (a), ---- + + 1/2 2 [ 1 1 ] t -C(x) (x -uo)1/2 (x -x.)1/2
O(e ) and lim x. = x-.

t-++ 00

In case (b),

-t = --- + + 0 €1/2 2 [ 1 1 J C(x) (x-v 0 ) 112 (x-x 1 ) 112 ( )
and lim x. =x.

t-+-00

SectionS

The membrane approximation.

(4.9)

The differential equation describing the radial motion of a membrane is given by [3), [6],

[7]' [START_REF] Jaunzemis | Continuum Mechanics[END_REF] The same equation is obtained using a membrane approximation up to the first order in 6 [ 11] . In this case x is the average of the tangential stretch divided by the initial thickness of the shell. The energy equation is

x 2 + V(x) = E, where 2 V(x) = -W(x 3 )-px 3 pRj is the potential energy.
The equilibrium solutions of(5.1) are the roots of the equation

2 , 3 R 2 W(x) = p. p 1 (5.2) (5.3) (5.4)
Results similar to Proposition (2.2), (2.3), (2.4) can be obtained by exactly the same method under suitable assumptions on W'(x 3 ). In (5.4)p is measured per unit of thickness.

Next we show that in the neighbourhood of an equilibrium solution for the membrane we can fmd equilibrium solutions for sufficiently thin shells, and also that dynamic solutions for the membrane are the uniform limits of solutions for the thick shell as 6 tends to zero.

Proposition (S.l)

We suppose that: (i) x 0 > 1 is a solution of (5.4) for a given p > 0, (ii) W"(x~) =I= 0, (iii) Hypotheses (A) and (B) hold for W' (x 3 ).

Then, there exists an interval (0, 6 I) and a function x = x(6) E C 2 (0, 6 d the values of which give equilibrium solutions for a shell with thickness parameter 6 corresponding to an external pressure P = 6p.

Proof.

We are going to apply the implicit function theorem to

h(y, 6) = g(~, 6 ) 0 Lety =x 3 • From h( 0) = lim g(y, 6 ) = ~ W'( ) y, 6-+o 6 pRj y '
(5.5) and with (i), it follows that h(x~, 0)-p = 0.

It is easy to check that d 2 , lim -h(y, 6) = -2 yW (y).

6-.o dy pR 1

In particular,

d 3 . d 2 3 , 3) dx 3 h(x 0 ,0} = hm -h(y,6)i:v=x3 = R 2 x 0 W (x 0 =I= 0. 6-.o dy • p 1
By the implict function theorem, there exists a 6 1 > 0 and a X E C 2 [0, 61) such that x(O) = x 0 , x(6) > 1 and

(5.6)

(5.7)

Hence, x =x( 6) is an equilibrium solution for a shell of thickness parameter 6 andP = p6.

Remark. The interval (0, 6 d can be extended up to a 6max• If W"(x~) = 0, the previous result for this x 0 cannot be guaranteed. where x(t) is the solution of (5.2) with x(O} = x 0 and energy E.

Proof.

Let x 0 > 0. We consider t 1 < T such that E-V(x(tt)) = 0. Then£-V(x) >O fort< t 1 •

The result for x(t), tE [0, tt],x(t, 6), tE [0, t 1 -e] e>O follows by considering (5.2}, {5.3}, (1.34}, (1.36} and also that 3 lim ~ V(x, 6) = V(x).

6-.o u

For t > t 1 , we consider the minus sign in the square root of (1.34) and (5.2). If there exists t 2 such that E-V(x(t 2 ))=0, then the result follows by the same arguments for x(t), t E [t 1 , t 2 ] and x(t, 8), t E [t 1 + e, t 2 -e]. We proceed similarly up tot= T. This also holds for x 0 < 0. o Remark. It follows from (5.4) that a necessary condition for the existence of more than one equilibrium solution for appropriate values of p is that W"(x 3 ) = 0

(5.9) for at least one x > 1.

In terms of a(>..) the above condition says that (2

) N = 1, M = 1 Q = 1, p. = i (3 = 1, v = 1
There is one root of (5.10) in (0, 1).

Remark: For this class of materials hypothesis (D2) cannot be satisfied.

Section 6

The differential equation describing the radial motion of the cylindrical shell is

• I

The hypotheses for f(x 2 , 5) are the same as (1.29) except that x 2 replaces x 3 .

The phase•plane analysis yields similar results to the spherical shell with the exception that under conditions (C2) and (Dl) for f(x 2 , 5) and the other usual hypotheses, the solutions exist for all time and are of the following type:

( 1) Periodic solutions under the same kind of initial data as in the spherical shell.

(2) Unbounded solutions, with exponential growth for large values of t, under the analogous class of initial data that give blow-up for the spherical shell.

The remaining results are similar with the number 2 in this case playing the same role as 3 for the spherical shell.

p

  = p(x, t)is an arbitrary hydrostatic pressure.In particular, forB given by (1.4), (1.8) becomes

( 1 )

 1 If (C1) and (D1) hold, then (a) n = {y,, ... ,Y:zr},r~ 1, (b) x 0 E [ y 28 -1 , y 28 ] , for some s . ;;;;;. r, (c) V(x,o).r;;;.EforallxE [Y:z 8 -1 ,Y:z.]. (2) If (C2) and (Dl) hold, then (a) n = {y4, ... ,Y:zr-d,r~ 1, (b)xo ~ y,.

  3). 1(b) We prove first thatx 0 E fY1 ,Y:zrl. Let us assume that x 0 fl. [y 1 , Y:zrl, then V(x 0 , o) > E, which is impossible by the energy equation. In particular x 0 E [Y:za-1 ,y 28 ] for a integers.;;;;;, r, which is determined by V( •, o), (x 0 ,x 0 ) and E. 1 (c) follows by similar arguments. 0 Remark (2.1) (A) If (C3) and (D 1) or (D3) hold, then if P ~ D the results (2) hold, while if P.;;;;;. C then the results (1) hold. (B) If (C2) and (D2) hold, there exist some values of E such that n = 1/>.

  t-+tmaxMoreover, if x 0 is sufficiently small, we have the following estimate for tmax tmax = !(£-Vor 1 ' 2 xg 12 + e(xo), where and M, N, {3 are given in hypotheses {02).

  The first case in (b) and (c) and (d.2) correspond to Proposition {2.1). The second case in (b) and (c) and (d.2) correspond to case (a) in the present Proposition. oRemark (2.3)We suppose that (C2) and (D2) hold, and (ii) and (iii) are as in (2.4). Then, any (xo, x 0 ) E!» corresponds to one of the following situations:

  differentiate equation (1.27) with respect to x 3 • Using (3.I), (3.2), we obtain the finite difference equation w(y)-w(y-2') = 'YeYg'(I +e:ll),w(-oo) = f.J..

  r wN(Y + r£1 = wN(Y) + ' Y L .e>'+i£g'(I + eY•I£) j=t and r-1 wN(y-r.Y") = wN(y)-'Y L e>'-i£g'(I + ey-Jt), j=O where r is a postitive integer.

  z-;;:,.-N+!?' and i=l I r = y [z + N] , (the brackets mean the integer part) then-N.;;;;z -r Sf'.;;;;-N + !?.

  with (3.4), {3.7) becomes r WN(z) = IJ. + ' Y L ez+(j-r)tg'(1 + ez+(j-r)t), j=l z";;t:-N+Y.

  ) = JJ.-'Y L ez+.t(r-i>g'(l +ez+.t(r-J>), z<:.-N i=O and r = [-~(N+z)+1].

K

  = K(x 3 ,N) = Suplg"(~)l {3.14) ~E [1, (;~~;N] and K = sup K(x 3 , N) = K(N). n ~N,

  3 ) = W(F) is convex and non-decreasing in each variable. Therefore, a(o:A.+(1-o:)s.t) = a(o:A.+(1-o:)s.t, (o:A.+(l-o:)s.tF112, (o:s.t+(1-o:)s.tF1'2)

(h 2

 2 h x1 x 1 = h(x!)h (xi) = --(h (xi)(x!)] = ..2"(xd + r(xi)h 2 (xi).

  (4.4) We Iety(xd = h 2 (x!). Then, we can write the previous equation as d dx/(xd == r(xt)y(xt) = ..2"(xt)

  ' [22]. (5.1) This equation is obtained by approximating (1.24) up to the first order in 6 [(1.22)]. For the thin shell, 6 is proportional to the thickness e of the shell: 6 = 3 (e/R d.

  least one solution for >.. < 1. It is interesting to note that the assumption of convexity of a(>..) does not contradict (5.10).We are going to give examples for a class of materials studied byOgden [20] satisfying a" > 0. Moreover, by choosing the appropriate growth conditions we can obtain two of the cases studied in Section 2.The materials studied by Ogden have the following isotropic stored energy functionN a(>..t, >..2, >..3) = L o: 1 (X~i + >..~1 + M'-3) i•l M + L 131 ((>.., >..2) 11 1 +(At >..3)"1 + (>..2>..3)"1-3). /=1For the radial deformation under study, M a(>..) = a(>..,>.. -112 , >.. -112 ) = L o: 1 (X~-'i + 2>.. -~-'il2 -3) 1•1 N + L 131 (2>.."112 + A -IIJ-3).

  Hypotheses (Cl) and (Dl) for W'(~):(I) N = 1, M = 0 0: = 1' ll = 3/2A root of (5.10) in (0, 1) is>..= 1/256 119 •(2) N = 3, M = 0 O:t = 0.153, 0:2 = -0.092, 0:3 = 0.5 Ill = 2.8, ll2 = 2.5, ll3 = 2.0 Two of the roots of (5.10) for>..< 1 are >.. 1 = 0.5, >.. 2 = 0.45.Hypotheses (Dl) and (C2) for W'(~):(1) N = 1, M = 0 Q = 1, p. = !A root of (5.10) in (0, I) is X= 1/256 119 •

x

  log ( 1 + : 2 ) x + [log (1 + : 2 ) -5 : x 2 ] x 2 + [(x 2 , 5) = P, where 1 JX 2 W'(~) dl: f(x2' 5) = R2 x>+6 1:-1 '>• p I 1+6 'i and R~ 5 = Rj -1 >O. The energy equation is !x 2 log (1 + : 2 )x 2 +F(x 2 ,5)-!Px 2 = E where •x> F(x 2 , 5) = ! I [(X, 5)dA.

  Now we are going to state the results of the phase-plane analysis of the initial-boundaryvalue problem. We fix 8 > 0 and P > 0. Several examples are shown at the end to illustrate some of the different cases that may occur.

	Proposition 2.2
	We suppose that
	(i) (C1) and (Dl) hold,
	(ii) g(x 3 , 8) has 2n critical points,
	(iii) the roots of g(x 3 , 8) =Pare
	(2.1)
	(2.2)

•)is the set c+(x( •))given by c+ = Ux(t). tE [0, tmax) {.X 1 EIR+, 1...:i...:2k+ I, k...:n}.

Acknowledgment

I would like to thank Dr. J. M. Ball for his help and suggestions in the writing of this paper.

List of symbols

8, e 00 J,j

2'

T, t, t


Tapan Kumar Hota

Satyajit Pramanik

Manoranjan Mishra

A General Approach to the Linear Stability Analysis of Miscible Viscous Fingering in Porous Media

We analyze the linear growth of the viscous fingering instability for miscible, non-reactive, neutral buoyant fluids using the nonmodal analysis (NMA). The onset of instability is obscured due to the continually changing base state, and the normal mode analysis is not applicable to the non-autonomous linearized perturbed equations. Commonly used techniques such as frozen time method or amplification theory approach with random initial condition using transient amplifications yield substantially different results for the threshold of instability. We present the classical non-modal methods in the short-time limit using singular value decomposition of the propagator matrix. Using the non-modal approach we characterize the existence of a transition region between a domain exhibiting strong convection and a domain where initial perturbations are damped due to diffusion. Further, at the early times the algebraic growth of perturbations is possible which suggest that NMA could play an important role in describing the onset of instability in the physical phenomenon involving VF.

Introduction

Flow of miscible fluids in a porous medium is found in a wide range of natural processes as well as industrial applications [START_REF] Homsy | Viscous fingering in porous media[END_REF]. These displacement flows feature an interfacial instability that occurs when a less viscous fluid invade into the more viscous one. Finger like patterns are resulted at the interface, hence this is called viscous fingering (VF) instability. This hydrodynamic instability which is also known as Saffman-Taylor instability [START_REF] Saffman | The penetration of a fluid into a medium or Hele-Shaw cell containing a more viscous liquid[END_REF] has been a subject of many theoretical and experimental studies [START_REF] De Wit | Viscous fingering of miscible slices[END_REF][START_REF] Mishra | Differences in miscible viscous fingering of finite width slices with positive or negative log mobility ratio[END_REF][START_REF] Nagatsu | Hydrodynamic Fingering Instability Induced by a Precipitation Reaction[END_REF][START_REF] Rana | Combined influences of viscous fingering and solvent effect on the distribution of adsorbed solutes in porous media[END_REF] for many years. Stability analysis of VF has been performed extensively over the years [START_REF] Homsy | Viscous fingering in porous media[END_REF][START_REF] Ben | A spectral theory for small amplitude miscible fingering[END_REF][START_REF] Kim | Linear stability analysis on the onset of the viscous fingering of a miscible slice in a porous media[END_REF][START_REF] Pramanik | Linear stability analysis of Korteweg stresses effect on the miscible viscous fingering in porous media[END_REF]. Traditionally, the modal analysis along with quasi-steady-state approximation (QSSA) has been used in different coordinate systems [START_REF] Kim | Linear stability analysis on the onset of the viscous fingering of a miscible slice in a porous media[END_REF][START_REF] Pramanik | Linear stability analysis of Korteweg stresses effect on the miscible viscous fingering in porous media[END_REF]. Although, at the early times QSSA in a self-similar coordinate system successfully predicts the unconditional stability of perturbations, it fails to describe the early period dynamics of the disturbances [START_REF] Hota | Non-modal linear stability analysis of miscible viscous fingering in a Hele-Shaw cell[END_REF]. Another approach which is more practical is the amplification theory (AT), in which the linear equations are solved as an initial value problem. However, this approach has a drawback of choosing representative initial condition. In this article, using non modal stability theory, we unify the frozen time method and AT to explore two important aspects of VF instability, namely, the effect of nonnormality of the governing non-autonomous linear operators and the optimal amplifications of the disturbances. The paper is organised as follows. In Sect. 2, the mathematical formulation of the physical model is presented. Section 3 describes the non-modal stability analysis and the numerical method to calculate the energy growth function G(t). In Sect. 4 we present the results and discussions. Finally, the conclusions are given in Sect. 5. Fig. 1: Schematic of the flow configuration with coordinate system. Initially the interface is flat (dotted line) and then a wave like infinitesimal perturbation is applied.

Mathematical Formulations

Consider a uniform rectilinear displacement of miscible fluids in a two-dimensional, homogeneous, porous medium with constant permeability. The fluids are assumed to be incompressible, miscible, non-reactive, neutrally buoyant and the dispersion is isotropic. The non-dimensional governing equations for the prescribed two-dimensional flow in a reference frame moving with velocity U are given by the following coupled nonlinear partial differential equations (PDEs) [START_REF] Tan | Stability of miscible displacements in porous media: Rectilinear flow[END_REF][START_REF] Homsy | Viscous fingering in porous media[END_REF],

∇ • u = 0, ∇p = -µ(c)(u + i), ∂c ∂t + u • ∇c = ∇ 2 c, (1)
where p is the dynamic pressure, u = (u, v) is the Darcy's velocity, c is the concentration of the solvent, i is the unit vector along the x (downstream) direction, and µ(c) is the fluid viscosity which depends exponentially on the solute concentration, i.e., µ(c) = exp(Rc) [START_REF] Tan | Stability of miscible displacements in porous media: Rectilinear flow[END_REF]. Here R = ln(µ 2 /µ 1 ) is the log-mobility ratio, where µ 1 and µ 2 correspond to the viscosity of the less and more viscous fluid, respectively. The coupled PDEs (1) are provided with the following initial and boundary conditions in the moving frame of reference [START_REF] Homsy | Viscous fingering in porous media[END_REF][START_REF] Pramanik | Linear stability analysis of Korteweg stresses effect on the miscible viscous fingering in porous media[END_REF]: Initial conditions:

u = (0, 0), c(x, y, t = 0) = 0, x ≤ 0 1, x > 0, ∀y. (2)
Boundary conditions:

u = (0, 0), ∂c ∂x = 0, |x| → ∞, stream-wise direction (3)
u is arbitrary,

∂c ∂y = 0 = ∂v ∂y , ∀x, span-wise direction (4)
The base state of the flow is assumed to be pure diffusion of the concentration along the axial direction that can be written as

u b = (u b , v b ) = 0, c b (ξ) = 1 2 1 + erf ξ 2 , µ b = µ b (ξ), (5)
where erf(z) = 2 √ π z 0 e -η 2 dη is the error function and ξ := x/ √ t is the similarity variable transformation. Introducing the infinitesimal perturbations in terms of the Fourier modes of the form (u , c )(ξ, y, t) = (u , c )(ξ, t)e iky and using standard procedure of linear stability analysis along with central difference scheme to discretize the spatial variable, the nonlinear coupled PDEs (1) can be written as [START_REF] Pramanik | Linear stability analysis of Korteweg stresses effect on the miscible viscous fingering in porous media[END_REF][START_REF] Hota | Non-modal linear stability analysis of miscible viscous fingering in a Hele-Shaw cell[END_REF]

dc dt = A(k, R, t)c , (6)
where A = A(k, R, t) is the stability matrix.

For such non-autonomous system (Equation ( 6)), there are two distinct regions in the evolution of disturbances [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF]. The first is the asymptotic long time region, and the second is the transient region. In mathematical sense, instability behaviors in the the first region are determined by the eigenvalues of the stability matrix A. In the second region, the transient growth may be quite substantial, such that the nonlinear region may be reached before the growth of eigenvalue mode. Thus, the eigenvalue approach is insufficient to study the stability analysis. Non-normal, autonomous systems, i.e., where the stability matrix A is not depending on time and the eigenvalues are not orthogonal, have been explored by many researchers (see [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF] and the references within). Recently, few works have discussed non-modal stability analysis with unsteady base-state flow [START_REF] Rapaka | Non-modal growth of perturbations in density-driven convection in porous media[END_REF][START_REF] Doumenc | Transient RayleighBnard-Marangoni convection due to evaporation: a linear non-normal stability analysis[END_REF][START_REF] Daniel | Optimal perturbations of gravitationally unstable, transient boundary layers in porous media[END_REF][START_REF] Hota | Non-modal linear stability analysis of miscible viscous fingering in a Hele-Shaw cell[END_REF]. But in the case of VF, many aspects of stability matrix A(k, t) remained unexplored compared to its autonomous counterpart.

Transient Behaviors and Non-modal Analysis

In the framework of NMA, there are two major and broader aspects which are of ample interest, i.e., the responses to external excitations and the transient energy growth of initial conditions. Mathematically, the former can be studied from the structure of " -pseudospectra" [START_REF] Trefethen | Spectra and Pseudospectra: the Behavior of Nonnormal Matrices and Operators[END_REF] which is given by Λ (A) = {z ∈ C : σ min (z -A) ≤ }, where σ min (A) denotes the smallest singular value of A, and 0 < 1; and the latter can be analysed from the energy growth function G(t) that identifies the optimal growth of energy at time t. In addition, it is often useful in hydrodynamic stability problem to know the initial growth rate of the energy growth function. This can be obtained from numerical abscissa of the stability matrix A [START_REF] Trefethen | Spectra and Pseudospectra: the Behavior of Nonnormal Matrices and Operators[END_REF]. Following [START_REF] Rapaka | Non-modal growth of perturbations in density-driven convection in porous media[END_REF][START_REF] Hota | Non-modal linear stability analysis of miscible viscous fingering in a Hele-Shaw cell[END_REF]] G(t) can be described as

G(t) = G(t, k, R) := max c 0 Φ(t 0 ; t)c 0 2 = Φ(t 0 ; t) = sup j s j (t), (7)
where Φ(t 0 ; t) is the propagator matrix or matrizant satisfying the matrix differential equation d dt

Φ(t 0 ; t) = A(k, R, t)Φ(t 0 ; t), (8)
and it is related to the concentration perturbation field by c (t) = Φ(t 0 ; t)c 0 with c (t 0 ) = c 0 being an arbitrary function. Here, s j 's are the singular values of Φ(t 0 ; t) and • 2 is the standard Euclidean norm. It is important to note that the pseudospectra and the growth function are dependent on the definition of norm. From mathematical and physical considerations, an appropriate measure of the disturbance is indispensable. In this article the standard Euclidean norm • 2 and associated inner product have been used. Further, as the velocity perturbation is slaved to the concentration perturbation, the optimal amplification G(t) has been calculated only for the concentration perturbation c [START_REF] Doumenc | Transient RayleighBnard-Marangoni convection due to evaporation: a linear non-normal stability analysis[END_REF]. For stability analysis we consider following quantities: the spectral abscissa and numerical abscissa are given by α(A) ≡ max{ (λ(A)}, η(A) ≡ max{λ(A + A T )/2}, respectively, where λ(A) denotes the eigenspectrum of A, denotes the real part, and A T denotes the transpose of the matrix A. In the present paper, the computational domain has been chosen to be [-50, 50] with step size 0.2. The initial value problem ( 8) is solved by Runge-Kutta fourth order method, and a matlab GUI EigTool [START_REF] Trefethen | Spectra and Pseudospectra: the Behavior of Nonnormal Matrices and Operators[END_REF] has been used to draw the pseudo-spectra of the stability matrix A. The detail of the numerical procedure can be found in [START_REF] Hota | Non-modal linear stability analysis of miscible viscous fingering in a Hele-Shaw cell[END_REF].

Results and Discussion

In stability analysis, the eigenspectrum α(A) is the principal aid to give an insight into how a system behaves. If the stability matrix A is non-normal (i.e., AA T = A T A), then the pseudospectra Λ (A) are likely to explain the system behavior better than the eigenvalues λ(A). Thus, Λ (A) will help us to understand the response to the external excitations in the parameter space k and R. Figure 2 shows the parabolic profile (the dashed line) of the numerical range, W (A) = {x T Ax : x = 1}, of the stability matrix A. It is clearly visible that the boundary of W (A) strictly contains the spectrum of A. This reflects that A can not be unitarily diagonalizable; or in other words, their eigenfunctions are not orthogonal. In the inset figures (Fig. 2), it is illustrated that the pseudospectra protrudes strongly into the right half-plane, which imply that the evolution process will be susceptible to large transient effects. This is due to the nonnormality of the stability matrix, and it signifies that the system is unstable, which is not captured by analysing the spectrum alone. This implies that at the initial time the energy of the disturbance is growing faster than what is anticipated from the spectral abscissa α(A). It can be conjectured that with increasing R, the boundary of numerical range becomes wider and the eigenvalues are more sensitive to the external forcing, as confirmed from the contours of pseduspectra. Next, we present the optimal amplification of the perturbations to understand the transient effects. For this purpose we plot the energy growth function G(t) for (R, k) ∈ {(3, 0.225), (4, 0.3), (5, 0.35)} in Fig. 3(a). This figure shows that at the initial period diffusion is dominated and each curve experiences a substantial energy growth. The onset of instability is the time when G(t) starts increasing (the first local minimum) and it is shown as black dot (•). Thus, NMA clearly distinguishes the domain where initial perturbations are damped or have no time to grow significantly due to diffusion, and a domain exhibiting strong convection. Further, to understand how the VF mechanism responses to external forcing, the spatial structures of the optimal input can be obtained from the singular value decomposition of the propagator matrix Φ(t 0 ; t) [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF][START_REF] Rapaka | Non-modal growth of perturbations in density-driven convection in porous media[END_REF][START_REF] Hota | Non-modal linear stability analysis of miscible viscous fingering in a Hele-Shaw cell[END_REF]. Using matlab curve fitting tool (cftool), with 95 % confidence bounds, we found empirical relationships between the onset of instability, t on , and R, as well as the critical wave number, k c (the minimum wave number at which the instability sets in), and R. The obtained results are shown in Fig. 3(b). It is shown that t on ∼ R -1.9 , i.e., nearly a inverse square of R, whereas k c ∼ R 0.82 This empirical relationships are in good agreement with Tan and Homsy [START_REF] Tan | Stability of miscible displacements in porous media: Rectilinear flow[END_REF], who measured t on ∼ R -2 , and k c ∼ R, for a step-like initial concentration profile. The non-normality of the stability matrix A could be the source of observed difference between the two cases.

Conclusion

A unified approach of linear stability analysis for the VF is presented. It is shown that at early times diffusion dominates, which causes the energy decay before it starts amplifying due to strong convection. For this purpose we have studied different components of the spectrum, such as: the spectral abscissa, boundary of numerical range,and -pseudospectra. To understand the transient growth of energy, we studied the optimal amplification of disturbances for various flow parameters. It can be concluded that, such approach not only presents a comprehensive stability analysis algorithm, but also explains the physical mechanism appropriately. It will be very interesting to apply this procedure for other unsteady base state problems such as VF in liquid chromatographic condition [START_REF] Rana | Combined influences of viscous fingering and solvent effect on the distribution of adsorbed solutes in porous media[END_REF], analysing the effect of precipitation reactions in CO 2 sequestration techniques [START_REF] Nagatsu | Hydrodynamic Fingering Instability Induced by a Precipitation Reaction[END_REF], or in understanding the effect of external forces, e.g., magnetic field [START_REF] Mishra | Influence of a simple magnetic bar on buoyancydriven fingering of traveling autocatalytic reaction fronts[END_REF] to VF to name a few.

Scientific Validation

This paper has been unanimously validated in a collaborative review mode with the following reviewers:

-Agota Toth, from University of Szeged -Denis Grebenkov -Soumya Banerjee, from Broad Institute of MIT and Harvard, Ronin Institute

Fig. 2 :

 2 Fig. 2: Pseudospectra of A for (a) R = 1, k = 0.08, t = 31, (b) R = 4, k = 0.3, t = 2.5. Black dots (•): eigenvalues; Black square ( ): numerical abscissa; Dashed line: boundary of the numerical range; Solid lines: contours from innermost to outermost representing levels from = 10 -2.5 to 10 -0.5 with increment 10 -0.5 . λ i , λ r are the imaginary and real part of eigenvalues, respectively. Inset figure shows the numerical abscissa, η(A) ( ) and spectral abscissa, α(A)(•).

Fig. 3 :

 3 Fig. 3: (a) Optimal amplification G(t) for different R and the corresponding critical wave numbers. The black dot denotes the onset.(b) The variation of the onset of instability t on and the critical wave number k c (inset), for different values of R determined from NMA. The circles (•) are simulation data and continuous lines are fitted curves.

Acknowledgment

The authors would like to thank the three reviewers for their helpful critics and suggestions. This helped us in improving the standard of the manuscript. S.P. gratefully acknowledges the National Board for Higher Mathematics, Department of Atomic Energy, Government of India for the financial support through a Ph.D. fellowship.