Slides

13 downloads 1497 Views 2MB Size Report
... Fluid Dynamics. Consider the initial-boundary value problem for the one ..... Computational Fluid Dynamics where and at each node e1 e2 e5 e6 e3 e7 e4 e8 .
Computational Fluid Dynamics

http://www.nd.edu/~gtryggva/CFD-Course/!

Computational Fluid Dynamics

Finite Volume Methods are used in the vast majority of CFD codes. However, there are other ways to compute the motion of fluids. Some of the more common alternatives include:!

Other! CFD Methods!

Spectral Methods! Finite Element Methods! Lattice Boltzmann Methods! Smoothed Particle Hydrodynamics! Vortex Methods and other particle methods!

Grétar Tryggvason! Spring 2011!

Here we will briefly examine those!

Computational Fluid Dynamics

Computational Fluid Dynamics

Consider the initial-boundary value problem for the one dimensional heat equation!

∂f ∂ 2 f = 0 < x < π, t ≥ 0 ∂t ∂ x 2 f (0,t) = f (π ,t) = 0 f (x,0) = f (x)

Spectral! Methods!



By standard technique (Fourier sine transform) the solution is!

Computational Fluid Dynamics

Assuming that the time integration is exact, the error by using only N terms is! f − f N ≈ e−N

2

t

for all time!

f (x,t) = ∑ an (t)sin nx n=1

an (t) = f n e−n t 2 π f n = ∫ f (x)sin nx dx π 0 2

Computational Fluid Dynamics

A spectral approximation to the heat equation is now simply! N

f (x,t) = ∑ an (t)sin nx n=1

This is faster than any power of N. Recall that the error for finite difference/finite volume methods is! O(h α ) ≈

1 Nα

Spectral methods are said to have “spectral” or “infinite order” accuracy! This property is not retained if the solution has discontinuities!

Sum to N only!!

Substituting into the heat equation we get:!

dan = −n 2 an ; an (0) = f n dt

(n = 1,2,…N)

Computational Fluid Dynamics

For nonlinear equations the use of spectral methods becomes much more involved! ∂f ∂f ∂ f + f =ν 2 ∂t ∂x ∂x 2

Expand the solution in a Fourier series! f ( x,t) =

∑ a (t)e

2π kx L

−i

k

|k|< K

We get the following equation for the coefficients:!

Computational Fluid Dynamics

Computing the nonlinear term can become expensive, particularly for three-dimensional systems. By using the Fast Fourier Transform, the computation can be made faster by taking the inverse transform and do the multiplication in real space. Then Fourier transform the product back into spectral space. The multiplication can, however, lead to aliasing problems where the resolved wave numbers are contaminated by unresolved ones.!

⎛ 2πl ⎞ 2 dal iπ + l ∑ al−m am = −ν ⎜ ⎟a ⎝ L ⎠ l dt L3 / 2 |m|< K |l -m|< K

Computational Fluid Dynamics

The fast Fourier transform! N

f j ≡ ∑ fˆl e

i

2π l j N

l=1

inverse! N

fˆl ≡ ∑ f j e

−i

2π lj N

l=1

Can be evaluated in 2N log2N operations!

Computational Fluid Dynamics

Although spectral methods can have very high accuracy, they can be difficult to apply to nonlinear problems, problems with complex boundary conditions, and complicated geometries.! In addition to Fourier series, several other basis function can be used for other geometries and to control the resolution near boundaries! Spectral elements methods, where a complex domain is broken into large blocks and each block is treated spectrally can be used treat somewhat complex boundaries !

The Cooley-Tukey algorithm!

Computational Fluid Dynamics

Finite Element Methods!

Computational Fluid Dynamics

Finite element methods are similar to spectral methods in that we expand the solution in terms of a known basis function. Unlike spectral methods, where the basis functions are defined globally over the whole computational domain, in the finite element method the basis functions are defined locally on each element.!

f ( x) = ∑ f I N I ( x) I

Element 1!

∂f ∂N I = ∑ fI I ∂x ∂x

Element 2!

Computational Fluid Dynamics

Computational Fluid Dynamics

