I'm looking for a bit of assistance in how to organize a FiPy script I am
trying to write. I am trying to model transient 1D diffusion into and out of
soil with an exponentially decaying boundary condition. I assuming a starting
concentration and that diffusion occurs for an initial using a low first-order
rate constant (-0.05/yr). After this period, an increased rate is applied
(-1.0/yr) for three years and then the rate decreases back to the initial value
(-0.05/yr).
The script I have is below. I'm not sure the way I am setting my left boundary
condition is correct. Any advice we be greatly appreciated.
from fipy import *
nx = 250
dx = 0.01
mesh = Grid1D(nx=nx, dx=dx)
# create a CellVariable and initialize it to zero
phi = CellVariable(name="Concentration", mesh=mesh, value=0.)
# create a diffusion equation where Deff = D0*tau*por/R
D0 = 0.0315576 # m**2/yr or 10**-9 m**2/sec
tau = 0.25
por = 0.35
R = 2.0
Deff = D0*tau*por/R
time = Variable()
D = Deff
eqI = TransientTerm() == DiffusionTerm(coeff=D)
## INITIAL PERIOD FROM 1970 TO 2013 ##
# apply boundary conditions
# variable concentration boundary assuming first-order attenuation
# and that release occurred in 1970;
valueLeft = 250. * numerix.exp(-0.05*time)
valueRight = 0
phi.constrain(valueLeft, mesh.facesLeft)
phi.constrain(valueRight, mesh.facesRight)
# create a viewer to see the results
if __name__ == '__main__':
viewer = Viewer(vars=phi, title="TCE Diffusion Low Permeability Sediment",
datamin=0., datamax=250.)
viewer.plot()
# assume 43 yrs (1970 to 2013) and time increment of 0.1
steps = 430
dt = 0.1
while time() <= 43:
time.setValue(time() + dt)
eqI.solve(var=phi, dt=dt)
if __name__ == '__main__':
viewer.plot()
if __name__ == '__main__':
raw_input("Time-dependent boundary condition. Press <return> to proceed...")
TSVViewer(vars=(phi, phi.grad)).plot(filename="TCE_2013.tsv")
## TREATMENT PERIOD FROM 2013 TO 2016 ##
# Change in boundary condition of left face
del phi.faceConstraints
valueLeft = phi[0] * numerix.exp(-1.0*(time-43))
phi.constrain(valueLeft, mesh.facesLeft)
phi.constrain(0., mesh.facesRight)
eqI = TransientTerm() == DiffusionTerm(coeff=D)
# assume additional 3 yrs (2013 to 2016) and time increment of 0.1
steps = 30
dt = 0.1
while time() <= 46:
time.setValue(time() + dt)
eqI.solve(var=phi, dt=dt)
if __name__ == '__main__':
viewer.plot()
if __name__ == '__main__':
raw_input("Time-dependent boundary condition. Press <return> to proceed...")
TSVViewer(vars=(phi, phi.grad)).plot(filename="TCE_2016.tsv")
## POST-TREATMENT PERIOD FROM 2016 TO 2020 ##
# Change in boundary condition of left face
del phi.faceConstraints
valueLeft = phi[0] * numerix.exp(-0.05*(time-46))
phi.constrain(valueLeft, mesh.facesLeft)
phi.constrain(0., mesh.facesRight)
eqI = TransientTerm() == DiffusionTerm(coeff=D)
# assume additional 4 yrs (2016 to 2020) and time increment of 0.1
steps = 40
dt = 0.1
while time() <= 50:
time.setValue(time() + dt)
eqI.solve(var=phi, dt=dt)
if __name__ == '__main__':
viewer.plot()
if __name__ == '__main__':
raw_input("Time-dependent boundary condition. Press <return> to proceed...")
TSVViewer(vars=(phi, phi.grad)).plot(filename="TCE_2020.tsv")
## POST-POP PERIOD FROM 2020 TO 2030 ##
# assume additional 10 yrs (2020 to 2030) and time increment of 0.1
f = open('Surface.txt', 'a')
steps = 100
dt = 0.1
while time() <= 60:
time.setValue(time() + dt)
eqI.solve(var=phi, dt=dt)
if __name__ == '__main__':
viewer.plot()
f.close()
if __name__ == '__main__':
raw_input("Time-dependent boundary condition. Press <return> to proceed...")
TSVViewer(vars=(phi, phi.grad)).plot(filename="TCE_2030.tsv")
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
