I am playing with Galerkin Least-Squares stabilization to solve advection diffusion problem in the context of the finite element method. This works very well for steady-state advection-diffusion equation, even in the High Péclet limit, when the problems are solved in the steady state.

The equation we solve is : $$\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = D \nabla^2 \alpha $$

with $ \mathbf{u}$ a velocity vector field which is known and $\alpha$ a scalar. $D$ is a diffusion coefficient. The weak form (without boundary terms) is : $$ \int \phi_i \frac{\partial \alpha}{\partial t} d\Omega + \int \phi_i \mathbf{u} \nabla \alpha d\Omega + \int D \nabla\phi_i \cdot \nabla\alpha d\Omega = 0$$

Adding GLS stabilization leads to the following: $$ \int \phi_i \frac{\partial \alpha}{\partial t} d\Omega + \int \phi_i \mathbf{u} \nabla \alpha d\Omega + \int D \nabla\phi_i \cdot \nabla\alpha d\Omega + \int (\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha - D \nabla^2 \alpha ) (\tau \mathbf{u} \cdot \nabla \phi_i) d\Omega = 0 $$

Where the additional term that is added effectively consists in the strong form of the equation. We define the stabilization coefficient $\tau$ as:

$$\tau = \left( \left(2 \frac{|\mathbf{u}|}{h} \right)^2+ \left(4 \frac{\mathbf{D}}{h^2} \right)^2 \right)^{-\frac{1}{2}}$$

where $h$ is the size of the cell which is calculated using the diameter of a sphere of equal volume to that of the cell. Consequently $\tau$ has units of time.

Solving the steady-state version of the equation : $$\mathbf{u} \cdot \nabla \alpha = D \nabla^2 \alpha $$

works perfectly fine for any Péclet number. Obviously the GLS term introduces artificial diffusion, but it ensures that the value remain strictly positive and bounded. As the mesh is refined, the solution gets better and better. My issue appears in transient simulations. There, I get loss of positivity and oscillations (overshoots) in my results. It is nothing too bad (nothing like regular Galerkin would give me), but for example, the advection of a step-function without diffusion leads to overshoot of alpha and the bounded character is not preserved.

I feel like something is missing for my GLS stabilization when my equations are transient, but I have not been able to find literature that discuss this problem exactly. What am I missing? Is it because of how I define my $\tau$? For information, I use BDF-type schemes (BDF1 or BDF2) for time integration.

Edit 1: Following the (Finite Element Method : Theory, Implementation and Application by Larson) I understand that a full GLS formulation would require me to add another term that would be :

$$ \int (\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha - D \nabla^2 \alpha ) (- D \nabla^2 \phi_i)$$

but this term is trivially zero for P1 elements, so I do not think it is a necessary component here. At least this is what I saw in the source below:

Larson page

Edit 2

  • Results obtained for steady-state advection-diffusion at Pe=1000 with a coarse mesh

High Péclet - Steady state

  • Advection of a sharp tanh function. Initial condition and after first step (note the overshoot)

Initial condition

After one time step

  • Advection of a step function. Initial condition and after first step (note the overshoot)

Initial step

First step

  • $\begingroup$ The method you are using is SUPG stabilisation, not GLS. Do you get oscillations for all the cases, or just when you have a discontinuous field, for example, square pulse? Did you try with a Gaussian function as the initial condition? It would help in answering if you can post a couple of pictures showing the oscillations. $\endgroup$
    – Chenna K
    Commented Apr 6, 2021 at 20:07
  • $\begingroup$ Well, I agree that a full GLS stabilization would entail adding a supplementary term which would use the laplacian of the test function, but in the case of P1 element, the laplacian of the test function is zero... I have added an edit to my question. $\endgroup$
    – BlaB
    Commented Apr 6, 2021 at 21:34
  • $\begingroup$ Not exactly. In GLS, we have to include the time derivative term as well. The contribution to the stiffness matrix due to stabilisation would be $\tau L^T L$. I don't think this will be of much help, anyway. If you can show the oscillations in the solution you are observing, then it help in identifying the problem. $\endgroup$
    – Chenna K
    Commented Apr 6, 2021 at 23:48
  • $\begingroup$ @ChennaK I added two examples of results. The results obtained for steady-state advection-diffusion at very high Péclet (1000, with a mesh that is far too coarse) and the results obtained with the advection of a tanh step. Consequently, the initial conditions are continuous. They are set by "setting" the nodal values. There is a slight overshoot which surprises me let's say... $\endgroup$
    – BlaB
    Commented Apr 7, 2021 at 3:02
  • 5
    $\begingroup$ These overshoots (and undershoots) are very common with SUPG and SUPG is often paired with a discontinuity capturing (DC) scheme to avoid them. For example see: Tezduyar, Tayfun E., and Y. J. Park. "Discontinuity-capturing finite element formulations for nonlinear convection-diffusion-reaction equations." Computer methods in applied mechanics and engineering 59.3 (1986): 307-325. There are plenty of more recent DC methods as well. $\endgroup$
    – wwfe
    Commented Apr 7, 2021 at 21:29


Your Answer

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