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

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?

  • 7
    $\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, 2022 at 11:39
  • $\begingroup$ The oscillations may be due to the FE mass matrix. Try using a lumped (diagonal) mass matrix. $\endgroup$
    – David T
    Sep 24, 2022 at 16:50
  • $\begingroup$ M is already diagonal. $\endgroup$
    – nickwinz
    Sep 26, 2022 at 16:11
  • $\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, 2022 at 16:33
  • 1
    $\begingroup$ @nickwinz, you have to solve the equations together since the two equations in $C$ and $\theta$ are coupled. That is, you will have a $2 \times 2$ block matrix system after discretisation. $\endgroup$
    – Chenna K
    Oct 9, 2022 at 13:47


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.