Hope everybody is safe during this COVID-19 pandemic!
I am having an issue numerically solving the following coupled diffusion equation with concentration dependent coefficients and need some suggestions on possible work arounds. The transport equation I wish to solve is:
$$ \theta^B\frac{\partial c_i^B}{\partial t} + \theta^D\frac{\partial c_i^D}{\partial t} =\theta^B\frac{\partial}{\partial x}\left(J_i^B\right) + \theta^D\frac{\partial}{\partial x}\left(J_i^D\right). $$
This an equation for diffusion of charged ions in a dual continuum system representative of charged porous media (e.g., membrane, clay). The superscript B and D refers to the properties at "bulk water" and "Donnan water" porosities, respectively. Considering a Donnan equilibrium, the concentrations and fluxes in these sub-domains are related as,
$$ c_i^D = c_i^Bexp\left(\frac{-z_iF\psi^D}{RT}\right) = c_i^D {BF}^{z_i} \\ J_i^B = -D_i^B\frac{\partial c_i^B}{\partial x} + \frac{z_iD_i^Bc_i^B}{\sum_j^N\left(z_j^2 D_j^B c_j^B\right)}\sum_j^N\left(z_jD_j^B\frac{\partial c_i^B}{\partial x}\right) \\ = \sum_j^N \left[\delta_{ijD_i^B} - \frac{z_iz_jD_i^BD_j^Bc_i^B}{\sum_j^N\left(z_j^2 D_j^B c_j^B\right)}\right]\frac{\partial c_i^B}{\partial x} \\ J_i^D = -D_i^D\frac{\partial c_i^D}{\partial x} + \frac{z_iD_i^Dc_i^D}{\sum_j^N\left(z_j^2 D_j^D c_j^D\right)}\sum_j^N\left(z_jD_j^D\frac{\partial c_i^D}{\partial x}\right) \\ = -D_i^D\frac{\partial ({BF}^{z_i} c_i^B)}{\partial x} + \frac{z_iD_i^D {BF}^{z_i} c_i^B}{\sum_j^N\left(z_j^2 D_j^D {BF}^{z_i} c_j^B\right)}\sum_j^N\left(z_jD_j^D\frac{\partial ({BF}^{z_i} c_i^B)}{\partial x}\right) \\ = \sum_j^N \left[\delta_{ijD_i^B}{BF}^{z_i} - \frac{z_iz_jD_i^BD_j^B{BF}^{z_i}{BF}^{z_j}c_i^B}{\sum_j^N\left(z_j^2 D_j^B {BF}^{z_i} c_j^B\right)}\right]\frac{\partial c_i^B}{\partial x} $$
Using these relationships, the above transport equation can be discretized with finite difference method, and using a backward Euler scheme, it can be represented in the following in the discretized form:
$$ \theta^B(\textbf{c}_i^{B,t+dt} -\textbf{c}_i^{B,t}) + \theta^D(\textbf{c}_i^{B,t+dt}\textbf{BF}^{z_i,t+dt} -\textbf{c}_i^{D,t}) = \textbf{M}^{t+dt}\textbf{c}_i^{B,t+dt} \\ (\theta^B + \theta^D \textbf{BF}^{z_i,t+dt} - \textbf{M}^{t+dt}) \textbf{c}_i^{B,t+dt} = (\theta^B\textbf{c}_i^{B,t} + \theta^D\textbf{c}_i^{D,t}) \\ $$
Where, $\textbf{M}$ contains all the cross-coupled flux entries for the finite difference method. Now this system of algebraic equations can be solved using the direct matrix solver. I have implemented these calculations in MATLAB and solve the above system of equations by "blackslash" function of matlab. However, the system of equations is highly nonlinear because both $\textbf{M}$ and $\textbf{BF}^{z_i}$ are dependent on the concentration $\textbf{c}$. To determine the so called "Boltzmann factor" (BF), two additional constraints, related to the electroneutrality of Donnan and bulk porosities, need to be considered:
$$ \text{in bulk water: }\sum_i^N z_ic_i^B = 0\\ \text{in Donnan water: }\sum_i^N z_ic_i^D + Q = 0 \\ \sum_i^N z_i{BF}^{z_i}c_i^B + Q = 0 $$ Here Q is the surface charge in the Donnan water and can be considered a constant in this case.
My problem is related to the effective solution of the above polynomial to determine BF. So far, I have implemented a Picard iterative loop to solve the system of equations. For a particular time step, the steps inside the loop is:
while loop
step 1: Calculate BF by the solving $\sum_i^N z_i{BF}^{z_i}c_i^B + Q = 0$
step 2: Calculate flux entries including diffusion parameters (entries of $\text{M}$)
- step 3: Form the matrix $\text{M}$
- step 4: Solve for concentration
end
Iteration continues until the concentration converges to a threshold. However, this approach does not and the iterative scheme is not able to update both $\textbf{c}$ and $\textbf{BF}$ ! The above steps seem to work without any error but the Boltmann factor stays the same and does not change as the concentration changes. So far, I have tried two matlab functions (fzero, and roots) to solve $\sum_i^N z_i{BF}^{z_i}c_i^B + Q = 0$, where:
fzero crashes as it sometimes produces the complex or undesired root (especially for the cubic and higher order polynomial cases)
roots works without crashing, and it is possible to manually screen the positive root (BF is always positive), but still it does not properly update BF in the iteration.
I wonder if this problem with BF is related to the fact that c and BF appear as multiplication form in the equation?
Please note that, for a system containing only one symmetrical electrolyte (charge, z = 1 and -1), $\sum_i^Nz_i{BF}^{z_i}c_i^B + Q = 0$ reduces to a quadrtic equation. In this condition, analytical solution for BF exists. If I plug in this solution in the iterative scheme, the calculation works and both concentration and Boltzmann factor is properly updated.
But for systems containing multivalent ions (z_i>=2), $ \sum_i^Nz_i{BF}^{z_i}c_i^B + Q = 0$ takes the form of cubic or higher order polynomials. In such cases, analytical solution does not exist and the numerical solution as used with the matlab functions does not work :(
Can anybody suggest any tips or tricks to resolve this issue? Sorry for the long post but I would really appreciate any suggestions as I am stuck!
Thank you