STABILIZED VOLUME AVERAGING FINITE ELEMENT METHOD FOR FLOW IN POROUS MEDIA AND BINARY ALLOY SOLIDIFICATION

1. STABILIZED VOLUME AVERAGING FINITE ELEMENT METHOD FOR FLOW IN POROUS MEDIA AND BINARY ALLOY SOLIDIFICATION Nicholas Zabaras and Deep Samanta Materials Process Design and Control Laboratory Sibley School of Mechanical and Aerospace Engineering188 Frank H. T. Rhodes Hall Cornell University Ithaca, NY 14853-3801 Email: zabaras@cornell.edu URL: http://www.mae.cornell.edu/zabaras/ Materials Process Design and Control Laboratory

2. OUTLINE OF THE PRESENTATION • Overview • Complexities in solidification simulation • Review of the previous related work • Salient features of the current model • Volume averaged governing transport equations • Stabilized equal-order velocity-pressure formulation for porous media flows. • Computational techniques and solution methodology • Applications to combined thermal, fluid and solutal transport in porous media and solidification of a binary mixture • Conclusions Materials Process Design and Control Laboratory

3. OVERVIEW Non-equilibrium effects Mass Transfer Phase Change Solidification Simulation Fluid flow Shrinkage Heat Transfer Deformation Microstructure evolution Materials Process Design and Control Laboratory

4. COMPLEXITIES IN SOLIDIFICATION SIMULATION 1. Morphological and micro-structural complexities 2. Variable time and length scales involved 3. Different physical phenomenon on different scales • fluid flow, convective-conductive heat transfer, • macro-segregation, solid movement and deformation • on macroscopic scale • interdendritic flow, latent heat release, nucleation • and microstructure formation on microscopic scale 4. Formation of two phase mushy zone 5. Diffuse solid-liquid interface difficult to model 6. Non-equilibrium effects on global and local scales Materials Process Design and Control Laboratory

5. Mushy zone liquid ~ 10-2 m solid q (a) Macroscopic scale ~ 10-4 – 10-5m liquid (b) Microscopic scale solid SCHEMATIC OF THE DIFFERENT LENGTH SCALES IN A TYPICAL SOLIDIFICATION PROCESS Materials Process Design and Control Laboratory

6. PREVIOUS WORK ON SOLIDIFICATION AND POROUS MEDIA TRANSPORT SIMULATIONS • Continuum mixture approach (Incropera and co-workers) – Finite difference approach using SIMPLER algorithm • Double diffusive and natural convection in porous media (Nithiarasu et al.) – Volume averaged finite element method (fractional step method for fluid flow problems) • Solidification phenomenon at high Ra number melt flows (Heinrich, Poirier et al.) – Finite element method using penalty based approach for fluid flow • Simulation of darcy flows (Hughes and Masud) – Stabilized finite element approach using mixed and equal order velocity-pressure elements. • Volume-averaged two phase model for solidification and porous media (Beckermann et al.) – Finite difference approach Materials Process Design and Control Laboratory

7. SALIENT FEATURES OF THE CURRENT MODEL 1.Single domain solidification model based onvolume averaging 2.Single set of transport equation for mass, momentum, energy and species transport. 3. No need to track the boundaries between phases 4. Single grid and single set of boundary conditions. 5. Volume averaging tracks solid, mushy and liquid regions 6. Solidification microstructures not modeled. Microscopic governing transport equations General form: Mass : Momentum : Energy : Species : Materials Process Design and Control Laboratory

8. VOLUME AVERAGING PROCESS Volume- averaging process Macroscopic governing equations Microscopic transport equations Volume averaged transport equation Microscopic deviation term = Interfacial flux term due to transport = Interfacial term due to phase change = Materials Process Design and Control Laboratory

9. for both heat and solute transport and and ASSUMPTIONS AND DIMENSIONLESS VARIABLES Assumptions in single-phase model Important dimensionless variables • Microscopic deviation terms neglected Prandtl number Pr • Interfacial terms not modeled in single phase Lewis number Le Darcy number Da • Interfacial drag term in momentum modeled • using Darcy law with isotropic permeability Thermal Rayleigh number RaT Solutal Rayleigh number RaC • Material properties uniform (μ, k etc.)in an • averaging volume dVk but can globally vary Conductivity Ratio Rk • Negligible solutal diffusion in solid phase • (Ds = 0) Capacity Ratio Rc • Solid phase stationary (vs = 0) Materials Process Design and Control Laboratory

10. VOLUME AVERAGED DIMENSIONLESS GOVERNING EQUATIONS Initial conditions : Boundary conditions : Materials Process Design and Control Laboratory

