Transcript Slide 1

23.3.2006
ISCM-20
Pinhas Z. Bar-Yoseph
Computational Mechanics Lab.
Mechanical Engineering, Technion
Copyright by PZ Bar-Yoseph ©
dy
 f  y, t 
dt
t u   j f  g
j
 t u   c F  u    v F  u, u   0
The time interval I   0, T 
t
is partitioned into subintervals I   t , t
n
n
n 1

where t n and t n 1 belong to an ordered partition
of time levels 0  t 0  t1 
 t N  T.
t n 1
tn
In
t n 1
Each time subinterval is represented by one spectral element.
Within each spectral element, the dependent variables are
expressed in terms of p-th order Lagrangian interpolants
through the Legendre-Gauss-Lobatto points.
dy
 f  y, t  , t  0
dt
0
y  0  y
Discontinuous Galerkin Method (DGM)
 T dy

T
w
,

f
y
,
t

w
, y



n
dt

I
t
n
0
where
y
  y   y acts as a stabilizer operator
Jump
operator

v  lim v  x  z 
z 0
dT
 T  0, t  t 0
dt
T t 0   T 0
Method of Weighed Residuals (MWR)
dT
R
  T , t  t 0
dt
rT  T  t
n


T
n
t n1
  T R dt   T r
tn
T t t n
0
t n1

tn
 dT

T 
  T  dt   T T  t n   T n
 dt


rT

0
t t n
R
T  S  Space of piecewise polynomials of degree p  0
with no continuity requirements
across interelement boundaries:

S  v  L2  I  : v I n  P p  I n  , I n  I

Galerkin + Linear element
T  N nT n  N n 1T n 1
where T  T  t
n
K C M
e
e
n

e
1  1 1
 2 1
e
C  
,
M  


2  1 1
6 1 2
n

T 
e
U   n1  ,   t
T 
e
Continuous Galerkin Method (CGM)
0 
F  
0 
T  T  tn  ,

n
T
n 1
e
 3 n

T 
A Tn
2  3
Amplification
A
factor
A  1  A-Stable
Discontinuous Galerkin Method (DGM)


T  tn  
e
F 
,
 0 
T
n 1
1 0 
K K 

0 0
e
e
2 3    n
 n
 2
T AT
  4  6
A
A 1
lim A =0  Asymptotic annihilation 
 
L -Stable
Bar-Yoseph, Appl. Num. Math. 33, 435-445, 2000
Aharoni & Bar-Yoseph, Comp. Mech. 9, 359-374, 1992
DSM for Dynamic systems

H
0  q  p
t
tf


 
H
 Q   q dt
 p   p 
q


 
 p 0  p  t 0    q  t 0 
 q 0  q  t    p  t
0
0
0
H  Hamiltonian
q  Generalized displacement
p  Generalized momentum
Discontinuous element
t
n 1
 t 
t
p
n 1
, q
th
n time step
n
t
n
t
n
t n 1
n1
n
p , q
p t
n
n
 , q t 
n
Nonlinear Spatio-Temporal Dynamics of a
Flexible Rod
Plat & Bar-Yoseph, 27th Israel Conf.
Mech. Eng. 683-685, 1998
Bar-Yoseph, Appl. Num. Math. 33, 435-445, 2000
Nave, Bar-Yoseph & Halevi, Dynamics. & Control. 9, 279-296, 1999
The unicycle system, presents an
example of inherently unstable system
which can be autonomously controlled
and stabilized by a skilled rider
uW -required to maintain the
unicylce’s upright position
uT
-required to maintain
lateral stability
M f -the friction torque
is assumed to be
dependent on the
yew rate  only
The adaptive technique performed
very well for all stiff systems that we
have experienced with (convection,
radiation and chemical reactions),
and is competitive with the best
Gear-type routines
 t u   j f  g,
j
j  1(1)ndim
Linear Eq's. (scalar wave eq.)
Nonlinear Eq's.
DGM's milestones papers:
Reed & Hill, " Triangular mesh methods for the neutron transprt equation",
Proc. Conf. Math. Models &
Comutational Techniques for Analysis of Nuclear Systems, Michigan, 1973 .
Lesaint & Raviart, " On a finite element method for solving the neutron transprt equation",
Mathematical Aspects of Finite Elements in PDE's, 1974
u
u
c
 0,
