# Spurious oscillations in solving diffusion problems using finite elements

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?

• 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. Sep 22, 2022 at 11:39
• The oscillations may be due to the FE mass matrix. Try using a lumped (diagonal) mass matrix. Sep 24, 2022 at 16:50
• M is already diagonal. Sep 26, 2022 at 16:11
• 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. Sep 26, 2022 at 16:33
• @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. Oct 9, 2022 at 13:47