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))

enter image description here

  • 6
    $\begingroup$ You are solving using a multi-step semi-explicit method. Since you use an explicit scheme for $\theta$ in the second equation, the stability of the overall scheme is limited by the parameters $\lambda$ and $\kappa$. If you solve the coupled equations in an implicit manner, you will not have such stability issues. $\endgroup$
    – Chenna K
    Sep 22 at 11:39
  • $\begingroup$ Thanks for your response. If I've understood you correctly the instability results from the fact that I'm taking θ from the previous iterative step in the equation for dθ/dt, correct? The trouble is that makes the final equation extremely complicated to solve. $\endgroup$
    – nickwinz
    Sep 26 at 16:33
  • $\begingroup$ Yes, your understanding is correct. It is not that complicated to solve. Give it a try. If you want unconditional stability, i.e., be able to use large time steps, then there is no alternative. $\endgroup$
    – Chenna K
    Sep 27 at 19:57

1 Answer 1


The oscillations may be due to the FE mass matrix. Try using a lumped (diagonal) mass matrix.

  • $\begingroup$ M is already diagonal. $\endgroup$
    – nickwinz
    Sep 26 at 16:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.