A Second Order Finite Volume Method for a Multi-material Heat - - PowerPoint PPT Presentation

A Second Order Finite Volume Method for a Multi-material Heat Equation

• n Cartesian Grids:

Application to Stefan problems.

ECCOMAS 2012 | Manuel LATIGE, Gerard GALLICE, Thierry COLIN September 10-14 2012

Scientific Background

General problem

The solid part can be composed by several materials with different properties One of the solids undergoes a phase change (fusion)

We need to ...

handle several interfaces. Impose the jump conditions across the interfaces

Goals ...

Develop a spatial second-order Finite Volume method with a compact stencil for the heat equation in 2D.

CEA, INRIA | September 10-14 2012 | PAGE 1/18

Heat equation on multi-material domain

Heat equation on multi-material domain

ρ1,2C1,2∂tT − ∇. (K1,2∇T) = 0 on Ω1,2/Γ. (1) [T]Γ = T|Ω1 − T|Ω2 = 0, (2)

• K ∂T

∂n

• Γ

= K1 ∂T ∂n

• Ω1

− K2 ∂T2 ∂n

• Ω2

= 0, (3) where T the temperature field, K1, ρ1 and C1 (respectively K2, ρ2 and C2) the scalar thermal conductivity, the mass density and the heat capacity in Ω1 (respectively in Ω2).

Elliptic equation with variable coefficients for the spatial discretization.

−∇. (K1,2∇T) = 0 on Ω1,2/Γ, (4) [T]Γ = T|Ω1 − T|Ω2 = h, (5)

• K ∂T

∂n

• Γ

= K1 ∂T ∂n

• Ω1

− K2 ∂T2 ∂n

• Ω2

= g. (6) where the functions h and g represent the general jump conditions.

CEA, INRIA | September 10-14 2012 | PAGE 2/18

Stefan problem

A non linear model : the unknowns are the temperature field and the location of the interface.

Assumptions :

Density remains constant and the thermophysical properties (conductivities and heat capacities) are different for each phase. No convection in the liquid phase. The interface thickness is null.

Stefan problem

ρsol,liqCsol,liq∂tT − ∇. (Ksol,liq∇T) = 0 on Ωsol,liq/Γ (7) T = Tfusion on Γ, ⇐ Dirichlet boundary condition decoupling the two subdomains (liquid and solid),

• K ∂T

∂n

• Γ = ρL

V . n, ⇐ Stefan condition (energy balance) governing the evolution of the interface, where Tfusion the fusion temperature, L the latent heat and V the interface velocity.

CEA, INRIA | September 10-14 2012 | PAGE 3/18

Outline

1 Elliptic equation with variable coefficients

Discrete form for a control volume Mixed Finite Volume – Finite Element Method Numerical results

2 Stefan problem

Stefan problem algorithm Dirichlet boundary condition Numerical results

3 Conclusion

CEA, INRIA | September 10-14 2012 | PAGE 4/18

Discrete form for a control volume

Equations :

−∇. (K1,2∇T) = 0 on ωi,j ⊂ (Ω1,2/Γ), [T]Γ = T|Ω1

ωi,j − T|Ω2 ωi,j = h,

• K ∂T

∂n

• Γ

= K1 ∂T ∂n

• Ω1

ωi,j

− K2 ∂T2 ∂n

• Ω2

ωi,j

= g.

Discrete form for the control volume ωi,j

−

• m=I,IV
• s=1,2
• eM

s

k1/2∇T.n dl =

• ωi,j

f dS +

• Γ∩ωi,j

g dl, (8) [T]Γ = T|Ω1

ωi,j − T|Ω2 ωi,j = h,

(9) where eM

s , s = 1, 2, the two boundary edges of ωi,j lying inside the dual cell M.

CEA, INRIA | September 10-14 2012 | PAGE 5/18

Mixed Finite Volume – Finite Element Method

Finite Element Approach

Polynomial functions T h are used in each dual cell to approximate the temperature

• field. Those polynomials depend of the four corner values of the dual cell and the

jump conditions across the interface.

⇓

Finite Volume Approach

The flux is analytically calculated on each edges with the polynomial functions.

⇓

✞ ✝ ☎ ✆

Evaluation of

• eM

s k1/2∇T.

n dl, s = 1, 2.

CEA, INRIA | September 10-14 2012 | PAGE 6/18

Without interface case : regular dual cell

