Geophysical Ice Flows: Analytical and Numerical Approaches Will - PowerPoint PPT Presentation

Geophysical Ice Flows: Analytical and Numerical Approaches Will Mitchell University of Alaska - Fairbanks July 23, 2012 Supported by NASA grant NNX09AJ38G

Ice: an awesome problem ...velocity, pressure, temperature, free surface all evolve

Outline ◮ I. Introduction to viscous fluids ◮ II. Exact solutions ◮ III. Finite element solutions

Stress: force per unit area A tornado sucks up a penny. At any time: ◮ The fluid into which n points exerts a force on the penny ◮ Force / area = stress ◮ The stress vector is a linear function of n ◮ In a Cartesian system: stress = σ · n ◮ σ is the Cauchy stress tensor Quiz: Suppose there is no p ≥ 0 such that σ · n = − p n . Physical interpretation?

Decomposition of stress ◮ In a fluid at rest, σ · n = − p n , so σ = − pI ◮ In general, choose p = − Trace( σ ) / d , so σ = − pI + τ where τ has zero trace. ... this defines pressure p and deviatoric stress τ

Strain rate ◮ let u be a velocity field ◮ the gradient of a vector is the tensor ( ∇ u ) ij = ∂ u j ∂ x i ◮ define Du = 1 2( ∇ u + ∇ u T ) 2 ∂ u 1 ∂ u 2 + ∂ u 1   ∂ x 1 ∂ x 1 ∂ x 2 ◮ in 2D: Du = 1   ∂ u 1 + ∂ u 2 2 ∂ u 2 2   ∂ x 2 ∂ x 1 ∂ x 2 ◮ Du is the strain rate tensor. True or False: “Since Du is a derivative of velocity, it measures acceleration.”

Constitutive Laws: Newtonian How does a fluid respond to a given stress? ◮ For Newtonian fluids ( e.g. water) a linear law: τ = 2 µ Du The proportionality constant µ is the viscosity .

Constitutive Laws: Glen’s For glacier ice, a nonlinear law. ◮ define � � 1 1 � τ � = 2Tr( τ T τ ) and � Du � = 2Tr( Du T Du ) ◮ assume � Du � = A � τ � n ◮ the law is either of τ = ( A � τ � n − 1 ) − 1 Du τ = A − 1 / n � Du � (1 − n ) / n Du A is the ice softness , n ≈ 3 is Glen’s exponent.

Stokes Equation What forces act on a blob occupying a region Ω within a fluid? ◮ body force, gravity: � Ω ρ g ◮ force exerted by surrounding fluid: � � ∂ Ω σ · n = Ω ∇ · σ Force = rate of change of momentum ρ g + ∇ · σ = ∂ � � ρ u ∂ t Ω Ω D � = Dt ρ u . Ω � � D � < 10 − 15 so � � � � In glaciers, Fr = Dt ρ u � : � p � � � ρ g + ∇ · σ = 0 .

Incompressible Stokes System (two versions) � ρ g + ∇ · σ = 0 ∇ · u = 0 ∇ · σ = ∇ · τ + ∇ p = 2 µ ∇ · ˙ ǫ − ∇ p = µ ∇ · ( ∇ u ) + µ ∇ · ( ∇ u T ) − ∇ p ∇ · ( ∇ u ) = ∂ ∂ u j = ∂ ( ∇ · u ) = 0 ∂ x j ∂ x i ∂ x i ∇ · ( ∇ u T ) = ∂ ∂ u i = ∆ u ∂ x j ∂ x j � − µ △ u + ∇ p = ρ g ∇ · u = 0

The Biharmonic Equation ◮ for 2D, incompressible flow: u = ( u , 0 , w ) and ∇ · u = 0 ◮ there is a streamfunction ψ such that ψ z = u , − ψ x = w . ◮ take the curl of the Stokes eqn � � ∇ × − µ △ u + ∇ p = ρ g to get the biharmonic equation ψ xxxx + 2 ψ xxzz + ψ zzzz = 0 or △△ ψ = 0 . Quiz: give an example of a function solving the biharmonic eqn.

Ice: an awesome problem ...velocity, pressure, temperature, free surface all evolve

Slab-on-a-slope: a tractable problem � � g = ( g 1 , g 2 ) = ρ g sin( α ) , − cos( α ) ...no evolution

Stokes bvp find a velocity u = ( u , w ) and pressure p such that −∇ p + µ △ u = − g on Ω ∇ · u = 0 on Ω u (0 , z ) − u ( L , z ) = 0 for all z u x (0 , z ) − u x ( L , z ) = 0 for all z u = f on { z = 0 } w = 0 on { z = 0 } w x + u z = 0 on { z = H } 2 w zx − ( u xx + u zz ) = g 1 /µ on { z = H } where f ( x ) = a 0 + � ∞ n =1 a n sin( λ n x ) + b n cos( λ n x ).

biharmonic bvp find a streamfunction ψ such that △△ ψ = 0 on Ω ψ z (0 , z ) − ψ z ( L , z ) = 0 for all z ψ xz (0 , z ) − ψ xz ( L , z ) = 0 for all z ψ x (0 , z ) − ψ x ( L , z ) = 0 for all z ψ xx (0 , z ) − ψ xx ( L , z ) = 0 for all z ψ ( x , 0) = 0 for all x ψ z ( x , 0) = f for all x ψ zz ( x , H ) − ψ xx ( x , H ) = 0 for all x 3 ψ xxz ( x , H ) + ψ zzz ( x , H ) = − g 1 /µ for all x . (1)

