Transcript Course Notes 1

The Australian Computational Earth
Systems Simulator
(ACcESS)
Multi-Scale Behaviour in the GeoScience I: The Onset of Convection
and Interfacial Instabilities
by
Hans Mühlhaus
Overview
What is Geodynamics
Mantle Convection, Spreading, Folding, Landscape
Evolution, Earthquakes, Volcanoes
The onset of Natural Convection:
Linear instability analysis
Numerical Simulations
Nusselt Number
Remarks on Weakly Non-Linear
Analysis
Galerkin methods
Earth dynamics
continental
crust
convergent
plate boundary
transform
plate boundary
divergent
plate boundary
oceanic
spreading ridge
trench
trench
convergent
plate boundary
strato volcano
shield volcano
oceanic
crust
island arc
hot spot
subducting plate
lithosphere
asthenosphere
TERRA MESH
Earth dynamics
BENEFITS
continental
crust
convergent
plate boundary
transform
plate boundary
divergent
plate boundary
oceanic
spreading ridge
trench
trench
convergent
plate boundary
strato volcano
shield volcano
oceanic
crust
island arc
hot spot
subducting plate
lithosphere
asthenosphere
2nd talk: Volcano modelling
Montserrat, West Indies
3rd talk: Particle Processes
–Spherical particles
–Selection of contact physics:
–Non-rotational and rotational dynamics
–Friction interactions
–Linear elastic interactions
–Bonded interactions
4th talk: modeling of geological folds
Director evolution
n i  Wijn n j
Wijn  Wij  ( Dki kj  Dkj ki )
 ij  nin j
n : the director of the anisotropy
W, Wn: spin and director spin
D, D’: stretching and its deviatoric part
Governing Equations
Navier Stokes Equations:
v
 (  v  v)  g  p    (v  (v)T )
t
 0 
g   
 g 
Heat Equation:
T
c p (  v  T )  k 2T
t
d
 v
dt
v(t,x) is the velocity vector, T is the temperature,  (4000kg/m2) is the
density,  (1021Pas) is the viscosity, cp (10-3WK-1s/kg) is the specific heat
at constant pressure and k (4Wm-1K-1) is the thermal conductivity
Governing Equations, cont.
Temperature dependence of density:
   0 (1   p (T  T0 ))
p ( 3 10-5K-1 ) is the thermal expansion
coefficient, T0 (288K) is the surface
temperature
Simplified convection model:
T=T0 , v2=0
T,1 , v1=0
x2
H
L
x1
T=T1 , v2=0
T,1 , v1=0
Consider a square planet……
Governing Equations, cont.
Nondimensionalisation:
x  H~
x,
H 2 0c p ~
t
t,
k



v
H~
v,
[t ]
x2 ~
T  T0  (T1  T0 )(1   T )
H
[t ]
Insertion and dropping tildes…..
1 v
1
(  v  v)  Ra (1  x2  T ) g  p    (v  (v)T )
Pr t
g
T
(
 v  T )   2T
t
02 gc p (T1  T0 ) H 3
( 106  108 )
Raleigh Number: Ra 
k
Prandtl Number: Pr 
 / 0
( 0.25 1018 )
k /( 0c p )
Relevant limit in Geophysics:
Pr  
Governing Equations, cont.
Stream function :
v1 

 , 2
x2
v2   ,1
and
In this way we satisfy the incompressibility constraint div v=0 identically.
Insertion into the velocity equations and the heat equation, assuming infinite
Prandtl number gives:
4  RaT,1  0
and
T,t  (,1T, 2  , 2T,1 )  ,1  2T
Dropping nonlinear terms and insert the “Ansatz”
(,1 ,T )  sin nx2eimx1 t
(m2  n2 2 ) 2 ,1  m2 RaT
, gives:
and
(  m2  n 2 2 )T  ,1
m2 Ra  (m2  n 2 2 )3

(m2  n 2 2 ) 2
Thus
Marginal instability:
 0
(m 2  n 2 2 ) 3
Ra 
m2
H
H
m
,n
L2 / 2
L1 / 2
For m=1 we obtain:
Ra min
27 4

, L1min  2 2 H
4
L1
((
Ra 
H
4H 2
)   2 )3
L1
4H 2
(
)
L1
Finite Element Approximations…
Ra=104, mesh: 128X128
The Nusselt Number
Definition:
Nu 
 (Tv
V
2
 kT, 2 )dV
  kT, 2 dV
total heatflux

flux by conduction only
V
It can be shown that for zero normal velocity b.c.’s and fixed top and bottom
Temperature
Nu  1 
1
RaV
1
T
T
(

v

(

v
)
)
:
(

v

(

v
)
)dV
V 
2 

Dimensionless MechanicalPower
The Nusselt Number
Hint for derivation of Nu-Power relationship
Ra (1  x2  T )e2  p    (v  (v)T )
form
 (.)  v dV
V
And apply Gauss’s Theorem. For the given b.c.’s it follows that
1
Ra  (1  x2  T )v2 dV   (v  (v)T ) : (v  (v)T )dV
2V
V
Finite Element Approximations…
Ra=104, mesh: 128X128
Nusselt Number
Finite Element Approximations…
Ra=105, mesh: 128X128
Finite Element Approximations…
Ra=105, mesh: 128X128
Nusselt Number
Finite Element Approximations…
Ra=106, mesh: 128X128
Finite Element Approximations…
Ra=106, mesh: 128X128
Nusselt Number
Finite Element Approximations…
Ra=107, mesh: 128X128
Finite Element Approximations…
Ra=107, mesh: 128X128
Nusselt Number
Galerkin Method
We consider the ansatz:
  A(m, t ) sin x2 cos mx1
