The Stefan Problem and the Exact Solution for the Two Phase Stefan Problem

1. The Stefan Problemand the Exact Solution for the Two Phase Stefan Problem DT Project

2. First Things First • We desire to formulate Two-Phase Stefan Model of Melting and Freezing • Non linear model • “Moving Boundary Problem” • Contains an unknown which is the region to be solved…..

3. Why Stefan Problem??? • The formulation of the Stefan Problem is a foundation on which more complex models can be built.

4. First (cont’d) • Note: • The phase changing process is governed by the conservation of energy. • The unknowns are the temperature field and the location of the Interface. • Involves a phase change material(PCM) with constant density( ), latent heat( ), melt temperature( ), Specific heats( ), and Thermal conductivities( ).

5. Physical Assumptions • Conduction only • Constant latent heat (L) • Fixed melting temperature( ) which is according to the phase change material (PCM) • Interface thickness is 0 and it is a sharp front; it separates the phases

6. Assumptions (cont’d) • Thermophysical properties are different for each phase • Conductivities ( ) • Specific heats ( ) • Density remains constant

7. Other Assumptions • Nucleation and supercooling are assumed to be not present • Surface tension and curvature is insignificant

8. Only Conduction • Conduction of Heat • Temperature • Heat(enthalphy) • Heat Flux • Characterizes phases

9. Heat Equation -Heat conduction equation (one space dimension): -well-posed -Heat equation: where

10. The Two-Phase Stefan Problem • A slab, , initially solid at temperature , is melted by imposing a hot temperature at the face and keeping the back face, insulated (all parameters constant).

11. Let’s Find a Solution • The solution of the Stefan Problem is T(x,t) and X(t)!!!!

12. Mathematical Model • PDE for • Interface(t>0) • Initial Condition • Boundary Condition

13. Be Exact! • In order to explicitly solve the Two-Phase Problem we need to assume the slab is semi-infinite. • Physical Problem: • We want to melt a semi-infinite slab, , initially solid at a temperature , by imposing temperature , on the face . • The “alphas” are different for each phase. All parameters constant.

14. Mathematical Model • Heat Equations for • Interface(t>0) • Stefan condtion: • Initial Condition • Boundary Condition

15. Two-Phase Neumann Solution • We derive the Neumann Solution • We use the similarity variable, ,and seek the solution for both for the liquid and for the solid. • Seek a solution for X(t) in the form:

16. Temperature • Temperature in the liquid region at t>0: • Temperature in the solid region at t>0:

17. Neumann similarity solution of the 2-phase Stefan Problem for the interface

18. Transcendental Equation • There is a different Transcendental Equations that incorporates the TWO Stefan numbers (one for each phase). • Trans. Equation: • Stefan Numbers and parameter “v”:

19. Newton!!! • Newton's method is an algorithm that finds the root of a given function. -Where F(x) = 0. • The fastest way to approximate a root.

20. What now??? • We can use the Newton algorithm to solve the transcendental equation for the Stefan problem • We do this in order to find a unique root lambda and therefore a unique similarity solution for each

21. Newton Algorithm • Guess x0 • Take a Newton Step where 3. Terminate if -(Want TOL to be very small so convergence should be noticed) -Stuck. Better guess

22. Approximate the Root for Newton • For the 2-phase problem there is a good approximation which can be x0 in the Newton program. Approximation of Lambda:

23. Glauber’s Salt Input Values -Tm = 32; -CpL = 3.31; -CpS = 1.76; -kL = .59e-3; -kS = 2.16e-3; -rho = 1460; -Lat = 251.21; -Tinit = 25; -Tbdy = 90;

24. Water Input Values - Tm = 0; - CpL = 4.1868; - CpS = .5; - kL = .5664e-3; - kS = 2.16e-3; - rho = 1; - Lat = 333.4;

