examples.convection.exponential1D.mesh1D¶
Solve the steady-state convection-diffusion equation in one dimension.
This example solves the steady-state convection-diffusion equation given by
with coefficients \(D = 1\) and \(\vec{u} = 10\hat{\imath}\), or
>>> diffCoeff = 1.
>>> convCoeff = (10.,)
We define a 1D mesh
>>> from fipy import CellVariable, Grid1D, DiffusionTerm, ExponentialConvectionTerm, Viewer
>>> from fipy.tools import numerix
>>> L = 10.
>>> nx = 10
>>> mesh = Grid1D(dx=L / nx, nx=nx)
>>> valueLeft = 0.
>>> valueRight = 1.
The solution variable is initialized to valueLeft:
>>> var = CellVariable(mesh=mesh, name="variable")
and impose the boundary conditions
with
>>> var.constrain(valueLeft, mesh.facesLeft)
>>> var.constrain(valueRight, mesh.facesRight)
The equation is created with the DiffusionTerm and
ExponentialConvectionTerm. The scheme used by the convection term
needs to calculate a Péclet number and thus the diffusion term
instance must be passed to the convection term.
>>> eq = (DiffusionTerm(coeff=diffCoeff)
... + ExponentialConvectionTerm(coeff=convCoeff))
More details of the benefits and drawbacks of each type of convection
term can be found in Numerical Schemes.
Essentially, the ExponentialConvectionTerm and PowerLawConvectionTerm will
both handle most types of convection-diffusion cases, with the
PowerLawConvectionTerm being more efficient.
We solve the equation
>>> eq.solve(var=var)
and test the solution against the analytical result
or
>>> axis = 0
>>> x = mesh.cellCenters[axis]
>>> CC = 1. - numerix.exp(-convCoeff[axis] * x / diffCoeff)
>>> DD = 1. - numerix.exp(-convCoeff[axis] * L / diffCoeff)
>>> analyticalArray = CC / DD
>>> print(var.allclose(analyticalArray))
1
If the problem is run interactively, we can view the result:
>>> if __name__ == '__main__':
... viewer = Viewer(vars=var)
... viewer.plot()
FiPy