Daniel,
The code I was referring to is the code used in the paper titled " Diffusion
under temperature gradient: A phase-field model study " by Mohanty, Guyer and
Sohn. I found an archived conversation from this mailing list "
http://article.gmane.org/gmane.comp.python.fipy/3064/match=thermal+gradient "
where Jonathan Guyer mentioned digging up the FiPy code used for the paper and
was wondering if that would still be possible.
I have attached my code as well, in response to your other email. The code is
quite temperamental as I have combined elements from quite a few sources but it
should converge to a solution as it is written right now, albeit an unphysical
one.
Thank you for your time and interest,
Luke Johnson
----- Original Message -----
From: "Daniel Wheeler" <daniel.wheel...@gmail.com>
To: "Multiple recipients of list" <fipy@nist.gov>
Sent: Tuesday, April 29, 2014 2:52:29 PM
Subject: Re: Question about Diffusion Under Temperature Gradient
On Mon, Apr 28, 2014 at 7:00 PM, Johnson, Luke A
<lukejohn...@neo.tamu.edu> wrote:
Hello,
I am trying to create a relatively simple thermomigration model in FiPy for
a final project in my Kinetics class and have been using a paper titled
Diffusion Under Temperature Gradient as a reference. I can get my system to
converge and segregate if I assume constant values for the Chemical Mobility
and Mobility of Thermotransport but am having convergence problems when the
mobilities become functions of temperature and composition. I think my code
is right but the problem is that I am using unrealistic values for the
parameters not explicitly listed in the paper. Could you provide the
thermodynamic system parameters you used and maybe a workflow for the code
so I can check with my work?
Unfortunately, I am not sure what you're referring to. Which code /
example are you asking about?
--
Daniel Wheeler
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
### Summary: This script will create a model with which we can study thermal migration in single and double phase binary alloys. The thermal term is coupled to the equation via "deGroot". The most important factors are initial composition, atomic mobility and heat of transport.
### Context: Should be a fully contained script with few imported modules.
### phi = phase variable
### psi = transformation variable necessary to write coupled equation
### T = temperature variable
### DATE MODIFIED: 4/22/2014
#=====================================================================================================================
### STATE OF THE CODE: Grid has been drawn and the variables have been laid out over it. Temperature dependent thermo and chemical mobilities are coded, just need values to make it realistic. Boundary conditions for temperature and mass have just been introduced as well. Phase Field equations are coded and the variables are initialized as repeating zeros. Not sure if the 3rd equation is correct (because it changes the temperature profile, which should be an imposed part of the problem) but introducing it allowed the solver to at least try to solve the problem and plot some results.
### WORK ON NEXT:
### -input reasonable variable for comparison to Tunde's paper (added: Apr 21)
#=====================================================================================================================
from fipy import *
# boilerplate initialization function which also sets the number of cells
#=========================================================================
if __name__ == "__main__":
nx = ny = 50
else:
nx = ny = 50
# setting up the mesh and laying necassary variables over that mesh (T variable may not be necessary)
#=========================================================================
dx=.1
dy=.1
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)
phi = CellVariable(name=r"$\phi$", mesh=mesh) # concentration_b
psi = CellVariable(name=r"$\psi$", mesh=mesh) # intermediate variable for use in CahnHilliard equation implementation
T = CellVariable(name=r"$T$", mesh=mesh) # temperature
#Setting the temperature across the mesh to have a linear gradient in the x direction
T_hot=1500.
T_cold=1000.
T_gradient=(T_hot-T_cold)/(nx-1)
for height in range(0,ny):
index=height*nx
T[index]=T_cold
for dist in range(1,nx):
T[index+dist]=T[index+dist-1] + T_gradient
#T[:]=900
phi[:] = GaussianNoiseVariable(mesh=mesh, mean=0.5, variance=0.1).value
if __name__ == "__main__":
viewer = Viewer(vars=(phi)) # , datamin=0., datamax=1.)
# Setting up equations for chemical mobility and mobility of thermotransport respectively (THESE ONLY NEED TO BE SWITCHED ON IF WE'RE LOOKING AT TWO PHASE ALLOYS, SINGLE PHASE CAN BE APPROXIMATED AS CONSTANT ATOMIC MOBILITY AND HEAT OF TRANSPORT)
#=========================================================================
Vm = 10**(-5)
rho = 1/Vm
R=8.314
Beta=10**(-5)
l=10**(-3) #meters
#element a properties
Qa = -30 # heat of transport for element a
D_a = 10**(-12) # diffusivity of a
beta_ao = 5 # mobility constant for element b
beta_a = beta_ao*numerix.exp(-Qa/(R*T)) # atomic mobility of element a
C_a = 24 # specific heat of a
Qa_tilde = C_a + D_a*T # heat of transport for element a
#element b properties
Qb = 30 # heat of transport for element b
D_b = 10**(-11) # diffusivity of b
beta_bo = 5 # mobility constant for element b
beta_b = beta_bo*numerix.exp(-Qb/(R*T)) # atomic mobility of element b
C_b = 25 # specific heat of b
Qb_tilde = C_b + D_b*T # heat of transport for element b
Mc = 1#rho*phi*(1 - phi)*(phi*beta_a + (1 - phi)*beta_b) # temperature dependent chemical mobility (Tunde's paper)
Mq = 1#rho*phi*(1 - phi)*(beta_a*Qa_tilde - beta_b*Qb_tilde) # temperature dependent mobility of thermotransport (Tunde's paper)
# Setting up the free energy functional and its derivative
#=========================================================================
delta_f = R*(T_hot+T_cold)/2 # energy barrier between the two neighboring minima
kappac = 0.2 # gradient energy coefficient (Tunde's paper)
alpha = 600# some "positive constant" (Long-Qing Chen) (200 slow but separates... go to the 1000's to see it fast) ((affects the temp dependence of the energy well))
Tm = delta_f/R # median temperature
f = 4*delta_f*(-(phi**2)/2 + phi**4/4) + (15*alpha/8)*(phi - 2*phi**3/3 + phi**5/5)*(T - Tm) # free energy function (Long-Qing Chen paper)
f_deriv = (phi**2 - 1)*((phi**2 - 1)*(15*alpha/8*(T - Tm)) + phi*4*delta_f) # first derivative of free energy function
f_2deriv = 4*phi**3*(15*alpha/8*(T - Tm)) + 3*phi**2*4*delta_f - 4*phi*(15*alpha/8*(T - Tm)) - 4*delta_f # second derivative of free energy function (do this by hand and input)
#L_ab=9.6
#f = R*T*(phi*numerix.log(phi)+(1-phi)*numerix.log(1-phi))+L_ab*phi*(1-phi)
#f_deriv = R*T*(-numerix.log(1-phi)+numerix.log(phi))+L_ab*(-2*phi+1)
#f_2deriv = R*T/(phi-phi**2)-2*L_ab
#non-dimensionalizing parameters
Mq_= CellVariable(mesh=mesh, value=Mq*Vm/(delta_f*Beta))
Mc_= CellVariable(mesh=mesh, value=Vm*Mc/Beta)
f_=Vm*f/delta_f
f_deriv_=f_deriv*Vm/delta_f
f_2deriv_=f_2deriv*Vm/delta_f
kappac_=kappac*Vm/(delta_f*l**2)
# writting the actual thermomigration governing equation (split into two because of the 4th order aspect of the CahnHilliard situation
# fourth order because we've got coupled thermal and concentration gradients so twice as many gradients and they're nested inside one another
#=========================================================================
eq1 = TransientTerm(var=phi) == DiffusionTerm(coeff=Mc_,var=psi) - (Mq_.faceValue * (T.faceGrad)/T.faceValue).divergence
eq2 = ImplicitSourceTerm(coeff=1,var=psi) == ImplicitSourceTerm(coeff=-f_2deriv_, var=phi) - f_2deriv_*phi + f_deriv_ - DiffusionTerm(coeff=kappac_, var=phi)
eq = eq1 & eq2
# Boundary and Initial Conditions (MIGHT NOT BE NECESSARY >>>> BOUNDARY CONDITIONS HAVE NO EFFECT ON ""SOLUTION"" THAT FiPy IS OUTPUTTING (Apr 22)
#=========================================================================
# concentration BCs
BC_left = FixedFlux(faces=mesh.facesLeft, value=0)
BC_right = FixedFlux(faces=mesh.facesRight, value=0)
BC_top = FixedFlux(faces=mesh.facesTop, value=0) # specify no flux through the top of the domain
BC_bottom = FixedFlux(faces=mesh.facesBottom, value=0) # specify no flux through the bottom of the domain
# temperature BCs
#BC_Tleft = FixedFlux(faces=mesh.facesLeft, value=0)
#BC_Tright = FixedFlux(faces=mesh.facesRight, value=0)
#BC_Ttop = FixedFlux(faces=mesh.facesTop, value=0) # specify no flux through the top of the domain
#BC_Tbottom = FixedFlux(faces=mesh.facesBottom, value=0) # specify no flux through the bottom of the domain
BCs=(BC_left,BC_right,BC_top,BC_bottom)
# MAY NEED TO INCLUDE HIGHER ORDER BOUNDARY CONDITIONS IN THIS SINCE WE HAVE A PHI^5 TERM IN THE FREE ENERGY FUNCTIONAL
#steps = 50
#for step in range(steps):
#eq.solve(dt=0.01)
#viewer.plot()
dexp = -5
elapsed = 0.
if __name__ == "__main__":
duration = 100
else:
duration = 100
while elapsed < duration:
dt = min(100, numerix.exp(dexp))
elapsed += dt
dexp += 0.03
eq.solve(boundaryConditions=BCs,dt=dt)
if __name__ == "__main__":
viewer.plot()_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