t
x
IC ' s
 x, t ]x0 , x1[
u  x, 0   f  x 
BC ' s
 ]0, T [
f  x
u  x0 , t   u  x1 , t 
where
c0
x
u  x
Exact solution
Continuous Galerkin
Classical artificial diffusion
Discontinuous Galerkin
x

t
3
G3
G1
 u  0 in G   I
u  u
ν
G

2
on  
w udG 

G3


w u d  0
3
x
p0



w u d  


3
u3 


3

u


w u


3

wu

 n  ν d   0
3
 n  ν d    n  ν d 

3
1
2
1
u 
u 
u2
 1  2
 1  2
3
 n  ν d   
First Order
Upwind FD scheme
x
t
f
G
p 1
e
t
x
t
x
t
n
n 1
j 1
j
j 1
x
DGM
u 
 u
e
e
e
ˆ
ˆ
w

c
d


w
u
d


w
cu
d

t
x 0
e  t x 
e
e
G
t
x
n
n 1
n
D

c
D

L

c
R
U

L
U

c
R
U
 t
x
i
i
j
0
j
0
j 1
n indicates the time slab
j indicates element number within the time slab
where
N
e
T N
Dt   W
d  , Dx   W
d e
t
x
Ge
Ge
T
Li 
ˆ T Nd  e , L 
W
t
0

te
ˆ T Nd  e
W
t

te
ˆ T Nd  e , R  c W
ˆ T Nd  e
Ri  c  W
x
0
x

ex
ex
Space-Time Discontinuous Approximations
u  S  Space of piecewise polynomials of degree p  0
with no continuity requirements
across interelement boundaries:


S  v  L2  G  : v G  P p  nG e  ,  nG e  G
Conventional el.
Bar-Yoseph & Elata, IJNME, 29, 1229-1245, 1990
"Gauss-Lagrange" el.
Basis functions are po int wise orthogonal
at the integration points.
ˆ = N and bilinear element:
For W = W
 2

x  1
Dt  xDt 
12  2

 1
1 2 1 
2 1 2 
,
1 2 1 

2 1 2 
 2

t  2
D x  tD x 
12  1

 1
2 1 1 
2 1 1 
1 2 2 

1 2 2 
4

x  2
Li  xLi 
12  0

0
2 0 0
4 0 0 
,
0 0 0

0 0 0
0

x 0
L 0  xL 0 
12 0

0
0 4 2
0 2 4 
0 0 0

0 0 0
4

t  0
R i  tR i 
12  2

0
0 2 0
0 0 0 
,
0 4 0

0 0 0
0

t 0
R 0   tR 0 
12 0

0
4 0 2
0 0 0 
2 0 4

0 0 0
Von Neumann Analysis
U
n
j 1
U 
n
j
j 1
U
n
0
where
  exp  i 2 / k  ,
k  L / x,
wave
number

i  1
wave
length

Dt   D x  Li   R i   R 0 1 U 0n  L 0 U 0n 1
A
Amplification
matrix

 D t   D x  L i   R i   R 0
1

1
L0
Bar-Yoseph & Elata & Israeli,
IJNME, 36, 679-694, 1993;
Golzman & Bar-Yoseph (Project)
Bar-Yoseph & Elata & Israeli,
IJNME, 36, 679-694, 1993;
Golzman & Bar-Yoseph (Project)
Bar-Yoseph & Elata, IJNME, 29, 1229-1245, 1990
Adaptive
Refinements
Note : This problem can be solved by use of tilted elements,
where two elements are sufficient to reconstruct the exact solution.
Fischer & Bar-Yoseph, IJNME, 48, 1571-1582, 2000
Advanced CAD Visualization Methods
Adaptive Level of Details
Technique for Meshing
Morphing between Meshes at Different Times
CGM
Conforming elements.
Elemental Contributions
are assembled to generate
the GLOBAL set of equations.
DGM
Elements are discontinuous.
The size of the LOCAL equations is equal to
e
the ndof
inside the corresponding element.
The element matrices can be inverted by
using symbolic manipolator once and for all!
CGM
"Conforming" elements
DGM
"Non-conforming"
elements
The degree of the
approximating polynomial
can be easily changed
from one element to the other.
Space-Time Discontinuous Approximations
u
t
x
DGM remains compact
with high-order polynomial basis
(essential for unstructured mesh computation)
Discontinuous
SPECTRAL ELEMENTS
•Gauss-Lobatto nodes are clustered near element
boundaries and are chosen because of their
interpolation and quadrature properties.
• Mass lumping by nodal quadrature.
• Exponential rate of convergence.
•The increase in the ndof due to the discontinuity
at the interelement boundaries is balanced in high
order elements.
DGM can easily handle adaptivity strategies
since refinement or unrefinement of the mesh
can be achieved without taking into account
the continuity restrictions
typical of conforming FEM (needs transition elements).
x
h - p type of convergence
can be easily implemented,
since no continuity requirement
is a priori imposed on the
test and trail space of functions.
t
x
Parallel adaptive DGFE computation.
t
DGM are highly parallelizeable.
14
15
12
9
6
1
V
13
10
7
3
16
8
11
5
4
2
x
The compact (nearest neighbor) scheme
minimizes interelement communication.
in  , j  11M
t u   j f j  0
Flux Splitting

 i
f f f
i
i
w ,  u  
t
j f
j

e
N
N
  w , [[ f i]]
e
i 1
 w , [[u]]
e
t
2 i1
0
Bar-Yoseph,Comput. Mech., 5, 145-160, 1989
  w , [[ f i]]
i 1
w V h ;
e
2 i
, j  11M
Miles Rubin (2005)
Nonlinear Wave Eq.
2
u

 2u
2
,
c

2
2
X
t
where



2 
c


c 2  0 1 
2
1   
u  
 
 1 
  X  
E
c02  0
0
Flux splitting for non homogeneous f i u


f  i f i i 0
f Au u


f  i f i i 0
f Au u
Traper & Bar-Yoseph (Project)
 2 u f
  0 , x  (0, L)
2
t
x
where:


1 
 
u2   1 


 1  u2  
f   c 2

0
1 




 u1


 u1 
u 
u 2 
The effective wave speed:
1
c c
1 
2
In a matrix form:

1
1 0   u1   0  c02
1 
0 1 t u   

  2 
 1
2
0



1 

 1  u 2 
2





1 
   u1  0
2
 1  u  
 
2

 x 
u2  0


0
• The Jacobian matrix:

1
2
0

c
0
Au   
1 

 1



1 

2
 1  u  
2



0
• The eigenvalues:
1   1  2u2  u2 2   
1  c0
1   1  u2 
1   1  2u2  u2 2   
2  c0
1   1  u2 
• The corresponding eigenvectors:


v1   c0


1   1  2u2  u2 2    

1   1  u2 
1


v 2  c0




1   1  2u2  u2 2    

1   1  u2 
1




2
1


1    1  2u2  u2  
1 2   1  u2
 c0
u1  c0
u2 

1   1  u2 
2 1   1  u2  
2

f  f1  

2
1    1  2u2  u2     1  u2  
1
1

u2 
2
  2 u1  2 c0
1  u2  1  2u2  u2  








2
 1


1    1  2u2  u2  
1 2   1  u2
 c0
u1  c0
u2 

1   1  u2 
2 1   1  u2  
 2

f  f2  

2
1    1  2u2  u2     1  u2 
1
1



u

c
u
1
0
2
2


2
2


1


u
1

2
u

u


2
2
2







2
1

1    1  2u2  u2  
 c0

1   1  u2 
2

A  A1  
1



2


2
 1

1    1  2u2  u2  
 c0

1   1  u2 
 2

A  A2  
1



2



1    1  2u2  u2 2     1  u2 

1  u2  1  2u2  u2 2  
1
  1  u2
 c02
2 1   1  u2 
1
c0
2








1    1  2u2  u2 2     1  u2 

1  u2  1  2u2  u2 2  
1
  1  u2
 c02
2 1   1  u2 
1
 c0
2




Traper & Bar-Yoseph (Project)
 x, t 
Displacement
  0.5, L  1, c0  1, 0  0.2
Velocity
  0.5, L  1, c0  1,0  0.2
Strain
  0.5, L  1, c0  1,0  0.2
• Tcr -Time for breakdown [Lax (1964)]:
 

Tcr  2  max K,u2 0 , xx 0
1
K  c0

1  u2 2
1 
K ,u2 0   c0

1 
maxu 2 , x 0   0.2
2
L2
21   L2
 3.0396
Tcr 
2
0.2c0
1
Coarse (& Uniform) grid
bilinear
biquadratic
u1
x
Velocity at t = 3 sec
  0.5, L  1, c0  1,0  0.2
u2
bilinear
biquadratic
x
Strain at t = 3 sec
  0.5, L  1, c0  1,0  0.2
Explicit Vs. Implicit schemes
The discontinuous approximation
can capture shock waves
and other discontinuities
with accuracy.
Bar-Yoseph et al., JCP, 119, 62-74, 1995
Bar-Yoseph & Moses, IJNMHFF, 7, 215-235, 1997
Runge-Kutta Discontinuous Galerkin (RKDG ; LDGM)
 t u   c F  u    v F  u, u   0
 t u   F  u, u   0
where
F  u, u   c F  u   v F  u, u 
d
 u u d     u F  u, u  n d      u  F  u, u  d   0

dt e
e
e
 H  numerical
flux
Cockburn& Shu, JCP, 84, 90, 1989; Basi & Rebay, JCP, 131,267-279, 1997
Runge-Kutta Discontinuous Galerkin (RKDG)
 t u   c F  u    v F  u , u   0
(1)
D  u  0

c
v

u


F
u


F  u, D   0


 t
(2)
(3)
(2)
 D D d    D
e
e
u n d 
D
H
numerical
flux
   D  u d   0
e
(4)
(3)
d
u u d    u

dt e
e

c
F  u  n d      u  e F  u  d 
e
 H  numerical
flux
c



   u v F u, D n d       u  v F u , D d   0
e
e
 H  numerical
flux
v
Note : When evaluating boundary integral of (5) along e ,


the flux terms, u n, c F  u  & v F u, D , are not uniquely defined
due to the discontinuous approximation.
c
c
H
c
F u  ,
H 1 2  v F


D
H 1 2   u   u  n,

u,  D  v F



u,  D  n.
(6)
(5)
Semi-discrete method
dU
M
 K U  0
dt
This system of ODE's is integrated
with a Runge-Kutta method.