Hi,
My ultimate goal: a numerical model very similar to doi:10.2298/SOS1001033N (Computer simulation of rapid solidification with undercooling: A case study of spherical ceramics sample on metallic substrate). The full text pdf of that paper is easily found online. I'm interested in solidification of a pure metal with interface velocity a function of interface temperature (thermally undercooled solidification of a pure metal). My knowledge of solidification, Python, and FiPy is quite basic. I just discovered FiPy a few weeks ago. I want to start by considering the level set method in FiPy. Start with the basic 1D Stefan problem: solidification of a pure material on a finite domain. I know nothing about the phase field method and whether it is applicable to the ultimate goal in paragraph 1, in which interface velocity is a specified function of interface temperature. I found a thread from 2007 https://www.mail-archive.com/fipy@nist.gov/msg00320.html, but the code of Daniel Wheeler is too advanced for me to understand, and does not run unmodified in current FiPy. Without being able to comprehend this code in 1D, extension to 2D cylindrical coordinates seems unattainable. I looked at the level set example of a 'simple trench system:' http://www.ctcms.nist.gov/fipy/examples/levelSet/generated/examples.levelSet.electroChem.howToWriteAScript.html I tried to imitate this example (the metal ion Cm part) with the scripts below, but am stuck as to how to couple the interface velocity into my problem as an unknown. In my script below, the interface velocity appears to remain at the initial specified value, since it apparently never gets updated or coupled into the system. It would seem necessary to have the Stefan condition as the coupling equation, but I am unable to look at the FiPy manual and FiPy Python code and identify the necessary commands. I looked for the equivalent of the Stefan condition (not expressed as a source term, but expressed in terms of gradients of metal ion concentration on either side of the interface) in the code for the example Simple Trench System, but was unable to identify it. The Stefan condition is currently represented in my script as a source term, but I do not see how to couple Vinterface into the unknowns. It seems like one approach could be to replace Vinterface in the source term with an expression involving gradients of temperature on either side of the interface, but I am stuck as to how to write this in the language of FiPy in both 1D and eventually 2D (first Cartesian 2D then cylindrical 2D). I would love to see code that is can be obviously generalized to 2D Cartesian then 2D Cylindrical, instead of being specialized only to 1D. The code in https://www.mail-archive.com/fipy@nist.gov/msg00320.html might be 1D-only and/or non-obvious to translate to 2D, for example. Thanks MY MAIN PYTHON SCRIPT # from IPython import embed #then put embed() in your script at a point where you want to be at ipython prompt # take the approach of having variables with conventional names e.g. temperature, then set to 1.0 etc for nondimensional # this lets you choose different nondimensionalization schemes # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( nx = 20 L=1. cellSize = L/nx #there is no use for a graded mesh since grid is fixed and interface is moving mesh = Grid1D(nx=nx, Lx=L) x=mesh.cellCenters # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( numberOfCellsInNarrowBand = 10 narrowBandWidth = numberOfCellsInNarrowBand * cellSize cflNumber=.2 # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( #level set method distance variable, initialized to -1. distanceVar = DistanceVariable(name = 'distance variable',mesh = mesh,value = -1.,hasOld = 1) #TODO I think the interface has to advance a little into the domain to establish where it is? #alternative is to have extra cells to the left of the actual x=0, as in simpleTrenchSystem.py distanceVar.setValue(1., where=(x > cellSize)) distanceVar.calcDistanceFunction(order=2) # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( #thermal diffusivity alphaSolid=1. alphaLiquid=1. # #thermal conductivity # kSolid=1. # kLiquid=1. # #mass density # rhoSolid=1. # rhoLiquid=1. #specific heat cpSolid=1. cpLiquid=1. latentHeat=1. Tm=.5 #melting T temperature = CellVariable(name="temperature",mesh=mesh,value=1.,hasOld=False) #left boundary is step change in temperature valueLeft=0. temperature.constrain(valueLeft, mesh.facesLeft) #right boundary is zero flux, which is default BC in FiPy # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( #TODO this scheme is not going to work. interface velocity starts out with initial guess. how is it updated or linked to anything? #simple trench electro problem: interface velocity was a function of metal ion bulk concentration and theta interfaceVelocVar=CellVariable(name='interface velocity',mesh=mesh,value=0.1) #TODO initializing V at zero results in errors advectionEquation = TransientTerm() + AdvectionTerm(interfaceVelocVar) # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( energyEquation=buildEnergyEquation(temperature=temperature,distanceVar=distanceVar,interfaceVelocVar=interfaceVelocVar,cellSize=cellSize,alphaSolid=alphaSolid,alphaLiquid=alphaLiquid,latentHeat=latentHeat,cpSolid=cpSolid,cpLiquid=cpLiquid,Tm=Tm) # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( viewer = Viewer(vars=(temperature,),datamin=0., datamax=1.) viewer.plot() raw_input() levelSetUpdateFrequency = int(0.8 * narrowBandWidth / (cellSize * cflNumber * 2)) #TODO where does this come from? numberOfSteps=80 for step in range(numberOfSteps): if step>5 and step % 5 == 0 and viewer is not None: viewer.plot() raw_input() if step % levelSetUpdateFrequency == 0: distanceVar.calcDistanceFunction(order=2) distanceVar.updateOld() distanceVar.extendVariable(interfaceVelocVar, order=2) dt = cflNumber * cellSize / interfaceVelocVar.max() advectionEquation.solve(distanceVar, dt = dt) energyEquation.solve(temperature,dt=dt) MY SUBROUTINE SCRIPT TO BUILD ENERGY EQUATION from fipy.terms.implicitSourceTerm import ImplicitSourceTerm from fipy.terms.transientTerm import TransientTerm from fipy import FaceVariable, CellVariable, DiffusionTerm def buildEnergyEquation(temperature,distanceVar,interfaceVelocVar,cellSize,alphaSolid=1.0,alphaLiquid=1.0,latentHeat=1.0,cpSolid=1.,cpLiquid=1.,Tm=0.): r''' #based on metalIonDiffusionEquation.py from NIST #solve the one dimensional solidification problem from materials processing class, #9/28/1999 lecture: solidification in finite domain; page 105 in class notes #level set method ''' # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( mesh = distanceVar.mesh # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( # http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html #trying to adapt the position-dependent diffusivity example to distanceVar-dependent diffusivity D = FaceVariable(mesh=mesh, value=alphaSolid) #https://www.mail-archive.com/fipy@nist.gov/msg00854.html distVarAtFaceCenters=distanceVar.getHarmonicFaceValue() D.setValue(alphaLiquid, where=(distVarAtFaceCenters > 0.)) # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( cp=CellVariable(name='specific heat',mesh=mesh,value=cpSolid,hasOld=0) cp.setValue(cpLiquid, where=(distanceVar > 0)) # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( #TODO is it OK that interfaceVelocVar is a CellVariable? in metalIonDiffusion.py, the default interface V was a scalar #here is where the Stefan condition is expressed as a source term #dividing by solution variable because of ImplicitSourceTerm: see FiPy manual coeff = interfaceVelocVar * latentHeat * distanceVar.cellInterfaceAreas / (mesh.cellVolumes) / cp / temperature # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( #need to fix interface T at solidification T # http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-internal-boundary-conditions mask1=((distanceVar >= -cellSize/2.) & (distanceVar<cellSize/2.)) #TODO will this always select a single cell center? is syntax correct? ipython largeValue=1e+10 # (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((( #TODO or the == version? #TODO check signs transientCoeff=1. eq = TransientTerm(transientCoeff) - DiffusionTerm(D) - ImplicitSourceTerm(coeff) + ImplicitSourceTerm(largeValue*mask1) - largeValue*mask1*Tm return eq END OF MESSAGE
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
