# Solve 1D saturation equation using FEM

I want to solve the following 1D equation to model the time evolution of saturation in a porous medium using finite elements (I'm a beginner). The permeability $$K$$ and the pressure $$p$$ depend on the saturation $$s$$. The initial condition in the field $$\Omega = (0,l)$$ is given as: $$s(t=0)=0.1, \text{ in } \Omega$$ $$s[0] = 0.5; \text{left boundary}$$

$$\frac{\partial s}{\partial t} = - \frac{1}{\phi} \nabla \cdot u \\ u = \frac{K(s)}{\mu} \nabla p(s)$$

For the time-development, I discretized the time derivative using finite differences (explicit Euler). I'm aware that this is not optimal, because a linear system needs to be solved to obtain the new timestep, but I think that it should work for a preliminary solution.

I used linear shape functions and the standard Galerkin approach (shape function = weight function). My approximate solution is: $$s = \sum_i \tilde s_i N_i$$ Where $$N$$ is the linear shape function. For the discrete weak formulation, with the exponent $$n$$ denoting the time, I get:

$$MS^{n+1} = MS^{n} + \frac{\Delta t}{\phi}\frac{K(S^n)}{\mu}Dp^n$$ With the mass matrix: $$m_{ij} = \int_\Omega N_iN_jd\Omega$$ and the stiffness matrix: $$d_{ij} = \int_\Omega \nabla N_i \nabla N_j d\Omega$$.

However, the approach seems to be faulty because the solution doesn't converge and it oscillates strongly at the saturation jump. I had a look at different problems like the heat equation but found it to be different from my system. Especially because the saturation dependencies of $$K$$ and $$p$$ are strongly non-linear (exponential functions).

Could it be that the Standard Galerkin approach is not applicable? I would be very grateful for any hints.

PS: I could successfully simulate the case using a particle-based simulation (SPH).

• What a batty problem :p. Just so you know, explicit euler doesn't require solution of a matrix. Implicit euler does. Which one are you using? – EMP May 22 at 14:41
• Haha thanks :D I know that explicit Euler doesn't require solving a matrix system and I'm not sure how to call my procedure. I solved my rhs with the values of the current saturation, so I think it's explicit. Actually I saw a similar problem occurring for the heat equation and there is a way around this called mass lumping, but it seems to be quite elaborate – murcielagos May 22 at 15:11
• @murcielagos Explicit-in-time schemes like what you're using usually require very small timesteps and the errors you're describing sound like that might be the problem. Using implicit Euler is usually the remedy but since your problem is nonlinear that makes life harder. If it were me I'd back off to writing a solver for the linearized problem (the heat equation) first to work out the kinks, then introduce one nonlinearity (either $p$ or $K$), and finally use both of them. – Daniel Shapero May 22 at 16:18
• Thanks @DanielShapero! I'm not sure about the explicit time step, because it worked for the SPH solver and I also tried tiny steps. I will try the implicit scheme though and trace the error by first simplifying the problem :) – murcielagos May 22 at 18:08
• What are the expressions for $p$ and $K$? What are the values of $\psi$ and $\mu$? – knl May 28 at 12:51