Here, the temperature field is approximated by a Q1 polynomial T h : T(x, y) ≈ T h(x, y) = α0 + α1x + α2y + α3xy in ̟M Tcr , r = 1, 4, are the corner values defined at (xcr , ycr ). ̟M the dual cell domain.

Determination of the four unknowns coefficients αγ, γ = 1, 4, of the bilinear representation.

T h (xcr , ycr ) = Tcr r = 1, 4, (10) = ⇒    α0 . . . α3    = A

−1

   Tc1 . . . Tc4    . (11)

Laplacian discretization in the case ∆x = ∆y :

• ωi,j ∇.(k ∇T)dS = k

[

1 4Ti−1,j+1

+

1 2Ti,j+1

+

1 4Ti+1,j+1

+

1 4Ti−1,j

− 3Ti,j +

1 4Ti+1,j

+

1 4Ti−1,j−1

+

1 2Ti,j−1

+

1 4Ti+1,j−1

] . (12)

• E. S¨

uli, Convergence of finite volume schemes for poisson’s equation on nonuniforme meshes, SIAM Journal on Numerical Analysis 28 (5) (1991) 1419-1430

CEA, INRIA | September 10-14 2012 | PAGE 7/18

With interface case : two geometric configurations

Type I interface :

We want ...

the same processing for the 2 geometric configurations. to avoid singularities in the definition

• f the polynomial coefficients.

a continuous processing of the interface during its movement. Type II interface : Here, we choose T h in P2 polynomial space : T(x, y) ≈ T h(x, y) = T h

1 (x, y)

• ̟M Ω1 + T h

2 (x, y)

• ̟M Ω2 where

T h

ς (x, y) = βς 0 + βς 1x + βς 2y + βς 3xy + βς 4x2 + βς 5y 2 , ς = 1, 2.

= ⇒ 12 unknowns coefficients βς

γ, ς = 1, 2, γ = 0, ..., 5.

• M. Oevermann, R. Klein. A cartesian grid finite volume method for elliptic equation with variable coefficients and

embedded interfaces. Journal of Computationnal Physics, 219 :749-769, 2006.

CEA, INRIA | September 10-14 2012 | PAGE 8/18

With interface case : parametrization of the interface

4 relations thanks to the corner values : T h (xcr , ycr ) = Tcr r = 1, 4.

5 relations thanks to the interface parametrization and the jump conditions :

We define s such as

• OM = s

τ , ∀M ∈ Γ. For ς = 1, 2 : T h

ς (s)

= βς

0 +

βς

1

βς

2

t . τs +

• βς

3τxτy + βς 4τ 2 x + βς 5τ 2 y

• s2.

We obtain for the jump conditions ∀s :

• T h

Γ

= T h

1 (s) − T h 2 (s)

• k ∂T h

∂n

• Γ

= k1 ∂T h

1

∂n − k2 ∂T h

2

∂n = a + b s + c s2, = d + e s, where {a, b, c, d, e} are given by the functions h and g (here e is chosen null).

3 closure relations to ensure a continuous processing and avoid singularities :

βς

4 = βς 5,

ς = 1, 2, (13) find

• β1

0, . . . , β1 5, β2 0, . . . , β2 5

• which minimizes
• ̟M

1

• β1

5 +

• ̟M

2

• β2

5,

(14) where ̟M

ς = ̟M Ως, ς = 1, 2.

CEA, INRIA | September 10-14 2012 | PAGE 9/18

Numerical results

We have developed a Finite Volume method with a compact 9-point stencil for 2D Elliptic equations with variable coefficients.

Grid L2 Order L∞ Order 64 × 64 1.0772e-03 4.6612e-04 128 × 128 3.1098e-04 1.79 1.1993e-04 1,96 256 × 256 8.6661e-05 1,84 2.7800e-05 2,11 512 × 512 2.5887e-05 1,74 7.4453e-06 1,90 1024 × 1024 7.9134e-06 1,70 2.3704e-06 1,65 64 × 192 8.5133e-04 5.0004e-04 128 × 384 2.4627e-04 1,79 1.2730e-04 1,97 256 × 768 6.4761e-05 1,92 3.1171e-05 2,03 512 × 1536 1.7885e-05 1,85 8.5367e-06 1,87

Table: Convergence results for the solution u in the L2 and L∞-norm on two different sets of grids with K1 = 10 and K2 = 1.

CEA, INRIA | September 10-14 2012 | PAGE 10/18

Numerical results

