JDS19
commented
3 years ago Correction: The flux term is $ N_i = -D\nabla c_i - DF/(RT)z_ic_i\nabla \phi_2 $ and thus the replaced equations in the code are:
# Constants in equations
beta = ct*F/R/T
mig = D*F/R/T
bv = a*i0/F/por
# Equation for first species
eq1 = TransientTerm(var=c1) == DiffusionTerm(coeff=D,var=c1) \
# add migration
+ DiffusionTerm(coeff=c1*z1*mig,var=phi2)
# subtract Faradaic term
- bv*c1/c1_0*numerix.exp(z1*bv*(phi_1-phi2-phi_0)) \
- ImplicitSourceTerm(coeff=bv*(-c2/c2_0),var=None)
eq2 = TransientTerm(var=c2) == DiffusionTerm(coeff=D,var=c2) \
+ ImplicitSourceTerm(coeff=-bv/c2_0,var=c2)
# z2=0 for ^^
+ ImplicitSourceTerm(bv*c1/c1_0*numerix.exp(z1*beta*(phi_1- \
phi2 - phi_0)),var=None)
# Equation relating flux to current for electrolyte potential
# Can't simplify to one Diffusion term because
# of z1 multiplying both flux terms
eq3 = DiffusionTerm(coeff=c1*mig,var=phi2) + mig*c1.grad.dot(phi2.grad) \
# Add diffusion
# + D*z1*c1.faceGrad.divergence <--- don/t use this?
+ DiffusionTerm(coeff=z1*D*c1,var=None)
# Add Faradaic contributions,
# can't use ImplicitSourceTerm because it's nonlinear in phi
+ bv*(c1/c1_0*numerix.exp(z1*beta*(phi_1-phi2-phi_0)) \
-c2/c2_0)The question stays the same on if I'm including these nonlinear terms in each equation correctly.
guyer
I wrote this problem on another site but realized this site is probably the better place for it. I'm solving a nonlinearly coupled system of PDEs. Here's the post, with equation format used from Jupyter Notebook formatting (I'm new to this site, so if it's easier for me to add images of the equations, I can do that): The system I am solving is described as: $ \frac{\partial c_1}{\partial t} = -\nabla N_1 - \nabla i/F $ \ $ \frac{\partial c_2}{\partial t} = -\nabla N_2 + \nabla i/F $ \ $ \sum{z_i \nabla N_i} = \nabla i/F $ \ where $ N_i = -D\nabla c_i - Dz_ic_i\nabla \phi_2 $ \ and $ \nabla i = ai0/por(c_1/c_1^0 e^{z_1ctF/(RT)(\phi_1 - \phi_2 - \phi_0)} - c_2/c_2^0 e^{z_2ctF/(RT)(\phi_1 - \phi_2 - \phi_0)}) $
a,i0,F,D,phi_1,z_1,ct,R,T, c_1^0, and c_2^0 are constants phi_2, c_1, and c_2 are time-dependent, time-evolving variables to solve for. In the problem, z_2 is 0 which simplifies the flux and exponent terms for c_2.
In implementing the flux terms in my code given by N, the chain rule is used to better couple the system and try to stay implicit:
$\nabla (c \nabla \phi) \equiv \nabla c \cdot \nabla \phi + c \nabla^2 \phi$
For simplicity, start with 0 flux boundary conditions on the bounds of x on [0,L]. The initial conditions are:
$ c_1(x,0) = 0.1 $ \ $ c_2(x,0) = 0. $ \ $ \phi_2 (x,0) = \phi_0 $
When I run the attached code, I get "RuntimeError: Factor is exactly singular" I've messed around with this code a bit, and in some formulations of this problem, with very small time changes, the code will solve, but the values don't change. I've looked into changing the iterations or tolerance of the solver, but it either doesn't change the final solution or the code doesn't run. Another issue is the long runtime because dt seems limited by the diffusion coefficient. Here's the code. Any help is greatly appreciated, as I've tried making my own solver but found that FiPy can greatly reduce my number of lines of code so it may be easier to de-bug.
Python is version 3.8.8 used in Jupyter Notebook. FiPy is version 3.4.2. Thank you again for any input.