Re: stefan problem level set implementation

[Daniel Wheeler]

Chris, You need to check out trunk to run this. Also, can we take this discussion over to the list. Updated Stefan problem:

#!/usr/bin/env python

from fipy import *


boxSize = 1.
nx = 100
dx = boxSize / nx
Kliquid = 1
Ksolid = 1.4
cpliquid = 1
cpsolid = 1
density = 1
latentHeat = 1
liquidBoundTemp = 1
solidBoundTemp = 0
steps = 300
tolerance = 0.1
dt = Variable(1.0)

mesh = Grid1D(nx=nx, dx=dx)
X = mesh.getCellCenters()[...,0]

temperature = CellVariable(name="Temperature", mesh=mesh, hasOld=1)
extensionVelocity = CellVariable(name="Extension Velocity", mesh=mesh)
distanceVar = DistanceVariable(name='level set variable', mesh=mesh, value=-1, hasOld=1)
distanceVar.setValue(1, where=X > boxSize / 2.)
distanceVar.calcDistanceFunction()

K = Ksolid + (distanceVar < 0) * (Kliquid - Ksolid)
cp = cpsolid + (distanceVar < 0) * (cpliquid - cpsolid)

def faceFunc(cell1, cell2, disVar=None):
    value = numerix.minimum(cell1, cell2)
    id1, id2 = mesh._getAdjacentCellIDs()
    disVar1 = numerix.take(disVar,id1)
    disVar2 = numerix.take(disVar,id2)
    return numerix.where(disVar1 * disVar2 < 0, 0, value)
    
Kface1 = K.getFaceValue(faceFunc, disVar=distanceVar)
Kface2 = K.getFaceValue(lambda a, b: numerix.minimum(a,b))

##velocity =  (Kface2 * temperature.getFaceGrad()).getDivergence() * dx/ latentHeat
velocity =  - K * temperature.getGrad().dot([1,]) / latentHeat

advectionEq = buildHigherOrderAdvectionEquation(advectionCoeff=extensionVelocity)

from fipy.viewers import make
viewer = make(vars=(distanceVar,temperature), limits={'datamin':-1., 'datamax':1.})

class CellInterfaceAreas(CellVariable):
    def __init__(self, distanceVar):
        CellVariable.__init__(self, mesh=distanceVar.getMesh())
        self.distanceVar = self._requires(distanceVar)

    def _calcValue(self):
        return self.distanceVar.getCellInterfaceAreas()

interfaceAreas = CellInterfaceAreas(distanceVar)

sourceVelocity = (distanceVar - distanceVar.getOld()) / dt / distanceVar.getGrad().getMag()
temperatureEq = TransientTerm(density * cp) == ImplicitDiffusionTerm(Kface1) +  latentHeat * cp * velocity * interfaceAreas / mesh.getCellVolumes()

from fipy.boundaryConditions.fixedValue import FixedValue
BCs=(FixedValue(mesh.getFacesLeft(),solidBoundTemp),FixedValue(mesh.getFacesRight(),liquidBoundTemp))

dt.setValue(1.e-5)
for step in range(steps):
    distanceVar.updateOld()
    temperature.updateOld()
    advectionRes = 1e+10
    temperatureRes = 1e+10
    totalSweeps = 20
    sweep = 0
    print 'step',step
    dt.setValue(dt * 2.0)
##    while max(advectionRes, temperatureRes) > tolerance and sweep < totalSweeps:
    extensionVelocity.setValue(velocity)
    distanceVar.extendVariable(extensionVelocity)
## 	print
## 	print 'extensionVelocity',extensionVelocity
    print 'dt',dt
    print 'max(extensionVelocity)',max(extensionVelocity)
    print 'min(extensionVelocity)',min(extensionVelocity)
    dt.setValue(min(dx / max(abs(extensionVelocity + 1e-10)) / 10, dt.getValue()))
    advectionEq.solve(distanceVar, dt=dt)
    temperatureEq.solve(temperature, dt=dt, boundaryConditions=BCs)
##    raw_input()
##    print 'advectionRes',advectionRes
##    print 'temperatureRes',temperatureRes
## 	print 'velocity',velocity
## 	print 'interfaceAreas',interfaceAreas
## 	print 'temperature',temperature
## 	print 'distanceVar',distanceVar
##        sweep += 1     


    viewer.plot()


raw_input()

Cheers On Jan 22, 2007, at 9:52 PM, Chris Calebrese wrote: Hello,   I’m working with Dan Lewis at RPI and am currently trying to solve a simple 1D Stefan problem using the level set method to account for interface movement.  The eventual goal is to build up a bubble growth model that will be more complex (heat flow, mass flow, exothermic reaction, etc), but right now I am working understanding how to solve a simpler problem so that the framework for the more complex problem is sound.   The problem I am looking at is the melting or solidification of “ice” (T=0 C at the interface, no density change in during the phase transition, heat flow only).  My difficulty lies in maintaining the interface temperature at 0 C.  The heat flow equation I am solving has a source term that evaluates to zero everywhere except at the interface (with a width 2* epsilon, where epsilon is spacing distance).  The velocity of the interface, the delta function for the source term, and the heat flow equation are below.   Rho = density k = thermal conducitivity/(heat capacity * density) L = heat of fusion     from Numeric import cos delta = CellVariable(name = "delta", mesh = mesh, value=0) delta.setValue((1 + cos(3.14159265*distVar()/epsilon))/(2*epsilon), where=((distVar() > -epsilon) & (distVar() < epsilon)))   velocity=-Temp.getGrad()[i-1] * kliquid/(rho * L)+Temp.getGrad()[i]*ksolid/(rho *L) (The index i is the index of the cell closest to the interface in the solid, and i-1 the index of the cell closest to the interface in the liquid.)   eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=DiffusionCoeff)-velocity*L*delta/cp   L is the latent heat of fusion, and cp is the heat capacity (taken as the same value in the solid and liquid for now).  The coefficient for the diffusion term is dependent on the phase of the material.  The source term value is set to absorb the amount of heat absorbed by melting of the interface in a given time step.     When I run the code the source term will absorb heat, but the temperature of the interface rises above T=0, and continues to do so as the script runs.  I’m stuck right now as how to maintain this internal boundary condition without artificially resetting the interface temperature after each time step.  Could you offer any advice on how to implement this so the correct behavior is maintainted?  If you need more info on solution method or code I can provide it.  Thank you for your help.   Chris Calebrese -- Daniel Wheeler