11. NUMERICAL SCHEME FOR FLUID FLOW • Stabilized equal-order velocity-pressure formulation for fluid flow • Derived from SUPG/PSPG formulation • Additional stabilizing term for Darcy drag force incorporated • Darcy stabilizing term necessary for convergence in even basic problems Galerkin formulation for the fluid flow problem Materials Process Design and Control Laboratory

12. NUMERICAL SCHEME FOR FLUID FLOW Stabilized formulation for the fluid flow problem Advection stabilizing term Pressure stabilizing term Diffusion stabilizing term Darcy drag stabilizing term Materials Process Design and Control Laboratory

13. Find such that for all are defined as where and EQUAL ORDER FE FORMULATION Finite element function spaces for velocity and pressure : the following holds : Materials Process Design and Control Laboratory

14. STABILIZING PARAMETERS FOR FLUID FLOW Stabilizing parameters Stabilizing terms advective viscous Darcy pressure continuity • Convective and pressure stabilizing terms • modified form of SUPG/PSPG terms • Darcy stabilizing term obtained by least • squares, necessary for convergence • Viscous term with second derivatives • neglected • A fifth continuity stabilizing term added • to the stabilized formulation Materials Process Design and Control Laboratory

15. for pressure terms modified UNIFIED APPROACH TOWARDS STABILIZING PARAMETERS • Sub-grid scale or multi-scale approach for determining stabilizing parameters • Single expression to account for all effects • Limiting behavior similar to one discussed before Unified expressions for stabilizing parameter Darcy dominant regime Diffusion dominated flows Convection dominated flows for 2D problems including transient effects Materials Process Design and Control Laboratory

16. NUMERICAL SCHEME FOR THERMAL AND SOLUTAL TRANSPORT • Consistent SUPG formulation for thermal and solutal transport • Non-dimensional liquid enthalpy hl expressed as • Gradient of temperature expressed as • Use of consistent instead of lumped matrices Governing Equations Materials Process Design and Control Laboratory

17. NUMERICAL SCHEME FOR THERMAL AND SOLUTAL TRANSPORT Petrov-Galerkin shape function Weak formulation for energy equation: Weak formulation for solute equation: Materials Process Design and Control Laboratory

18. THERMODYNAMIC RELATIONS Non-dimensional parameters Materials Process Design and Control Laboratory

19. ε = and θ = and updating liquid concentration as Cl = Cl = Ce = 0 and ε = • if hSolidus < h < he (region 3), solidification at eutectic point, θ = θeutectic, • Lever rule assumption can be replaced by Scheil rule, ε = UPDATE FORMULAE FOR THERMODYNAMIC QUANTITIES • Supplementary relationships between enthalpy, liquid concentration, volume fraction and temperature • Phase diagram divided into four regions • if h > hLiquidus (region 1), pure liquid region. θ determined from thermodynamic relation, • Cl = 1.0 and ε = 1.0 • if he < h < hLiquidus (region 2), mushy zone, θ and ε determined iteratively by solving • if hSolidus < h (region 4), solid region, θ determined from thermodynamic relations, Cl = Ce and ε = 0.0 Materials Process Design and Control Laboratory

20. SOLUTION METHODOLOGY All fields known at time tn • Multi-step predictor-corrector • method for energy and solute • equations • Standard Gauss-elimination used • for both energy and solute transport n = n +1 Solve for the enthalpy field (energy equation) • Solution typically obtained in • 2 – 3 steps except at initial times Solve for velocity and pressure fields (momentum equation) • Newton – Raphson method for • fluid flow problem Solve for the concentration field (solute equation) • Preconditioned BICGSTAB • algorithm employed for fluid flow Yes • LU factorization is done for few • time steps only Is the error in temperature, liquid concentration and liquid volume fraction less than tolerance Solve for temperature, liquid concentration and volume fraction (Thermodynamic relations) • Line search employed to ensure • global convergence. No Materials Process Design and Control Laboratory

21. NUMERICAL EXAMPLES 1. Double diffusive convection in a constant porosity medium Governing equations Momentum Transport: Energy Transport: Solute Transport: Physical parameters Thermal Rayleigh number, RaT = 2x108 Solutal Rayleigh number, RaC = -1.8x108 Darcy number, Da = 7.407x10-7 Prandtl number, Pr = 1.0 Lewis number, Le = 2.0 Liquid volume fraction, ε = 0.6 Materials Process Design and Control Laboratory

22. TRANSIENT FE SOLUTION (a) Temperature (b) Concentration (c) Streamfunction Materials Process Design and Control Laboratory

23. STEADY-STATE FE SOLUTION (b) Streamlines (a) Finite element mesh - 50x50 Q4 elements (c) Isotherms (d) Iso-concentration lines Materials Process Design and Control Laboratory