Linear basis functions! (1) N i−1

For element 1:!

N

N i(2)

(1) i

i −1 element 1! i

For element 1:!

N i(1) = For element 2:!

N i(2) =

N

element 2!

x − x i−1 ; Δx i

(1) N i−1 =

xi − x Δx i

x i+1 − x ; Δx i+1

(2) N i+1 =

x − xi Δx i+1

N i(1) =

(2) i+1

i +1

For element 2:!

N i(2) =



Ω

W

∂ ⎛ ∂f ⎞ ⎜ k ⎟dx = ∂ x ⎝ ∂x ⎠



∂ ⎛ ∂f ⎞ R = ⎜k ⎟ − q ∂x ⎝ ∂ x ⎠

∂ ⎛ ∂f ⎞ ⎜ k ⎟dx = ∂ x ⎝ ∂x ⎠



W(x): Weighting function!

W (x)R(x)dx = 0 Ω

Ω



Ω

Ω

I

q W dx

For node i:! i +1

∑ ∫

j =i −1

fj

i −1

k

q Wdx



Ω

∂ ⎛ ∂f ⎞ W ⎜ k ⎟dx = ∂x ⎝ ∂x ⎠



Ω

q Wdx

i +1 ∂ N j ∂ Ni dx − ∫ qN i (x)dx = 0 ∂x ∂x i −1



Ω

q Wdx

Computational Fluid Dynamics

Take weighting function equal to the interpolation function!

i +1

Ω

Cancel for adjacent elements!

f ( x) = ∑ f I N I ( x)

Ω

(2) ∂N i+1 1 = ∂x Δx i+1

⎡ ⎛ ∂f ⎞⎤ ⎡ ⎛ ∂f ⎞⎤ ⎛ ∂f ⎞ ∂W dx + ⎢W ⎜ k ⎟⎥ − ⎢W ⎜ k ⎟⎥ = ⎜k ⎟ ⎝ ∂x ⎠ ∂x ⎣ ⎝ ∂x ⎠⎦2 ⎣ ⎝ ∂x ⎠⎦1

q Wdx

Galerkin Method!





⎛ ∂f ⎞ ∂W dx + ⎜k ⎟ ⎝ ∂x ⎠ ∂x

Computational Fluid Dynamics

⎛ ∂ f ⎞ ∂W ∫Ω ⎝⎜ k ∂ x ⎠⎟ ∂ x dx =

x − xi Δx i+1

Evaluating the second term:!

The weighted residual is:!

W

(2) N i+1 =

Integrate by parts!

The residual is:!

Ω

x i+1 − x ; Δx i+1

Computational Fluid Dynamics

∂ ⎛ ∂f ⎞ ⎜k ⎟ = q ∂x ⎝ ∂ x ⎠



xi − x Δx i

(1) ∂N i−1 −1 = ∂x Δx i

∂N i(2) −1 = ; ∂x Δx i+1

Example: Elliptic Problem!

Or:!

(1) N i−1 =

∂N i(1) 1 = ; ∂x Δx i

Computational Fluid Dynamics



x − x i−1 ; Δx i

i +1

i +1

∑ f ∫k

j =i −1

j

i −1

i +1 ∂ N j ∂ Ni dx − ∫ qN i (x)dx = 0 ∂x ∂x i −1

integrating!

∂N i(1) 1 = ; ∂x Δx i (1) ∂N i−1 −1 = ∂x Δx i ∂N i(2) −1 = ∂x Δx i+1

∂N 1 = ∂x Δx i+1 (2) i+1

ki +1/2

fi +1 − fi f −f q + 4qi + qi +1 − ki −1/2 i 2i −1 − i −1 =0 Δx 2 Δx 6

where!

ki +1/2

1 = Δx

i +1

∫ k dx i

Notice that a finite difference method would give the same result except with qi on the RHS!

∂N i(1) 1 = ; ∂x Δx i (1) ∂N i−1 −1 = ∂x Δx i ∂N i(2) −1 = ∂x Δx i+1 (2) ∂N i+1 1 = ∂x Δx i+1

Computational Fluid Dynamics

Finite Element approximation of an unsteady problem:!