T  B(m, t ) sin x2 cos mx1  C(m, t ) sin 2x2
Insert into:
1 2
  (  RaT )dx dx
,1
1 2
1
2
0
and
0 0
2
(
T

(

T


T
)




  ,t ,1 ,2 ,2 ,1 ,1 T )T dx1dx2  0
0 0
We obtain:
Ram 2
a  a  2
ab
2 3
(m   )
and
2
2
Ram
3
b   2b  2
a
(m   2 )3
Rayleigh-Taylor Instabilities
x2
x1
Consider two fluids of
different densities,
the heaviest above the
lightest. An horizontal
interface separates the two
fluids.
This situation is unstable
because of gravity.
Effectively, if the interfaces
modified then
a pressure want of balance
grows. Equilibrium can be
found again tanks to surface
tension that's why there is a
competition between surface
tension and gravity. Surface
tension is stabilizing instead
gravity is destabilizing for this
configuration.
x1
RT fingers evident in the Crab Nebula
Benchmark Problem
Rayleigh-Taylor
instability
benchmark
Method
mesh
γ0
(vrms)max
reached at
t=
van Keken
coarse
0.01130
0.003045
212.14
van Keken
fine
0.01164
0.003036
209.12
Particle-incell
1024 el
0.01102
0.003098
222
Particle-incell
4096 el
0.01244
0.003090
215
Level set
5175 el
0.01135
0.003116
215.06
Linear Instability Analysis
Ground state:
vi  0,
x2
n
h2
p  const  1 gx2 , x2  0 and
p  const   2 gx2 , x2  0
h1
Equilibrium to be satisfied in ground state at
time t=t0 and at t0+dt:
Stokes equation:
h1
 p    (v  (v)T )  g  0


S ( p,v)
Continued Equilibrium:
d 1/ 2
1/ 2
1/ 2
S  S , t  S , k vk  0

dt
0
Linear Instability Analysis
with ui  vi ,t and P  p,t
x2
n
h2
we obtain
h1
 P,i1/ 2  1/ 2 (ui1,/j2  u1j/,i2 ), j  0
h1
Or, considering the
incompressibility constraint: u1j/, 2j  0
 P,i1/ 2  1/ 2ui1,/jj2  0
Linear Instability Analysis: Boundary and interface
conditions
x2
n
h2
We assume that the velocities are zero on
top and bottom; on the sides we assume
that v1=0 and s120, i.e. symmetry
boundary.
On the interface the velocities as well as
the natural boundary termsh have to be
continuous.
Natural b.c.’s: replace () , j in
1
h1
(P  ij  (u
1/ 2
1/ 2
1/ 2
i, j
 u )), j  0
1/ 2
j ,i
By n j . The vector
ti  (P1/ 2 ij 1/ 2 (ui1,/j2  u1j/,i2 ))n j
as well as its time derivative have to be continuous on the interface
Linear Instability Analysis: Boundary and interface
conditions
x2
n
h2
We assume that the velocities are zero on
top and bottom; on the sides we assume
that v1=0 and s120, i.e. symmetry
boundary.
On the interface the velocities as well as
the natural boundary termsh have to be
continuous.
Natural b.c.’s
d
((s ij n1j  s ij2 n 2j )dA0 )  0
dt
Result:
1
h1
x1 :  1 (u11, 2  u12,1 )   2 (u12, 2  u22,1 )  0
Or:
x2 :  1u12, 2   2u22, 2  ( P1  P 2 )  (  1   2 ) gv2  0
Exercises
1. The normal component of the surface velocity of a 2D half
plane x2  0 reads v2  V cos kx1 . The half plane is occupied by
a Stokes fluid with the viscosity . The normal stress is obtained as
s 22  Q cos kx1 . Show that Q  2kV .
2. Consider a Rayleigh-Taylor instability problem involving 2
infinite half-planes; i.e. h1k , h2k  1 . The normal velocity of the
Interface reads v2  Vet cos kx1 . Show that  
 ( 1   2 ) g
.
2(1   2 )k
Hint: Use the relationship Q,gt 0  2kU , U  V,t  U for
a gravity free Half-plane and note that Q,gt 0  Q,t  gV
.
Exercises (folding)
w
Exercise (folding) cont.
Exercise (folding) cont.
Excercise 4:
Solve convection equations
   RaT,1  0
4
and
T,t  (,1T, 2  , 2T,1 )  ,1  2T
using the perturbation expansion
  1   2  .... T  T1   T2  ....    t
2
2
2
Ra  Ra 0   2 Ra 2
up to terms of order  3 .
Hint:
Insert into pde’s, collect coefficients of , 2 and 3
Individual coefficients must be equal to zero. Get
1, 2 etc
3. 3 and T3 are a little bit harder to get since the pde’s
contain inhomogeneous terms which are proportional
to (1, T1)
1.
2.
Perturbation solution for the weakly nonlinear
problem
4.
We require a solubility condition (Fredholm’s
alternative) In the present case this just means
that the coefficient of the resonant inhomogeneous
term must vanish. The coefficient has the form of an
ode which can be written as:
a
0
m2 Ra 2
1 3
a,  2
a a
2 2
(m   )
2
Racrit
Ra