This is basically a follow-up of this question. I had tried to slightly simplify the actual problem and ask to get the basic idea/method and then implement it myself on my more complex PDAE system but am now running into some issues with the initial/boundary conditions so figured I will just ask the full, actual problem. Hopefully, it's not too large/involved.
I am trying to solve a PDAE system of 4 differential and 4 algebraic variables. I believe the formulation is OK since I have 4 initial conditions for the 4 variables, 4 Dirichlet BCs for the 4 variables, and 2 Neumann BCs for the 2 variables appearing as second-order space derivatives. The system is described in the pictures below (not sure why it was not formatting as math).
System [blue are constants, rest are f(t,x)]
First 4 PDEs of differential variables: $$ \frac{\partial C_{\mathrm{b}}}{\partial t}=-\frac{1}{\color{blue}{t_{\mathrm{f}} W}} \frac{\partial\left(C_{\mathrm{b}} F_{\mathrm{b}}\right)}{\partial x}+\color{blue}{D_{\mathrm{b}}} \frac{\partial^2 C_{\mathrm{b}}}{\partial x^2}-\frac{J_{\mathrm{w}} C_{\mathrm{p}}}{\color{blue}{t_{\mathrm{f}}}} $$ $$ \frac{\partial C_{\mathrm{p}}}{\partial t}=-\frac{1}{\color{blue}{t_{\mathrm{p}} W}} \frac{\partial\left(C_{\mathrm{p}}\left(F_{\mathrm{b} 0}-F_{\mathrm{b}}\right)\right)}{\partial x}+\color{blue}{D_{\mathrm{p}}} \frac{\partial^2 C_{\mathrm{p}}}{\partial x^2}+\frac{J_{\mathrm{w}} C_{\mathrm{p}}}{\color{blue}{t_{\mathrm{f}}}} $$ $$ \frac{\partial F_{\mathrm{b}}}{\partial t}=\left(-\color{blue}{W} J_{\mathrm{w}}-\frac{\partial F_{\mathrm{b}}}{\partial x}\right) \frac{F_{\mathrm{b}}}{\color{blue}{t_{\mathrm{f}} W}} $$ $$ \frac{\partial P_{\mathrm{b}}}{\partial t}=\left(-\color{blue}{b} F_{\mathrm{b}}-\frac{\partial P_{\mathrm{b}}}{\partial x}\right) \frac{F_{\mathrm{b}}}{\color{blue}{t_{\mathrm{f}} W}} $$
4 algebraic equations: $$ J_{\mathrm{w}}=\color{blue}{A_{\mathrm{w}}}\left(\left(P_{\mathrm{b}}-\color{blue}{P_{\mathrm{p}}}\right)-\color{blue}{R T_{\mathrm{b}}}\left(C_{\mathrm{w}}-C_{\mathrm{p}}\right)\right) $$ $$ J_{\mathrm{s}}=\color{blue}{B_{\mathrm{s}}} \exp \left(\frac{J_{\mathrm{w}}}{k}\right)\left(C_{\mathrm{b}}-C_{\mathrm{p}}\right) $$ $$ C_{\mathrm{w}}=C_{\mathrm{p}}+\exp \left(\frac{J_{\mathrm{w}}}{k}\right)\left(C_{\mathrm{b}}-C_{\mathrm{p}}\right) $$ $$ k=0.0953606 F_{\mathrm{b}}^{0.13} J_{\mathrm{w}}^{0.739} C_{\mathrm{b}}^{0.135} $$
Initial/Boundary conditions (for some reason wasn't rendering the last BC correctly in math mode):
The independent variables are $t\geq 0$ and $x \in [0,1]$. Here is what I tried.
Attempted Mathematica Solution
First, specify parameters:
tf = 0.0008; tp = 0.0005; L = 0.934; W = 8.4; R = \
0.082; Pp = 1; b = 8530; Aw = 9.5*10^-7; Bs = 8.5*10^-8; Db =
1.7*10^-9; Dp = 1.7*10^-9; Tb = 300; Tp = 300; Fb0 =
2.166*10^-4; Cb0 = 6.226*10^-3; Cp0 = 0; Pb0 = 5.83;
Then define equations and BCs:
eqns = {D[Cb[t, x], t] == -(D[Cb[t, x]*Fb[t, x], x]/(tf*W)) +
Db*D[Cb[t, x], x, x] - (Jw[t, x]*Cp[t, x])/tf +
NeumannValue[0, x == 1],
D[Cp[t, x], t] == -(D[Cp[t, x]*(Fb0 - Fb[t, x]), x]/(tp*W)) +
Dp*D[Cp[t, x], x, x] + (Jw[t, x]*Cp[t, x])/tf +
NeumannValue[0, x == 1],
D[Fb[t, x], t] == (-W*Jw[t, x] - D[Fb[t, x], x])*Fb[t, x]/(tf*W),
D[Pb[t, x], t] == (-b*Fb[t, x] - D[Pb[t, x], x])*Fb[t, x]/(tf*W),
Jw[t, x] == Aw*((Pb[t, x] - Pp) - R*Tb*(Cw[t, x] - Cp[t, x])),
Js[t, x] == Bs*E^(Jw[t, x]/kx[t, x])*(Cb[t, x] - Cp[t, x]),
Cw[t, x] == Cp[t, x] + E^(Jw[t, x]/kx[t, x])*(Cb[t, x] - Cp[t, x]),
kx[t, x] ==
0.095360572*Fb[t, x]^0.13*Jw[t, x]^0.739*
Cb[t, x]^0.135}; bcs = {Cb[t, 0] == Cb0, Cp[t, 0] == Cp0,
Fb[t, 0] == Fb0, Pb[t, 0] == Pb0, Cb[0, x] == Cb0, Cp[0, x] == Cp0,
Fb[0, x] == Fb0, Pb[0, x] == Pb0};
Finally, solve:
sol[t_, x_] =
NDSolveValue[{eqns, bcs}, {Cb[t, x], Cp[t, x], Fb[t, x], Pb[t, x],
Jw[t, x], Js[t, x], Cw[t, x], kx[t, x]}, {t, 0, 10}, {x, 0, 1}]
which gives this error:
NDSolveValue::femcnsd: The PDE coefficient -((W Cp[t,x] Jw[t,x]+Fb[t,x] (Cb^(0,1))[t,x]+Cb[t,x] (Fb^(0,1))[t,x]+tf W (Cb^(1,0))[t,x])/(tf W)) does not evaluate to a numeric scalar at the coordinate {5.,0.5}; it evaluated to -((0. +0.1 tf W)/(tf W)) instead.
I tried some tricks from various Q/A here, specifically forcing MethodOfLines and FiniteElement but then got problems with the BCs. So adding this to sol:
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement", "MinPoints" -> 41,
"MaxPoints" -> 41, "DifferenceOrder" -> 2}}}
gives
NDSolveValue::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable.
I also tried forcing $t$ to be the temporal variable by using "TemporalVariable" -> t but got
NDSolveValue::tvic: t cannot be used as the temporal independent variable because the conditions {Cb[0,x]==Cb0,Cp[0,x]==Cp0,Fb[0,x]==Fb0,Pb[0,x]==Pb0} for that dimension do not constitute sufficient initial conditions given at only one value of t.
NDSolveValue::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable.
PS. The 1st term in the 2nd PDE is technically $\frac{\partial\left(C_{\mathrm{p}} F_{\mathrm{p}}\right)}{\partial x}$ with $F_p=F_{b0}-F_b$ and $\frac{\partial F_{\mathrm{p}}}{\partial x} = - \frac{\partial F_{\mathrm{b}}}{\partial x}$ but I found the above to be less complex as it avoids a variable creation. I don't think it affects the actual solution.