∂f ∂f +U =0 ∂t ∂x

Computational Fluid Dynamics

∂f

∂f

f ( x) = ∑ f I N I ( x) I

Ω

Introduce the finite element representation!

∂f I

∑ ∂t ∫ N N

Weak formulation!

I

∂f

∫ ∂t W dx + ∫ U ∂x W dx = 0

Ω

∂f

∫ ∂t W dx + ∫ U ∂x W dx = 0

Ω

Ω

Introduce the finite element representation!

f ( x) = ∑ f I N I ( x) I

I

or!

Ω

∂f I

∑ ∂t

J

dx + U ∫ N I

M IJ + U ∑ K IJ = 0

I

M IJ =

∫ N I N J dx;

Ω

“Mass matrix”!

Computational Fluid Dynamics

Ω

dN J dx = 0 dx

I

K IJ =

∫N

Ω

I

dN J dx; dx

“Stiffness matrix”!

Computational Fluid Dynamics

In two-dimensions triangular (or tedrahedron) elements are frequently used! The mass matrix makes it complex to use simple explicit time integrators. This can be overcome by “lumping” the matrix!

Computational Fluid Dynamics

While finite elements dominate solid mechanics, their role in computational fluid dynamics has been secondary to finite volume methods. Indeed, most methodological progress has been first developed in the context of finite volumes and then implemented in the finite element framework. Some areas, such as simulations of viscoelastic flows have, however, relied heavily on finite elements for historical reasons.!

Computational Fluid Dynamics

Lattice Boltzmann Methods!

Computational Fluid Dynamics

Lattice Boltzmann methods !

Computational Fluid Dynamics

Lattice Boltzmann methods !

Lattice Boltzmann methods use “particles” on hexagonal grids where the particles move according to discrete rules. It can be shown that the macroscopic motion of the particles resembles the Navier-Stokes equations. It is now believed that these methods can yield results comparable in accuracy to other numerical methods!

Update the velocity distributions at each node by:! ∂ fγ 1 (eq) + eγ ⋅ ∇fv = fγ − fγ , γ = 0,1,, M ∂t ετ where! 2 9 3 2⎤ ⎡ fγ(eq) = wγ ρ ⎢1 + 3 eγ ⋅ u + eγ ⋅ u − ( u ⋅ u ) ⎥ 2 2 ⎣ ⎦

(

(

)

) (

)

and at each node! M

ρ = ∑ fγ γ =1

e2

e6

M

ρu = ∑ fγ eγ

e1

γ =1

e3 e7

Computational Fluid Dynamics

Lattice Boltzmann methods ! Proof!

(

Set! α = Δt eγ

))

(

)

α = Δt eγ Δx; β = Δx ε

γ

γ

γ

Δt eγ

ετ

( f (x,t) − f (x,t)) (eq ) γ

γ

1 (eq ) ( fγ (x,t ) − fγ (x,t )) τ

Summary:!

Δx ⎛ π π⎞ ⎛π⎞ ⎛π⎞ cos ⎜ ⎟ (γ − 1) + ,sin ⎜ ⎟ (γ − 1) + ⎟ , γ = 2, 4, 6, 8 ⎝ 2⎠ ⎝ 2⎠ Δt ⎜⎝ 4 4⎠

e6

wγ = 1 36 v = 2, 4, 6, 8

( f (x + Δte ,t ) − f (x,t )) −

)

by!x + Δteγ

Update the velocity distributions at each node by:! 1 fγ ( x + Δteγ ,t + Δt ) − f γ ( x,t ) = ( f γ(eq ) ( x,t ) − f γ ( x,t )) τ where! 2 ⎡ 9 3 2⎤ fγ(eq ) = wγ ρ⎢1+ 3(eγ ⋅ u ) + (eγ ⋅ u ) − (u ⋅ u ) ⎥ ⎣ ⎦ 2 2

Δx ⎛ ⎞ ⎛π⎞ ⎛π⎞ cos ⎜ ⎟ (γ − 1) ,sin ⎜ ⎟ (γ − 1)⎟ , γ = 1, 3, 5, 7 ⎝ 2⎠ ⎝ 2⎠ ⎠ Δt ⎝⎜