24. Width ratio determines extent of porous region 2. Natural convection in a fluid saturated variable porosity medium Governing equations Momentum Transport: Energy Transport: • Central liquid core surrounded by porous medium • Porosity varying uniformly from wall to core Physical parameters Thermal Rayleigh number, RaT = 1x106 Solutal Rayleigh number, RaC = 0 Darcy number, Da = 6.665x10-7 Prandtl number, Pr = 1.0 Wall porosity, εw = 0.4 = 0.3 Width ratio, Materials Process Design and Control Laboratory

25. STEADY-STATE FE SOLUTION (c) Streamlines (b) Isotherms (a) 50x50 finite element mesh (c) Streamlines (b) Isotherms (a) 100x100 finite element mesh Materials Process Design and Control Laboratory

26. TRANSIENT FE SOLUTION (a) Isotherms (b) Streamlines Materials Process Design and Control Laboratory

27. ( = 0.2) STEADY-STATE FE SOLUTION OF THE NATURAL CONVECTION PROBLEM Physical parameters Thermal Rayleigh number, RaT = 1x106 Darcy number, Da = 6.665x10-7 Prandtl number, Pr = 1.0 Wall porosity, εw = 0.4 = 0.2 Width ratio, (c) Streamlines (b) Isotherms (a) 80x80 finite element mesh Materials Process Design and Control Laboratory

28. 3. Solidification of a binary aqueous solution Physical parameters Thermal Rayleigh number, RaT = 1.938x107 Solutal Rayleigh number, RaC = -2.514x107 Prandtl number, Pr = 9.025 Heat conductivity ratio, Rk = 0.84 Heat capacity ratio, Rc = 0.576 Dimension slope of liquidus, m = 0.905 Darcy number, Da = 8.896x10-8 Lewis number, Le = 27.84 Initial and boundary conditions Temperature of hot wall, Thot = 311 K Temperature of cold wall, Tcold = 223K Initial temperature, T0 = 311K Initial concentration, C0 = 0.7 • Presence of a dendritic mushy zone during solidification Eutectic concentration, Ce = 0.803 • Thermal and solutal buoyancy forces opposing each other Eutectic temperature, Te = 257.75K • Highly irregular liquid interface Solutal flux on all boundaries = 0 (adiabatic flux condition) • Temperature and concentration fields distorted by • advection effects. Materials Process Design and Control Laboratory

29. FE SOLUTION AT DIFFERENT TIME STEPS (t* = 0.009) (c) Isotherms (d) Liquid concentration (b) Streamlines (a) Velocity vectors Materials Process Design and Control Laboratory

30. FE SOLUTION AT DIFFERENT TIME STEPS (t* = 0.018) (c) Isotherms (d) Liquid concentration (b) Streamlines (a) Velocity vectors Materials Process Design and Control Laboratory

31. FE SOLUTION AT DIFFERENT TIME STEPS (t* = 0.036) (c) Isotherms (d) Liquid concentration (b) Streamlines (a) Velocity vectors Materials Process Design and Control Laboratory

32. FE SOLUTION AT DIFFERENT TIME STEPS (t* = 0.071) (c) Isotherms (d) Liquid concentration (b) Streamlines (a) Velocity vectors Materials Process Design and Control Laboratory

33. FE SOLUTION AT DIFFERENT TIME STEPS (t* = 0.142) (c) Isotherms (d) Liquid concentration (b) Streamlines (a) Velocity vectors Materials Process Design and Control Laboratory

34. NUMERICAL COMPARISON AND MACROSEGREGATION Macrosegregation patterns at steady state (t = 0.142) Materials Process Design and Control Laboratory

35. TRANSIENT FE SOLUTION OF THE SOLIDIFICATION PROBLEM (a) Temperature (b) Liquid concentration (c) Liquid volume-fraction (d) Streamfunction Materials Process Design and Control Laboratory

36. CONVERGENCE STATISTICS FOR FLUID FLOW SOLVER Residual norm Iteration number Example 3 Example 2 Example 1 Materials Process Design and Control Laboratory

37. CONCLUSIONS AND FUTURE RESEARCH PLANS CURRENT ACHIEVEMENTS • Stabilized formulation for high Rayleigh number porous media and solidification problems • Single domain, enthalpy based models incorporating mushy zone • Porous media simulations with convergence studies • Solidification of a binary aqueous solution at high thermal and solutal Rayleigh numbers FUTURE RESEARCH INTERESTS • Incorporation of a unified stabilizing parameter for all regimes • 3D solidification simulation and development of a parallel simulator • Solidification of metal alloys with phase constituents of different densities • Sensitivity analysis and parametric studies • Magneto-convection and solidification in the presence of magnetic fields. • Marangoni convection and effects of thermo-capillary convection on solidification • Combined analysis of solidification and deformation, and modeling of residual stresses • Inverse design and control of solidification processes with mushy zones Materials Process Design and Control Laboratory