A Second Order Finite Volume Method for a Multi-material Heat Equation on Cartesian Grids: Application to Stefan problems . ECCOMAS 2012 | Manuel LATIGE , Gerard GALLICE, Thierry COLIN September 10-14 2012
Scientific Background General problem 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. The solid part can be composed by several materials with different properties One of the solids undergoes a phase change (fusion) CEA, INRIA | September 10-14 2012 | PAGE 1/18
Heat equation on multi-material domain Heat equation on multi-material domain ρ 1 , 2 C 1 , 2 ∂ t T − ∇ . ( K 1 , 2 ∇ T ) = 0 on Ω 1 , 2 / Γ . (1) [ T ] Γ = T | Ω 1 − T | Ω 2 = 0 , (2) � � � K ∂ T � = K 1 ∂ T − K 2 ∂ T 2 � � = 0 , (3) � � ∂ n ∂ n ∂ n � � Γ Ω 1 Ω 2 where T the temperature field, K 1 , ρ 1 and C 1 (respectively K 2 , ρ 2 and C 2 ) 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 . where the functions h −∇ . ( K 1 , 2 ∇ T ) = 0 on Ω 1 , 2 / Γ , (4) and g represent the [ T ] Γ = T | Ω 1 − T | Ω 2 = h , (5) general jump conditions. � � � � K ∂ T = K 1 ∂ T − K 2 ∂ T 2 � � = g . (6) � � ∂ n ∂ n ∂ n � � Γ Ω 1 Ω 2 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 , liq C sol , liq ∂ t T − ∇ . ( K sol , liq ∇ T ) = 0 on Ω sol , liq / Γ (7) T = T fusion on Γ, ⇐ Dirichlet boundary condition decoupling the two subdomains (liquid and solid), Γ = ρ L � � K ∂ T � V .� n , ⇐ Stefan condition (energy balance) governing the ∂ n evolution of the interface, where T fusion 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 : −∇ . ( K 1 , 2 ∇ T ) = 0 on ω i , j ⊂ (Ω 1 , 2 / Γ) , [ T ] Γ = T | Ω 1 � ω i , j − T | Ω 2 � ω i , j = h , � � � � K ∂ T = K 1 ∂ T − K 2 ∂ T 2 � � = g . � � ∂ n ∂ n � ω i , j ∂ n � ω i , j � � Γ Ω 1 Ω 2 Discrete form for the control volume ω i , j � � � � � − k 1 / 2 ∇ T . n dl = f dS + g dl , (8) e M ω i , j Γ ∩ ω i , j m = I , IV s =1 , 2 s [ T ] Γ = T | Ω 1 � ω i , j − T | Ω 2 � ω i , j = h , (9) where e M 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 s k 1 / 2 ∇ T .� n dl , s = 1 , 2. e M ✝ ✆ CEA, INRIA | September 10-14 2012 | PAGE 6/18
Without interface case : regular dual cell Here, the temperature field is approximated by a Q 1 polynomial T h : T ( x , y ) ≈ T h ( x , y ) = α 0 + α 1 x + α 2 y + α 3 xy in ̟ M Determination of the four unknowns coefficients α γ , γ = 1 , 4, of the bilinear representation. T h ( x c r , y c r ) = T c r r = 1 , 4 , (10) α 0 T c 1 − 1 . . . . T c r , r = 1 , 4, are the corner = ⇒ = A . (11) . . values defined at ( x c r , y c r ). α 3 T c 4 ̟ M the dual cell domain. Laplacian discretization in the case ∆ x = ∆ y : 1 1 1 � ω i , j ∇ . ( k ∇ T ) dS = k [ 4 T i − 1 , j +1 + 2 T i , j +1 + 4 T i +1 , j +1 + 1 1 . (12) 4 T i − 1 , j − 3 T i , j + 4 T i +1 , j + 1 1 1 4 T i − 1 , j − 1 + 2 T i , j − 1 + 4 T i +1 , j − 1 ] E. S¨ uli , Convergence of finite volume schemes for poisson’s equation on nonuniforme meshes, SIAM Journal on CEA, INRIA | September 10-14 2012 | PAGE 7/18 Numerical Analysis 28 (5) (1991) 1419-1430
With interface case : two geometric configurations We want ... Type I interface : Type II interface : the same processing for the 2 geometric configurations. to avoid singularities in the definition of the polynomial coefficients. a continuous processing of the interface during its movement. Here, we choose T h in P 2 polynomial space : T ( x , y ) ≈ T h ( x , y ) = T h � ̟ M � Ω 1 + T h � 1 ( x , y ) 2 ( x , y ) ̟ M � Ω 2 where � � 4 x 2 + β ς 5 y 2 , ς = 1 , 2. T h ς ( x , y ) = β ς 0 + β ς 1 x + β ς 2 y + β ς 3 xy + β ς ⇒ 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 ( x c r , y c r ) = T c r 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 T h = β ς 1 ς ( s ) 0 + .� τ s β ς 2 β ς 3 τ x τ y + β ς 4 τ 2 x + β ς 5 τ 2 s 2 . � � + y We obtain for the jump conditions ∀ s : k ∂ T h = k 1 ∂ T h ∂ n − k 2 ∂ T h � � 1 2 � T h � = T h 1 ( s ) − T h 2 ( s ) Γ ∂ n ∂ n Γ = a + b s + c s 2 , = 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) � � � � � � β 1 0 , . . . , β 1 5 , β 2 0 , . . . , β 2 � ̟ M � β 1 � ̟ M � β 2 find which minimizes 5 + 5 , (14) � � � � 5 1 2 ς = ̟ M � Ω ς , ς = 1 , 2. where ̟ M 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. L 2 L ∞ Grid Order 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 L 2 and L ∞ -norm on two different sets of grids with K 1 = 10 and K 2 = 1. CEA, INRIA | September 10-14 2012 | PAGE 10/18
Numerical results K 1 / K 2 = 10 − 1 Grid Order K 1 / K 2 = 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 L 2 -norm with two different ratios K 1 / K 2 , K 1 = 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 Interface velocity field On Ω sol , liq / Γ The interface velocity is ρ sol , liq C sol , liq ∂ t T − ∇ . ( K sol , liq ∇ T ) = 0 , define by the Stefan condition : T = T fusion on Γ . Both domains, solid and liquid ones, � � K ∂ T = ρ L � V .� n . are solved with a Dirichlet boundary ∂ n Γ condition on Γ. ւ տ Evolving of the interface Γ The Level Set method is used to evolve the interface. CEA, INRIA | September 10-14 2012 | PAGE 13/18
Recommend
More recommend
