I've been struggling with this problem for a while so I hope someone can help me here.
I'm trying to solve the McNabb-Foster equations for hydrogen diffusion in metals using a simple 1D finite element simulation in Python. This works fine under isothermal conditions, but I can't make it work for anisothermic conditions (i.e. increasing temperature). The equations are essentially a non-linear Poisson equation: $$ \frac{dC}{dt}=D ∇^2 C-N\frac{dθ}{dt} $$ where $$ \frac{dθ}{dt}=κC(1-θ)-λθ $$
Where C is the hydrogen concentration, D, κ and λ are temperature-dependent terms, N is constant and θ is the so-called trap occupancy (traps are sites with high binding energy). Using the backward difference method the first equation can be expressed in matrix form as: $$ (\frac{1}{∆t}[M]+D[K])\{C\}_{t+∆t}={F}_{t+∆t}+\frac{1}{∆t}[M]\{C\}_t+[M]N\Bigl\{\frac{dθ}{dt}\Bigl\}_{t+∆t} $$ Where [M] and [K] are the system matrices. So far, so good, but the tricky part is how to incorporate θ. For the isothermal case I have done this again using the backward difference method, but taking θ from the previous iterative step: $$ \Bigl\{\frac{dθ}{dt}\Bigl\}_{t+∆t}=\frac{\{θ\}_{t+∆t}-\{θ\}_t}{Δt}=κ\{C\}_{t+∆t}\{1-θ\}_t-λ\{θ\}_t $$ Combining the above two equations gives: $$ (\frac{1}{∆t}[M]+D[K]+[M]Nκ\{1-θ\}_t )\{C\}_{t+∆t}=\{F\}_{t+∆t}+\frac{1}{∆t}[M]\{C\}_t+[M]Nλ\{θ\}_t $$ Below is the iterative part of my code where I have tried to implement this. Basically, I'm solving the equations twice; once to find the new concentrations C and a second time to find the flux F, which is actually what I need.
The trouble is that I get (spurious?) oscillations at the boundary for reasonable input parameters (see image). These disappear for very small time increments, but it makes the code too slow (I'm running it in a least-squares routine so it needs to be fast). The mesh size has no influence.
Does anyone have any suggestions how I can improve the stability?
for i in range(1, Iter+1):
Clold = Cl # save Cl from last iteration
T = Tstart + i * dt * dTdt # calculate temperature
# calculate temperature dependent terms
Dl = D0 * np.exp(-El / R / T)
kappa = kappa0 * np.exp(-Et / R / T)
lambd = lambda0 * np.exp(-Ed / R / T)
# establish reverse difference equation
LHS = sysm + dt * Dl * sysk + dt * Nt * kappa * sysm * (1-theta)
RHS = np.dot(sysm, Clold) + dt * np.dot(sysm, (Nt * lambd * theta))
# apply constraints & boundary conditions
LHS[0, :] = 0.0
LHS[0, 0] = 1.0
RHS[0, 0] = 0.0
# calculate new Cl
Cl = np.linalg.solve(LHS, RHS)
# calculate flux matrix
LHS = sysm + dt * Dl * sysk + dt * Nt * kappa * sysm * (1-theta)
RHS = np.dot(sysm, Clold) + dt * np.dot(sysm, (Nt * lambd * theta))
F = np.dot(LHS, Cl) - RHS
# calculate new trap occupancies
theta = theta + dt * (kappa * Cl * (1 - theta) - lambd * theta)
# store flux for plotting
Fstore = np.append(Fstore, np.array(-F[0] / dt))
EDIT: using an implicit scheme the last equation would look like this:
$$ (\frac{1}{∆t}[M]+D[K]+[M]Nκ\{1-θ\}_{t+∆t} )\{C\}_{t+∆t}=\{F\}_{t+∆t}+\frac{1}{∆t}[M]\{C\}_t+[M]Nλ\{θ\}_{t+∆t} $$
So now you have a second unknown term in the equation. The new θ can be given in terms of the new C, i.e.:
$$ θ_{t+∆t}=\frac{∆tκC_{t+∆t}+θ_t}{1+∆tκC_{t+∆t}+∆tλ} $$
but this makes the first equation non-linear. Any suggestion how to solve this?