biharmonic bvp, subproblem: f = 0 u ( x , z ) = g 1 H µ z − g 1 2 µ z 2 ψ ( x , z ) = g 1 H 2 µ z 2 − g 1 6 µ z 3 − → w ( x , z ) = 0 ...this is Newtonian laminar flow , a well known solution.

biharmonic bvp, subproblem: f � = 0 strategy: ◮ separate variables: ψ ( x , z ) = X ( x ) Z ( z ) ◮ periodicity: take X ( x ) = sin( λ x ) + cos( λ x ) for λ = 2 π n L ◮ the biharmonic eqn reduces to an ODE: � λ 4 Z − 2 λ 2 Z ′′ + Z (iv) � 0 = △ 2 ( XZ ) = X ◮ for λ > 0 this gives Z ( z ) = a sinh( λ z ) + b cosh( λ z ) + cz sinh( λ z ) + dz cosh( λ z ) ◮ homogeneous bcs determine b , c , d in terms of a ◮ weighted sum gets the nonzero condition

Exact Solutions Horizontal Component of Velocity: ∞ λ n H 2 ( a n sin( λ n x ) + b n cos( λ n x )) u ( x , z ) = a 0 + g 1 H µ z − g 1 2 µ z 2 + � Z ′ n ( z ) n H 2 + cosh 2 ( λ n H ) λ 2 n =1 where n ( z ) = − 1 � � Z ′ H cosh( λ n H ) sinh( λ n ( z − H )) + λ n z cosh( λ n ( z − H )) + cosh( λ n H ) − λ n H sinh( λ n H ) � · cosh( λ n ( z − H )) λ n H 2 � + λ n z sinh( λ n ( z − H )) + λ n cosh( λ n z ) .

Exact Solutions Vertical Component of Velocity: ∞ λ 2 n H 2 � w ( x , z ) = n H 2 + cosh 2 ( λ n H )( b n sin( λ n x ) − a n cos( λ n x )) Z n ( z ) λ 2 n =1 where Z n ( z ) = sinh( λ n z ) − 1 H cosh( λ n H ) z sinh( λ n ( z − H )) � cosh( λ n H ) − sinh( λ n H ) � + z cosh( λ n ( z − H )) λ n H 2 H

Exact Solutions Pressure: p ( x , z ) = g 2 z − g 2 H ∞ λ 3 n H ( a n cos( λ n x ) − b n sin( λ n x )) � + 2 µ × n H 2 + cosh 2 ( λ n H ) λ 2 n =1 � sinh( λ n z ) − cosh ( λ n H ) � cosh ( λ n ( z − H )) . λ n H ... this is new.

The finite element method ◮ numerical approximation of p and u ◮ requires a mesh of the domain: ◮ leads to a system of linear equations A x = b

Variational Formulation: incompressibility Incompressibility: ∇ · u = 0. We seek a u ∈ H 1 (Ω) such that for all q ∈ L 2 (Ω), we have � q ∇ · u = 0 . Ω u is a trial function ; q is a test function .

Variational Formulation: Stokes Put σ = τ − pI in the Stokes equation: 0 = ∇ · τ − ∇ p + ρ g Dot with v ∈ H 1 and integrate over Ω. Integration by parts gives � � � � τ : ∇ v − p ∇ · v − n · σ · v = ρ g · v . Ω Ω ∂ Ω Ω More manipulation gives 1 � � � � ∇ u T + ∇ u � � ∇ v + ∇ v T � 2 µ : − p ∇ · v = ρ g · v . Ω Ω Ω

Variational Formulation E × L 2 such that for all ( v , q ) ∈ H 1 E 0 × L 2 we have Find ( u , p ) ∈ H 1 1 � � � � � ∇ u T + ∇ u � � ∇ v + ∇ v T � 2 µ : − p ∇· v + q ∇· u = ρ g · v . Ω Ω Ω Ω Still a continuous problem: ( u , p ) satisfy many conditions. Make a discrete problem using finite-dimensional spaces.

Pressure Approximation space Continuous functions that are linear on each triangle: a 16-dimensional space with a convenient basis.

Velocity Approximation space Continuous functions that are quadratic on each triangle: a 49-dimensional space (per component) with a convenient basis.

Implementation I ...based on [Jar08] but with an important difference 1 from d o l f i n import ∗ 2 #Set domain parameters and p h y s i c a l c ons ta nts 3 Le , He = 4e3 , 5e2 #length , h e i g h t (m) 4 alpha = 1 ∗ p i /180 #s l o p e angle ( r a d i a n s ) 5 rho , g = 917 , 9.81 #d e n s i t y ( kg m − 3) , g r a v i t y (m sec − 2) 6 mu = 1e14 #v i s c o s i t y (Pa sec ) 7 G = Constant (( s i n ( alpha ) ∗ g ∗ rho , − cos ( alpha ) ∗ g ∗ rho ) ) 8 #Define a mesh and some f u n c t i o n spaces 9 mesh = Rectangle (0 ,0 , Le , He , 3 , 3 ) 10 V = VectorFunctionSpace ( mesh , ”CG” , 2) #pw q u a d r a t i c 11 Q = FunctionSpace ( mesh , ”CG” , 1) #pw l i n e a r 12 W = V ∗ Q #product space 13 ””” Define the D i r i c h l e t c o n d i t i o n at the base ””” 14 def LowerBoundary ( x , on boundary ) : r e t u r n x [ 1 ] < DOLFIN EPS and on boundary 15 16 S l i p R a t e = E x p r e s s i o n (( ”(3+1.7 ∗ s i n (2 ∗ p i/%s ∗ x [ 0 ] ) ) \ /31557686.4 ”%Le , ” 0.0 ” ) ) 17 18 bcD = Di ri c hl et B C (W. sub (0) , SlipRate , LowerBoundary )