Grid K1/K2 = 10−1 Order K1/K2 = 1000 Order 64 × 64 7.4567e-04 2.9247e-03 128 × 128 1.7658e-04 2.07 9.5600e-04 1,61 256 × 256 4.1495e-05 2.09 2.6487e-04 1.85 512 × 512 1.0900e-05 1,93 3.1096e-05 3.09 1024 × 1024 2.7083e-06 2.01 1.1598e-05 1,42

Table: Convergence results for the solution u in the L2-norm with two different ratios K1/K2, K1 = 1

Results...

A compact 9-point method Robustness Second order accuracy

CEA, INRIA | September 10-14 2012 | PAGE 11/18

Outline

1 Elliptic equation with variable coefficients

Discrete form for a control volume Mixed Finite Volume – Finite Element Method Numerical results

2 Stefan problem

Stefan problem algorithm Dirichlet boundary condition Numerical results

3 Conclusion

CEA, INRIA | September 10-14 2012 | PAGE 12/18

Stefan problem algorithm

Initialisation of a time step

Temperature field and interface position.

ր

New temperature fields

On Ωsol,liq/Γ

   ρsol,liqCsol,liq∂tT − ∇. (Ksol,liq∇T) = 0, T = Tfusion on Γ.

Both domains, solid and liquid ones, are solved with a Dirichlet boundary condition on Γ.

տ ց

Interface velocity field

The interface velocity is define by the Stefan condition :

• K ∂T

∂n

• Γ

= ρL V . n.

ւ

Evolving of the interface Γ

The Level Set method is used to evolve the interface.

CEA, INRIA | September 10-14 2012 | PAGE 13/18

Dirichlet boundary condition

Problem to solve :

.

         ∇. (k1∇T) = f1, T = TΓ on Γ, T = w on ∂Ω1/Γ.

Asymptotic equivalent problem :

                         ∇. (k1,2∇T) = f1,2, [T]Γ = TΓ on Γ,

• K ∂T

∂n

• Γ = 0 on Γ,

T = w on ∂Ω1/Γ, T = 0 on ∂Ω2/Γ.

Here f2 = 0 and

✞ ✝ ☎ ✆

k2 = +∞ on Ω2.

− − −− − → − − −− − →

Implementation :

                         ∇. (k1,2∇T) = f1,2, [T]Γ = TΓ on Γ,

• K ∂T

∂n

• Γ = 0 on Γ,

T = w on ∂Ω1/Γ, T = 0 on ∂Ω2/Γ.

Here f2 = 0 and

☛ ✡ ✟ ✠

k2 = k1 min (∆x2, ∆y 2)

• n Ω2.

CEA, INRIA | September 10-14 2012 | PAGE 14/18

Numerical results : Infinite corner problem

Datas :

Tinitial = Tfusion Tfusion = 0 ◦C Tw = 10◦C ρ = 1000 kg/m3 L = 3350000 J/kg kwater = 0.6 W /m/K kice = 2.18 W /m/K Cwater = 4186 J/kg/K Cice = 2260 J/kg/K

Figure: The interface location at different time steps with 35 × 35 grid-points in the (x, y).

Similarity variables : ξx = x √ 4αt , ξy = y √ 4αt with α = kwater Cwaterρ

Figure: The interface location at different time steps with 35 × 35 grid-points in the (ξx, ξy).

CEA, INRIA | September 10-14 2012 | PAGE 15/18

Outline

1 Elliptic equation with variable coefficients

Discrete form for a control volume Mixed Finite Volume – Finite Element Method Numerical results

2 Stefan problem

Stefan problem algorithm Dirichlet boundary condition Numerical results

3 Conclusion

CEA, INRIA | September 10-14 2012 | PAGE 16/18

Conclusion

What have been done ...

A compact 9-point method, A second order accuracy method, An easy way to impose Dirichlet and Neumann boundary conditions, Use the method to solve Stefan problems .

What remains to do ...

= ⇒ realise a coupling between a multi-material solid phase and a liquid phase with convection.

CEA, INRIA | September 10-14 2012 | PAGE 17/18

CEA DAM SSPP LMFA Commissariat ` a l’´ energie atomique et aux ´ energies alternatives Centre de Saclay | 91191 Gif-sur-Yvette Cedex

• T. +33 (0)1 69 08 66 30 | F. +33 (0)1 69 08 66 30

´ Etablissement public ` a caract` ere industriel et commercial RCS Paris B 775 685 019