# 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 / dt)) • 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 at 11:39
• 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 at 16:33
• 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. Sep 27 at 19:57