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, privacy policy and cookie policy

Browse other questions tagged or ask your own question.