I used NDSolve to solve the 1D Laplace equation for phi with a Dirichlet boundary condition on the left boundary and a Neumann boundary condition on the right. The calculation for phi relies on the results of the calculation for mu; however, mu simply has Dirichlet boundary conditions on both the left and right. The working code for which is below:
Needs["NDSolve`FEM`"]
e = 1.60217662*10^-19;
sige = 3.37*10^-4;
sigi = 18;
F = 96485;
n = -0.02;
c = 1;
pO2 = 1.52*^-19;
Ie = -(2*F)*(c*pO2^n );
mu2 = -5.98*^-19;
l = 10*10^-6;
bcmu = {DirichletCondition[mu[x] == 0, (x == 0)],
DirichletCondition[mu[x] == mu2, (x == l)]};
solmu = NDSolve[{0 == D[D[mu[x], x], x], bcmu}, mu, {x, 0, l}];
Plot[mu[x] /. solmu, {x, 0, l}, PlotLabel -> "MuO2 vs. X-Position",
AxesLabel -> {"Position", "Mu"}]
bc1phi = DirichletCondition[phi[x] == 0, (x == 0)];
A = (Ie - sigi/(4*e)*(D[mu[x] /. solmu, x] /. x -> l))/(-sigi);
solphi = NDSolve[{D[D[phi[x], x], x] == 0 + NeumannValue[-A, x == l],
bc1phi}, phi, {x, 0, l}];
Plot[phi[x] /. solphi, {x, 0, l}, PlotLabel -> "Phi vs. X-Position",
AxesLabel -> {"Position", "Phi"}]
I would like to increase the complexity by moving to a 2D Laplace equation with the same boundary conditions on the left and right for both mu and phi but now zero flux boundary conditions on the top and bottom as well. To do so, I modified the code for 1D as such:
Needs["NDSolve`FEM`"]
e = 1.60217662*10^-19;
sige = 3.37*10^-4;
sigi = 18;
F = 96485;
n = -0.02;
c = 1;
pO2 = 1.52*^-19;
Ie = -(2*F)*(c*pO2^n );
mu2 = -5.98*^-19;
l = 10*10^-6;
h = 0.01;
meshRefine[vertices_, area_] := area > 10^-10;
mesh = ToElementMesh[
DiscretizeRegion[ImplicitRegion[True, {{x, 0, l}, {y, 0, h}}]],
MeshRefinementFunction -> meshRefine];
bcmu = {DirichletCondition[mu[x, y] == 0, (x == 0 && 0 < y < h)],
DirichletCondition[mu[x, y] == mu2, (x == l && 0 < y < h)]};
solmu = NDSolve[{Laplacian[mu[x, y], {x, y}] ==
0 + NeumannValue[0, y == 0 || y == h], bcmu},
mu, {x, y} \[Element] mesh];
bcphi = DirichletCondition[phi[x, y] == 0, (x == 0 && 0 < y < h)];
A = (Ie - sigi/(4*e)*(D[mu[x] /. solmu, x] /. x -> l))/(-sigi);
solphi = NDSolve[{Laplacian[phi[x, y], {x, y}] ==
0 + NeumannValue[0, y == 0 || y == h] +
NeumannValue[-A, x == l && 0 < y < h], bcphi},
phi, {x, y} \[Element] mesh];
Plot3D[phi[x, y] /. solphi, {x, 0, l}, {y, 0, h},
PlotLabel -> "Phi vs. X and Y Position",
AxesLabel -> {"X-Position", "Y-Position", "Phi"}]
From this code, I get an error "CompiledFunction: Could not complete external evaluation; proceeding with uncompiled evaluation." and I have to abort the evaluation or Mathematica stops responding. I believe the error has to do with the calculation of A, the value of the flux at the right boundary for phi, but I'm not sure how to fix this. In this case, the flux for phi should still be a constant independent of y because dmu/dx should be independent of y, but I would like to define the boundary as if dmu/dx at x=l could be a function of y. Does anyone know how to properly define A and the Neumann boundary condition so this works?