25. Glauber’s salt Example maximum number of iterations to be performed,20 tolerance for the residual,1.0e-7 Iterations is 1 xn = 5.170729e-001 fx = -2.740474e-001 Iterations is 2 xn = 5.207715e-001 fx = 1.747621e-002 Iterations is 3 xn = 5.207862e-001 fx = 6.881053e-005 Done. Root is x=5.207862e-001, with Fx=1.068988e-009, Iterations is n=4 LAMBDA= • 5.207861719133929e-001

26. Plot Lambda vs Stefan For small Stefan numbers from 0:5 and M=51

27. Portion of Exact Solution Code • function lambda = neumann2p(CpL, CpS, kL, kS, rho, Tm, Lat, Tbdy,Tinit) • format long e • %----------------------Input values--------------------------------------- • CpL = input('Enter specific heat: '); • CpS = input('Enter specific heat of solid: '); • kL = input('Enter thermal conductivity of liquid: '); • kS = input('Enter themal conductivity of solid: '); • rho = input('Enter density which is constant: '); • Tm = input('Enter the melting temperature of substance: '); • Lat = input('Enter Latent heat: '); • TL = input('Enter temperature at x=0: '); • dat2p; • %--------------------------Constants to derive----------------------------- • alphaL = kL/(rho*CpL); • alphaS = kS/ (rho *CpS); • dT = Tbdy - Tm; • dTa = Tm - Tinit; • v = sqrt(alphaL./alphaS); StS = (CpS*dTa)/ Lat; • StL = (CpL*dT)/Lat; • %----------------------------Newton----------------------------------- • x0 = 0.5 * (-StS/ (v*sqrt(pi)) + sqrt(2*StL + (StS/(v*sqrt(pi)))^2)); • lambda = transnewton2p(x0,20,1.0e-7,StS,StL,v); • function Xt = XofT(lambda,alphaL,t) • Xt = 2*lambda*sqrt(alphaL*t); • function TofXTL = TofXTL(lambda,alphaL,Tbdy,dT,x,t) • TofXTL = Tbdy-dT*(erf(x./(2*sqrt(alphaL*t)))./erf(lambda)); • function TofXTS = TofXTS(lambda,alphaS,Tinit,dTa,x,t,v) • TofXTS = Tinit +dTa*(erfc(x./(2*sqrt(alphaS*t)))./erfc(lambda*v))

28. Exact Front for Glauber’s Salt

29. Exact Front for Water

30. Exact Histories for Salt -5,10,15 from top to bottom

31. Exact Histories for Water

32. Exact Profiles for Salt -notice sharp turn at the melt temperature

33. Exact Profiles for Water

34. Enthalpy Method for Stefan 2-Phase Problem

35. Formulation/ Discretization of Stefan 2-Phase Problem • Discretize Control Mesh • Discretize Heat Balance • Discretize Fluxes • Discretize Boundary Conditions

36. Partition Control Volumes: Subdivide the regions into M intervals, or control volumes: with each Subregion associate a node . Volume of :

37. Create the control mesh:

38. Discretize Heat Balance Integrating heat balance eq’n over control volume and over time interval

39. Continue… Dividing out the A and integrating the derivatives yields: Assuming E is uniform and is small:

40. Discretized Enthalpy From Last Slide: Or

41. Discretize Fluxes Fourier’s Law: Approximate q discretely:

42. Discretize Boundary Conditions • Imposing Temperature from left: • Impose exact temperature at back face. • The left boundary flux at x=0: • The right boundary flux at x=l:

43. Liquid Fraction & “Mushy” •  is Solid •  is liquid •  is “mushy” Liquid Fraction:

44. Energy & Temperature Relation Solve for T:

45. Energy vs. Temperature Graph with Enthalpy Scheme:

46. Explicit Time Scheme Initial Temperature Known: • Initial Enthalpy Set Resistance and Fluxes from Initial Temperature and Enthalpy Update Enthalpies at Update Temperature Update Liquid Fraction

47. Flux and Resistance Drive Heat Flows

48. Update Enthalpy at Next Time Step(The Discretized Enthalpy)

49. Update Temperature at Next Time Step:

50. Update Liquid Fractions/Phases: