CODEWITHSUNDEEP
differential-equationsfipy
jso
jso

Reputation: 23

Is it possible to add source terms on outer boundaries only in FiPy?

Is it possible to solve a problem on the form

steady state conservation equation

where the source term is zero everywhere on the mesh except for at the boundaries, and how can this be represented in FiPy?

I have an equation that I have so far solved in 1D with the source term set to zero everywhere. However now I would like to add a sink term (sink) to the left hand side of the mesh (x=0). However, I can't get this to work exactly as I want.

As the sink is exactly at the boundary it would make sense to constrain the value at the face to be sink. However this still results in the all-zero solution, and I assume face values are not really supposed to be used this way. If I instead set the value in the first cell center to be sink, then I get a solution on the form I would expect. However now the sink is not actually at the boundary.

Is it possible to get the source/sink exactly at the boundary?

Here is a code sample, where the flux N is defined by a simple Fickian diffusion equation:

  L = 20                            # thickness
  nx = 100                          # number of mesh points
  dx = L/nx                         # distance between mesh points
  mesh = Grid1D(nx = nx, dx = dx)   # 1D grid of length L

  sink = -10

  c = CellVariable(mesh = mesh, name = r'$c_i$')
  c.constrain(100, mesh.facesRight)

  diffusionCoefficient = CellVariable(mesh = mesh, value = 5, name = r'$D_i$')

  sourceTerm = CellVariable(mesh = mesh, name = r'$S_i$')
  #sourceTerm.constrain(sink, mesh.facesLeft)
  sourceTerm[0] = sink

  eq = DiffusionTerm(coeff = -diffusionCoefficient, var = c) == sourceTerm

  eq.solve()

  viewer = Matplotlib1DViewer(vars = c)
  viewer.plot()
  input('Press return to continue')

Edit: mathematical expression of the problem and boundary conditions

The PDE describing diffusion of species i is

main PDE

and the source term is defined as

source term

describing a phase change from phase i to j on the boundary where cj is the concentration of the same species in the second phase (e.g. liquid/vapor). K is a rate constant. There is assumed to be no phase change other than on the boundary where the two phases meet. In my example code I just put a constant value as I was mostly wondering how to get the source term only on the boundary, I understand how to represent this source if it was on the entire mesh I think. The concentration of the second phase is calculated separately, on another mesh, with a similar equation. This should also have the same kind of source term, with opposite sign, on its right boundary, to ensure mass conservation.

Upvotes: 1

Views: 88

Answers (0)

Related Questions