From: Torsten Hennig on
> Hello,
> I am having trouble with diffusion with drift and the
> conservation of mass in pdepe. The model I am trying
> to simulate is described as follows:
>
> Within a container is a fixed number of particles in
> solution that is subject to diffusion and an external
> force. The external force is analogous to a uniform
> electric field that influences each particle in the
> same direction and magnitude no matter its position
> within the container. The variable of interest will
> be the concentration of particles at a given position
> and time. I want my simulation to show the steady
> state behavior of such a system. The conservation of
> mass should hold as no matter will enter or leave the
> container. I would expect that the mass will collect
> along the side of the container perpendicular to the
> applied force and form a gradient in the opposite
> direction of the force, due to the equilibrium
> between diffusion and the flow created by the force.
>
> To model this as simply as possible the container is
> reduced to a single dimension, an interval of fixed
> length [0,L].
>
> The concentration, u, at any given point and time
> will satisfy this equation,
> u_t = -D*u_xx + V*u_x
> where u_t is the partial derivative w.r.t. time and
> u_xx is the second partial derivative w.r.t. its
> position x, etc.
>
> There can be any arbitrary initial condition u0
>
> and boundary conditions... This is where I begin to
> get confused...
> I am uncertain about what boundary conditions I
> should use and how to enter these into MATLAB's
> pdepe. I believe I am supposed to arrange some kind
> of reflective boundary conditions, but I am uncertain
> how to apply this in practice. I do not want to fix
> the concentrations at the edges, because the
> concentrations should be allowed to change there. I
> first thought I should arrange a no flux condition
> and set the Neumann conditions to zero, but this is
> not the answer, my solutions never reach steady
> state. However, when I remove the "drift" term, i.e.
> set V=0, diffusion reaches a steady state and
> conservation of mass appears to hold.
>
> I then learned about reflective boundary conditions,
> but do not know how to implement them or if this is
> even correct.
>
> If anyone has an answer, could you please give me an
> explanation and the how to implement these boundary
> conditions in the pde format of "pl=?, ql=?, pr=?,
> qr=?"?
> Thanks,
> KWK

Your equation reads dc/dt = -v*dc/dx + D*d^2x/dx^2
or dc/dt = d/dx(-v*c+D*dc/dx).
Integrating both sides from x=0 to x=L gives
d/dt (int_{x=0}^{x=L} c dx) =
(-v*c+D*dc/dx)|x=L - (-v*c+D*dc/dx)|x=0.

So if you want no particles leave on both ends of your
container, the boundary condition is
-v*c+D*dc/dx = 0 at x=0 and x=L.

This is easily achieved by setting
f=-v*c+D*dc/dx and letting ql=qr=1 and pl=pr=0
in the boundary routine for pdepe.

Best wishes
Torsten.