w0 = 0 wγ = 4 9 γ = 1, 3, 5, 7

Δx

Replace! x

Computational Fluid Dynamics

Lattice Boltzmann methods !

The displacement vectors are!

The weights must be taken to be!

Δt eγ

fγ ( x + Δteγ ,t + Δt ) − f γ ( x,t ) =

1 (eq ) ( fγ (x,t ) − fγ (x,t )) τ

e0 = 0

eγ = 2

(

Δx = Δt eγ

And crossing out the common term yields!

Computational Fluid Dynamics

Lattice Boltzmann methods !

eγ =

β fγ ( x − Δteγ ,t ) − f γ(eq ) ( x − Δteγ ,t ) τ

taking! Δt = Δx = Δy = ε

Δx = Δt eγ

taking! Δt = Δx = Δy = ε And replacing! x by! x + Δteγ yields!

fγ ( x + Δteγ ,t + Δt ) − f γ ( x,t ) =

Δx; β = Δx ε

fγ ( x + Δteγ ,t + Δt ) − f γ ( x + Δteγ ,t ) = −

β − f γ ( x − Δteγ ,t ) − f γ(eq ) ( x − Δteγ ,t ) τ

(

)

fγ ( x,t + Δt ) − f γ ( x,t ) = −α f γ ( x,t ) − f γ ( x − Δteγ ,t ) −

Approximate the time derivative by a first order forward in time, the advection term by upwind and the collision by downwind, we get! fγ ( x,t + Δt ) − f γ ( x,t ) = −α f γ ( x,t ) − f γ x − Δteγ ,t where!

e8

e4

Computational Fluid Dynamics

Lattice Boltzmann methods !

∂f γ 1 + eγ ⋅ ∇f v = ( f γ(eq ) − f γ ), γ = 0,1,, M ∂t ετ

(

e5

e2

e5 e1

e3

and at each node! M

ρ = ∑ fγ γ =1

e7

e4

e8

e6

e2

e5

M

ρu = ∑ f γ eγ γ =1

e1 e3 e7

e4

e8

Computational Fluid Dynamics

Lattice Boltzmann methods !

Computational Fluid Dynamics

Lattice Boltzmann methods !

Advantages! Very easy to program–parallelization easy! Fast–in part because incompressibility is not enforced! Disadvantages! Difficult to implement new physics—must find a discrete equation that corresponds to the new physics! The flow is not completely incompressible!

By now it is understood that LBM are give results comparable to second order finite difference/volume methods! LBM has been extended to more complex flow, including multiphase flows with sharp interfaces! The figure shows a comparison of the unsteady rise of a buoyant bubble computed by LBM (left) and front tracking finite volume method (right).!

Computational Fluid Dynamics

From: K. Sankaranarayanan, I.G. Kevrekidis, S. Sundaresan , J. Lu, G. Tryggvason. A comparative study of lattice Boltzmann and front-tracking finitedifference methods for bubble simulations. Int. J. Multiphase Flow 29 (2003) 109–116.!

Computational Fluid Dynamics

Smooth Particle Hydrodynamics! Smooth Particle Hydrodynamics: the fluid mass is lumped into smoothed blobs that are moved using Newtonʼs second law directly, without an underlying mesh!

Smooth Particle Hydrodynamics! x2

x1

Computational Fluid Dynamics

Smooth Particle Hydrodynamics!

Computational Fluid Dynamics

Smooth Particle Hydrodynamics!

In SPH the fluid is modeled as a collection of smooth “blobs” or particles !

∫ A(r ')δ (| r − r ' |) dr '

We could select the kernal to be a Gaussian, although usually we desire a compact support!

And then approximate the delta function by a smoother “kernal”! A(r) =

∫ A(r ')W (| r − r ' |, h) dr '

⎛r2 ⎞ W (r, h) = 3/ 2 exp ⎜ ⎟ h2 ⎠⎟ ⎜⎝ π h 1

Where:!

∫ W (r, h) dr = 1

Approximate the integral by a sum over the smoothed particles:! A(r) =

We start by the relationship! A(r) =

()

& W (r, h) ⎯h→0 ⎯⎯ →δ r

x3

=

∫ A(r ')W (| r − r ' |, h) dr ' ρ(r ')

∫ A(r ') ρ(r ') W (| r − r ' |, h) dr '

mk ≈ ∫ ρ(r ')dr '

A = ∑ mk k W (| r − rk |, h) ρk k

For density we have! ρi = ρ(ri ) = ∑ mk k

ρk W (| ri − rk |, h) = ∑ mkW (ri − rk , h), ρk k

Computational Fluid Dynamics

Smooth Particle Hydrodynamics! The introduction of the smoothed kernal allows us to compute the gradient! A ∇A(r) = ∑ mk k ∇W (| r − rk |, h) ρk k

Computational Fluid Dynamics

Start with the momentum equation! du −1 = ∇P dt ρ

∇A(r) = ∑ mk k

For the pressure we write:!

The introduction of the smoothed kernal allows us to compute the gradient!

Ak ∇W (| r − rk |, h) ρk

⎛ P⎞ P 1 ∇P = ∇ ⎜ ⎟ + 2 ∇ρ ρ ⎝ ρ⎠ ρ

And using that! ∇ρ = ∑ mk ∇W (| r − rk |, h) k

P ⎛ P⎞ ∇ ⎜ ⎟ = ∑ mk k2 ∇W (| r − rk |, h) ⎝ ρ⎠ ρk k

The momentum equation is! ⎛P du i P⎞ = − ∑ mk ⎜ k2 + i2 ⎟ ∇W (| ri − rk |, h) dt ⎝ ρ k ρi ⎠ k

Computational Fluid Dynamics

Computational Fluid Dynamics

Similarly, for density we start with! d ρi = − ρ∇ ⋅ u dt

(

)

The energy equation is ! ∇A(r) = ∑ mk

i

k

Ak ∇W (| r − rk |, h) ρk

And using that! k

The mass conservation equation is!

(

)

k

mk

ρk

(u

k

⎛P ⎞ ⎛ P⎞ P ∇ ⋅ u = ∇ ⋅ ⎜ u⎟ − u ⋅ ∇ ⎜ ⎟ ρ ⎝ρ ⎠ ⎝ ρ⎠

dEi dt

Sometimes an alternative form is used! = ρi ∑

k

Ak ∇W (| r − rk |, h) ρk

= ∑ mk k

Pk

ρk2

(u

k

)

− u i ⋅ ∇W (| ri − rk |, h)

)

Summary of governing equations! d ρi m = ρi ∑ k u k − u i ⋅ ∇W (| ri − rk |, h) dt k ρk

)

Mass! W (s) =

Momentum!

dEi P = ∑ mk k2 u k − u i ⋅ ∇W (| ri − rk |, h) dt ρk k

Energy!

dri = ui dt

Advect particles!

)

Plus an equation of state!

P = P( ρ,e)

Computational Fluid Dynamics

Smooth Particle Hydrodynamics! Several different kernals have been used. The following is a popular one:!

⎛P du i P⎞ = − ∑ mk ⎜ k2 + i2 ⎟ ∇W (| ri − rk |, h) dt ⎝ ρ k ρi ⎠ k

(

k

− u i ⋅ ∇W (| ri − rk |, h)

Computational Fluid Dynamics

Smooth Particle Hydrodynamics!

(

∇ρ = ∑ mk ∇W (| r − rk |, h)

gives!

d ρi = ∑ mk u k − u i ⋅ ∇W (| ri − rk |, h) dt k

dt

∇A(r) = ∑ mk

Using that! ∇ρ = ∑ mk ∇W (| r − rk |, h)

( ρ∇ ⋅ u) = ∇ ⋅ ( ρu) − u ⋅ ∇ρ

d ρi

dE P = − ∇⋅u dt ρ

10 7π h2

⎧ 3 3 ⎪ 1 − s 2 + s3 0 ≤ s ≤ 1 2 4 ⎪⎪ ⎨ 1 2 1< s ≤ 2 ⎪ 4 (2 − s) ⎪ 0 2