Melting Ice

Numerical methods and mathematical models of Elmer
tetraeder
Posts: 44
Joined: 27 Sep 2012, 13:09
Antispam: Yes

Melting Ice

Post by tetraeder »

Hello community,

First i would like to say this is a really great forum - most of my "basic" questions have been already answered by other threads - Thanks for that!

I have the following problem:
We've done the following experiment in our labor: Consider a block of ice. Now you put another body (lets say it's a block of copper) with an internal heat source on top of the block of ice.

Now i would like to simulate this experiment with elmer. But with my previous tests (Steady State) i get too low temperatures inside the block of copper (experimental steady state: 30°C ; simulation steady state 4°C). Maybe it's because i've only considered conduction, no phase change of ice and no convection of air on top of the block of copper.
I've done a sketch in 2D, but the simulation is in 3D. The annotations are the BCs.

So there are some Questions about my idealization:

1) does the implementation of phase change significantly change my simulation results?

2) is it necessary to use transient simulation, when i implement phase change by spatial 2 (the docs mention 3 types for phase change: spatial 1, spatial 2, temporal - and it only says that temporal is needed for transient simulation).

3) does the implementation of convection significantly change my simulation results?


I hope, that someone can answer my questions.

Regards
Attachments
Ohne Titel.png
Ohne Titel.png (609.69 KiB) Viewed 7451 times
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Melting Ice

Post by mzenker »

Hi,

my guesses would be

1) Yes

2) Yes

3) Maybe

HTH,

Matthias
tetraeder
Posts: 44
Joined: 27 Sep 2012, 13:09
Antispam: Yes

Re: Melting Ice

Post by tetraeder »

Thank you Matthias.

OK when i implement transient simulation, what about the enthalpy? I set the enthalpy of ice (330000 J/kg), but there's a warning, that elmer also needs the enthalpy of all the other materials. Maybe i'm completely wrong with my .sif ...

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 5
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Transient
  Steady State Max Iterations = 1
  Output Intervals = 20
  Timestepping Method = BDF
  BDF Order = 1
  Timestep intervals = 20
  Timestep Sizes = 0.5
  Solver Input File = case.sif
  Post File = case.ep
End

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.67e-08
  Permittivity of Vacuum = 8.8542e-12
  Boltzmann Constant = 1.3807e-23
  Unit Charge = 1.602e-19
End

Body 1
  Target Bodies(1) = 1
  Name = "Body 1"
  Equation = 1
  Material = 1
  Initial condition = 2
End

Body 2
  Target Bodies(1) = 2
  Name = "Body 2"
  Equation = 1
  Material = 3
  Initial condition = 2
End

Body 3
  Target Bodies(1) = 3
  Name = "Body 3"
  Equation = 1
  Material = 2
  Initial condition = 1
End

Body 4
  Target Bodies(1) = 4
  Name = "Body 4"
  Equation = 1
  Material = 1
  Initial condition = 2
End

Solver 1
  Equation = Heat Equation
  Procedure = "HeatSolve" "HeatSolver"
  Variable = -dofs 1 Temperature
  Exec Solver = Always
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-10
  Linear System Preconditioning = ILU0
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 1
  Linear System Precondition Recompute = 1
End

Solver 2
  Equation = Result Output
  Output Format = Vtu
  Output File Name = case
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Exec Solver = Always
End

Equation 1
  Name = "heatEq"
  Phase Change Model = Spatial 2
  Active Solvers(1) = 1
End

Equation 2
  Name = "Output"
  Active Solvers(1) = 2
End

Material 1
  Name = "Copper (generic)"
  Heat expansion Coefficient = 16.5e-6
  Heat Conductivity = 401.0
  Sound speed = 3810.0
  Heat Capacity = 385.0
  Mesh Poisson ratio = 0.34
  Density = 8960.0
  Poisson ratio = 0.34
  Youngs modulus = 115.0e9
End

Material 2
  Name = "Water (frozen in Celcius)"
  Viscosity Model = Power law
  Enthalpy = 333000.0
  Viscosity = Variable Temperature; Real MATC "(2.0*3.0*1.916E03 * exp( -139.0E03/(8.314 *(tx+273.15))))^(-1.0/3.0)
  Critical Shear Rate = 1.0E-6
  Heat Conductivity = Variable Temperature; Real MATC "9.828*exp(-5.7E-03*(tx+273.15))"
  Heat Capacity = Variable Temperature; Real MATC "146.3+(7.253*(tx+273.15))"
  Viscosity Exponent = $ (1.0/3.0)
  Density = 910.0
End

Material 3
  Name = "Epoxid"
  Heat Conductivity = 0.5
  Heat Capacity = 1000.0
  Density = 1200.0
End

Initial Condition 1
  Name = "InitialIceTemp"
  Temperature = -20.0
End

Initial Condition 2
  Name = "InitialProbeTemp"
  Temperature = 20.0
End

Boundary Condition 1
  Target Boundaries(1) = 2
  Name = "iceTemperature_contactingSurface"
  Temperature = 0.0
End

Boundary Condition 2
  Target Boundaries(1) = 3
  Name = "iceTemperature_onTop"
  Temperature = 0.0
End

Boundary Condition 3
  Target Boundaries(1) = 4
  Name = "iceTemperature"
  Temperature = -20.0
End

Boundary Condition 4
  Target Boundaries(1) = 5
  Name = "heatSource_1"
  Heat Flux = 149605.6465
End

Boundary Condition 5
  Target Boundaries(1) = 6
  Name = "heatSource_2"
  Heat Flux = 149605.6465
