Message Boards Message Boards

Make a 2D Electro-diffusion Model with NDSolve?

GROUPS:

Set Up

I am attempting to make a 2D model of the diffusion of ions in solution due to concentration gradients and an electric field.

The equation that needs to be solved to describe the diffusion is

∂C/∂t=-∇∙(-D∇C-zμC∇ϕ)

C: Concentration of species in a given space

D:Diffusion Coefficient

z:Charge of a species molecule

μ:Mobility of a species molecule

ϕ:Electric potential at a given point

-

The code is as follows:

All constants are set equal to "1"

diffusionC = 1;
charge = 1;
mobility = 1;

The equation that describes the electric field is defined

\[CapitalPhi][t_, x_, y_] := -2.5*x + 5;

The gradients inside the parenthesis are defined

diffGrad = Grad[u[t, x, y], {x, y}];
elecGrad = Grad[\[CapitalPhi][t, x, y], {x, y}];

The terms inside the parenthesis are defined

equationP1 = -diffusionC*(diffGrad[[1]] + diffGrad[[2]]);
equationP2 = charge*mobility*u[t, x, y]*elecGrad[[1]];

The final gradient outside the parenthesis is defined

fullGrad = Grad[equationP1 + equationP2, {x, y}]

The full equation is defined with flux boundary conditions for the top, bottom, and right sides

equation = 
   D[u[t, x, y], t] + fullGrad[[1]] == 
   NeumannValue[0, y == 0 || y == 2] + NeumannValue[3, x == 2]

The area over which the equation is to be solved is defined

R = ImplicitRegion[0 <= x <= 2 && 0 <= y <= 2, {x, y}];

The initial condition and boundary condition for the left side are defined

ic = u[0, x, y] == 1
bc = u[t, 0, y] == 8 - 7*Exp[-1000*t]

The equation is solved

sol = NDSolve[{equation, ic, bc}, {u, \[CapitalPhi]}, {t, 0, 
    10}, {x, y} \[Element] R]

The solution is animated

ListAnimate@
Table[Plot3D[u[t, x, y] /. sol, {x, 0, 2}, {y, 0, 2}, 
PlotRange -> All, AxesLabel -> {x, y, u}], {t, 0, 10, .025}]


ListAnimate@
Table[Plot3D[\[CapitalPhi][t, x, y] /. sol, {x, 0, 2}, {y, 0, 2}, 
PlotRange -> All, AxesLabel -> {x, y, u}], {t, 0, 10, .025}]

Problem

Upon running the line to solve the equation, the following errors are presented.

NDSolve::femcscd: The PDE is convection dominated and the result may not be 
stable. Adding artificial diffusion may help. >>

NDSolve::ndsz: At t == 0.06490840231453537`, step size is effectively zero; 
singularity or stiff system suspected. >>

The solution given is:

enter image description here

which only solves up until the time given in the error rather than to "10" as defined in the code.

I request a solution to resolve the errors and get the model to solve properly for the desired time range. If additional information is required, please do let me know. Thank you for your time and consideration.

POSTED BY: Kali Ellison
Answer
2 months ago

Group Abstract Group Abstract