End
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Melting Ice

Post by mzenker »

Hi,

you need the enthalpy of all materials for which the equation with phase change applies.
If there are materials where you don't want phase change, you can define another equation for them:

Code: Select all

Body 1
  Target Bodies(1) = 1
  Name = "Body with phase change"
  Equation = 1
...
End

Body 2
  Target Bodies(1) = 1
  Name = "Body without phase change"
  Equation = 2
...
End

...

Equation 1
  Name = "heatEqwithPhase"
  Phase Change Model = Spatial 2
  Active Solvers(1) = 1
End

Equation 2
  Name = "heatEqwithoutPhase"
  Active Solvers(1) = 1
End

...
BTW, the enthalpy and its units are discussed here: viewtopic.php?f=3&t=2345.

HTH,

Matthias
tetraeder
Posts: 44
Joined: 27 Sep 2012, 13:09
Antispam: Yes

Re: Melting Ice

Post by tetraeder »

Thank you so much! Using two equations works fine. And thank you for the link. Now i'll try to get the right enthalpy work. If it works i'll let you know.

regards
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Melting Ice

Post by raback »

Hi

Just a quick comment on the physics: I would be surpriced if you'll get a factor of 8 by considering the phase change. Your contact temperature with the water should be pretty close to 0 still. Of course this may raise a little bit the temperature. However, your dT over the block is only 4 degs and this will not much be affected.

I would look for a possible explanation to the discrepancy rather in the given heat fluxes. Are you sure that that is what you get into the system? Is the contact good? Or do you rather have a internal heating by electric current? If so, I would first make a proper model for that neglecting the ice in the fist place. The phase change has only an effect in transient and if you reach steady state you could do without.

-Peter
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Melting Ice

Post by mzenker »

Fortunately I declared my guesses above as such - I admit that I thought about the physics too shortly ... ;)

Matthias
tetraeder
Posts: 44
Joined: 27 Sep 2012, 13:09
Antispam: Yes

Re: Melting Ice

Post by tetraeder »

Hey,
today i did the transient simulation and unfortunately there's only a small increase of temperature (about 8 °C)... but i did not reach steady state, so i think it will be nearly the same. But a transient simulation will be a next task anyway, so your answers still helped me - thank you.

The heating is by two electrical heat units (cylindrical: 8mm diameter; 50mm length) and each of them provide a power of 200 W around the lateral surface. The holes lengths are 47mm, so i calculated the effective heat flux to 149605.6465 W/m^2 instead of 159154.9431 W/m^2. But even i use 159154 W/m^2 there will be no significant changes in Temperature. In the experiment, the two heat sources were used with thermal compound (Wärmeleitpaste) to provide a good heat transfer into the block of copper.
Epy
Posts: 11
Joined: 12 Nov 2012, 02:29
Antispam: Yes
Location: CA, USA
Contact:

Re: Melting Ice

Post by Epy »

Hi,

Could you explain your setup a little more clearly? In your picture, are the temperatures you have indicated held constant, or is that what they start out at? When your experiment finishes, is the block of ice still ice? I'm going to assume that the true BCs here are the constant heat flux applied to the holes and the what I assume is ambient temperature of -20 C. The 0 C at the edge of copper comes from the assumption that ice is still touching copper and not water, yes?

Due to copper being an excellent heat conductor, you can assume that the temperature of the block is mostly uniform--the differences will be in the water/ice region. As Peter mentioned, phase change will only affect the transient, not the steady state.

My thoughts:
1) Once/if ice is melted at the ice/copper interface, you'll get water (fluid) and that situation will change from conduction to free convection. You'll then have: (a) convection between copper and air, (b) convection between copper and water layer, (c) convection between water layer outer surface and air, (d) conduction through water layer, (e) convection between water and ice, (f) conduction through ice block, and (g) convection between ice block and air.
2) From 1), if you're getting such a situation, it's possible that the majority of heat is between transferred between copper and air, and copper->water->water outer surface and the ice is mostly untouched. This is likely because the thermal conductivity of water is less than ice, meaning that heat builds up and stays in the water layer, especially with free (not forced) convection.
3) You say ice, but is it ice with air bubbles frozen in? Ice with no air bubbles? Depending on the amount of air present, the heat conductivity will decrease (since air has much less heat conductivity than ice).

It's my guess that:
a) there's a small water layer between the copper and ice, acting as an insulating layer and b) your ice has air bubbles in it, making it more insulating than pure ice. These would both increase the steady-state temp of the copper.
tetraeder
Posts: 44
Joined: 27 Sep 2012, 13:09
Antispam: Yes

Re: Melting Ice

Post by tetraeder »

Thank you Epy

okay i will try to explain the setup in more detail. We've started freezing a block of ice to a temperature of -20°C. Then we've put the block of copper on the upper surface to start the melting process (penetration) by supply the heat sources with electric current. But we've had to switch them off 3 times, because we were afraid of destroying the heat sources by the high temperatures. After the block of copper reached a certain depth, there was no more need to switch them off. Because there were cables on the upper part of the block of copper, we finished the experiment at a depth of 30 mm.

So what i want to do is a steady state "snapshot" at a depth of 30 mm.

I was hoping that there's a workaround of using this thin water layer, but apparently i should implement it. I guess i can simulate the air bubbles by changing the heat conductivity of the ice. But what is the best way to simulate that thin layer of water? Do i have to create a new mesh and define the water as a new body. What would be the thickness at ambient pressure?
Attachments
Ohne Titel.png
Ohne Titel.png (855.22 KiB) Viewed 7422 times
Post Reply