PSPG/SUPG

Asked by Marcos Samudio on 2012-04-02

Hello,

I am trying to implement a stabilised Coupled Navier-Stokes solver with the k-epsilon turbulence model.

I have formulated the PSPG and SUPG stabilisations in the weak form and got different expressions than those from NSCoupled in the CBC package. I am using CG 1 for both velocity and pressure.

I get the following:

nut_ = turbulent diffusion
nu = viscous diffusion
eps and esu are constants

FOR PSPG:
    h = CellSize(mesh)
    eps = eps*h
    F_pspg = eps*inner(grad(q), grad(p))*dx + eps*inner(grad(q), dot(grad(u), u_))*dx - <=pressure gradient and convective term
                      eps*inner(grad(q), f)*dx + eps*inner(grad(q), dot(grad(u), grad(nut_)))*dx <=body forces, non-constant viscosity term

    for bc in self.bcs:
       if bc.type() in ('ConstantPressure','Outlet','VelocityInlet','Slip'):
           F = F + eps*nu*inner(grad(q), dot(grad(u), n))*ds(bc.bid) <=from integration of viscous term. nut_ is 0 on walls( there is a term missing for other BCs, but I only have wall bcs right now...)

FOR SUPG:
        h = CellSize(mesh)
        esu = esu*h
        stab = dot(grad(v), u_)
        nu = kwargs.get('nuM', None)(0)

        F = esu*inner(stab, grad(p))*dx + esu*inner(stab, dot(grad(u), u_))*dx <=pressure gradient and convective term
                + esu*(nu+nut_)*inner(grad(stab), grad(u))*dx - esu*inner(stab,f)*dx <=integration by parts of viscous term, body forces
                + esu*inner(stab, dot(grad(u), grad(nut_)))*dx <=from integration by parts, with non-constant nut_

        for bc in self.bcs:
            if bc.type() in ('ConstantPressure','Outlet','VelocityInlet','Slip'):
                F = F - esu*nu*inner(stab, dot(grad(u), n))*ds(bc.bid) <=from integration of viscous term. nut_ is 0 on walls( there is a term missing for other BCs, but I only have wall bcs right now...)

There is also a minus sign in the PSPG formulation from NSCoupled, which I do not quite get...

When I try running the turbulent channel demo with my code the simulation runs for about 40 iterations and then I get a 'nan' crash. If I make dot(grad(u), grad(nut_) explicit, taking u=u_, the simulation goes on without problems.

Is there something wrong? Any help will be very appreciated!

Thanks,

Marcos

Question information

Language:
English Edit question
Status:
Solved
For:
CBC.PDESys Edit question
Assignee:
No assignee Edit question
Solved by:
Marcos Samudio
Solved:
Last query:
Last reply:

Hi Marcos,

On 2 April 2012 22:50, Marcos Samudio
<email address hidden>wrote:

> New question #192451 on CBC.PDESys:
> https://answers.launchpad.net/cbcpdesys/+question/192451
>
> Hello,
>
> I am trying to implement a stabilised Coupled Navier-Stokes solver with
> the k-epsilon turbulence model.
>
> I have formulated the PSPG and SUPG stabilisations in the weak form and
> got different expressions than those from NSCoupled in the CBC package. I
> am using CG 1 for both velocity and pressure.
>
> I get the following:
>
> nut_ = turbulent diffusion
> nu = viscous diffusion
> eps and esu are constants
>
> FOR PSPG:
> h = CellSize(mesh)
> eps = eps*h
> F_pspg = eps*inner(grad(q), grad(p))*dx + eps*inner(grad(q),
> dot(grad(u), u_))*dx - <=pressure gradient and convective term
> eps*inner(grad(q), f)*dx + eps*inner(grad(q),
> dot(grad(u), grad(nut_)))*dx <=body forces, non-constant viscosity term
>
>
This is different from NSCoupled. Did you try the original form? It works
for me. I derived it using two consecutive integrations by parts. How do
you arrive at the form above? Are you using the Laplace form (no grad(u).T
term) of NS? That equation is not valid for variable nut.

Also, I derive my equation be taking the divergence of NS and then multiply
by q. Some people I've seen simply multiply be grad(q). There's a ds term
in difference.

> for bc in self.bcs:
> if bc.type() in ('ConstantPressure','Outlet','VelocityInlet','Slip'):
> F = F + eps*nu*inner(grad(q), dot(grad(u), n))*ds(bc.bid)
> <=from integration of viscous term.
> nut_ is 0 on walls( there is a term missing for other BCs, but I only have
> wall bcs right now...)
>
> FOR SUPG:
> h = CellSize(mesh)
> esu = esu*h
> stab = dot(grad(v), u_)
> nu = kwargs.get('nuM', None)(0)
>
> F = esu*inner(stab, grad(p))*dx + esu*inner(stab, dot(grad(u),
> u_))*dx <=pressure gradient and convective term
> + esu*(nu+nut_)*inner(grad(stab), grad(u))*dx -
> esu*inner(stab,f)*dx <=integration by parts of viscous term, body forces
> + esu*inner(stab, dot(grad(u), grad(nut_)))*dx
> <=from integration by parts, with
> non-constant nut_
>
> for bc in self.bcs:
> if bc.type() in
> ('ConstantPressure','Outlet','VelocityInlet','Slip'):
> F = F - esu*nu*inner(stab, dot(grad(u), n))*ds(bc.bid)
> <=from integration of viscous term. nut_ is 0 on
> walls( there is a term missing for other BCs, but I only have wall bcs
> right now...)
>
> There is also a minus sign in the PSPG formulation from NSCoupled, which I
> do not quite get...
>

What do you mean?

>
> When I try running the turbulent channel demo with my code the simulation
> runs for about 40 iterations and then I get a 'nan' crash. If I make
> dot(grad(u), grad(nut_) explicit, taking u=u_, the simulation goes on
> without problems.
>

Is it the same nan problem as previously?

Mikael

>
> Is there something wrong? Any help will be very appreciated!
>
> Thanks,
>
> Marcos
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>

Marcos Samudio (marcossamudio) said : #2

Hello Mikael,

I have tried the original form from NSCoupled, and it seemed to work, but I could not derive it on my own. That is why I decided to try this form, which I got from my calculations.

You are right about the Laplace form, I will have to consider that, as viscosity is not constant!

I do not quite understand what you did to get the NSCoupled form. You say you took NS divergence and then multiplied by q? Because what I did was weighting the residual of NS by grad(q), just the way I saw in, for example, Donea-Huerta. Is there a reason for doing it that way? What are the advantages? And why do you weight by q instead of grad(q)?

About the PSPG sign, weighting the residual yields the same signs as in the NS equation for the PSPG terms, but in the NSCoupled formulation I see opposite signs. I suppose that comes from the NS divergence formulation you talked about, but I am not sure.

About the nan error, I would rather trying to reformulate PSPG before drawing out conclusions!

I will be feeding back soon and see if I can figure out what is going wrong!

Thanks again,

Marcos

Den Apr 3, 2012 kl. 3:30 PM skrev Marcos Samudio:

> Question #192451 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/192451
>
> Status: Answered => Open
>
> Marcos Samudio is still having a problem:
> Hello Mikael,
>
> I have tried the original form from NSCoupled, and it seemed to work,
> but I could not derive it on my own. That is why I decided to try this
> form, which I got from my calculations.
>
> You are right about the Laplace form, I will have to consider that, as
> viscosity is not constant!
>
> I do not quite understand what you did to get the NSCoupled form. You
> say you took NS divergence and then multiplied by q? Because what I did
> was weighting the residual of NS by grad(q), just the way I saw in, for
> example, Donea-Huerta. Is there a reason for doing it that way? What are
> the advantages? And why do you weight by q instead of grad(q)?

Taking the divergence of the equation is an exact mathematical manipulation of the equation. I'm from a physics background more than mathematics and it makes more sense to me than to simply multiply by (piecewise constant) grad(q). In the end the two approaches are different because the first approach gives the second plus a ds term:

inner(div(NS), q) = - inner(NS, grad(q))*dx + inner(NS*n, q)*ds

I believe the first approach is more correct, but as you say, other people are simply multiplying by grad(q).

There are just a few steps required in the derivation of a computable form of inner(div(NS), q)*dx (a few integrations by parts), but I haven't written it down anywhere electronically. If you struggle I can find the time to write it down.

>
> About the PSPG sign, weighting the residual yields the same signs as in
> the NS equation for the PSPG terms, but in the NSCoupled formulation I
> see opposite signs. I suppose that comes from the NS divergence
> formulation you talked about, but I am not sure.

I guess so, there's a sign difference in the form above.

Mikael

>
> About the nan error, I would rather trying to reformulate PSPG before
> drawing out conclusions!
>
> I will be feeding back soon and see if I can figure out what is going
> wrong!
>
> Thanks again,
>
> Marcos
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Marcos Samudio (marcossamudio) said : #4

Ok, I see what you did. Then you get the following:

PSPG_stabilization = inner(NS, grad(q))*dx = - inner(div(NS), q) + inner(NS*n, q)*ds

which is the term you add to the weak form, right?

But when doing integration by parts of the divergence term, the signs are reversed, and so you would get positive signs again. I do not know if I am missing something or do not quite get it...

Ok, I did the derivation in detail. See attachment. There are probably some
errors.

I actually never had time to go any further with this, I was just happy it
worked at one stage. It would be very interesting to verify the method
though. I'm not sure anyone has done it exactly like this before. What are
your current plans with this?

Mikael

On 3 April 2012 20:25, Marcos Samudio
<email address hidden>wrote:

> Question #192451 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/192451
>
> Status: Answered => Open
>
> Marcos Samudio is still having a problem:
> Ok, I see what you did. Then you get the following:
>
> PSPG_stabilization = inner(NS, grad(q))*dx = - inner(div(NS), q) +
> inner(NS*n, q)*ds
>
> which is the term you add to the weak form, right?
>
> But when doing integration by parts of the divergence term, the signs
> are reversed, and so you would get positive signs again. I do not know
> if I am missing something or do not quite get it...
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>

Marcos Samudio (marcossamudio) said : #6

Mikael,

I do not see any links to attachments...

My goal is to simulate turbulent flow with the SUPG and PSPG stabilizations. I have been struggling with the variational forms today and did not get good results...I will give it a try tomorrow.

I forgot, attachments don't work. I'm adding the file derivation.pdf to the
repository instead.

Mikael

On 4 April 2012 05:30, Marcos Samudio
<email address hidden>wrote:

> Question #192451 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/192451
>
> Marcos Samudio posted a new comment:
> Mikael,
>
> I do not see any links to attachments...
>
> My goal is to simulate turbulent flow with the SUPG and PSPG
> stabilizations. I have been struggling with the variational forms today
> and did not get good results...I will give it a try tomorrow.
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>

Marcos Samudio (marcossamudio) said : #8

The term you are adding to the weak form is inner(div(NS), q), am I right?

I was thinking you were using PSPG stabilization but it looks like you are not. Adding inner(grad(q), NS, as I was trying to do, would reverse signs in the formulation.

How did you come up with this? I have not seen this before...at least if I am right in my assumptions!

On Wednesday, 4 April 2012, Marcos Samudio wrote:

> Question #192451 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/192451
>
> Marcos Samudio posted a new comment:
> The term you are adding to the weak form is inner(div(NS), q), am I
> right?

Yes. With FEM manipulations (int by parts) to obtain a form that works with
CG=1.

>
> I was thinking you were using PSPG stabilization but it looks like you
> are not. Adding inner(grad(q), NS, as I was trying to do, would reverse
> signs in the formulation.

But adding -inner(grad(q), NS) would not. The sign is arbitrary I guess.
Could be significant for iterative solvers.

>
> How did you come up with this? I have not seen this before...at least if
> I am right in my assumptions!

Taking the divergence of NS is a common way to obtain an equation for the
pressure. To use it for stabilizing FEM solvers with variable coefficient I
have not seen anywhere else.

By the way, when you're multiplying NS with grad(q) I don't think the
transient term will be zero, as it is when you take the divergence.

Mikael

>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>

I do not know if the sign is arbitrary...things go wrong when you reverse the stabilization parameter sign (eps).

Things are much clearer now, as I finally get what you are doing! I am trying to write the PSPG stabilization, and will submit my results so we can discuss it if you like.

About the transient term, when integrated by parts, inner(q_,i , u_i,t) yields velocity divergence variation (null) and a boundary term which is zero if boundaries are steady on time.

Marcos

Mikael

I have transcribed the derivation of the PSPG terms I am using. It is available online at https://docs.google.com/open?id=0B9qsTlUZxTpdNXVwejhXaW5SU1dMNjRXU2ZwbVh0QQ

Anyway, this is not working properly, at least with constant nu (setting nu_T to zero) but I do not know why.

I am going to start working in the SUPG stabilization as well, see if I can make it work (for constant nu it seems to be already working).

Thanks,

Marcos

Mikael,

I got my PSPG and SUPG stabilizations working today (only for low Re though, I think I must add SUPG to k and epsilon transport as well, to make them more stable).

I found what the problem was. In the Steady_NS_form definition, you did
F = inner(v, NS)*dx - inner(q, div(u))*dx

When you added the stabilization term you defined, inner(div(NS), q), the result was:

F = inner(v, NS)*dx - inner(q, div(u))*dx - inner(grad(q), NS)*dx + inner(NS*n, q)*ds

By doing this, you formulate a consistent PSPG stabilization (it is like using -q instead of q as weighting function), which explains why reversing the sign of my stabilization terms got things working. I did not have any boundary terms in my particular case (lid driven cavity), so the boundary integral did not affect the results.

The SUPG and PSPG formulations I am using are these:
PSPG: https://docs.google.com/open?id=0B9qsTlUZxTpdbi1nQVU3QnlUb1d6WnRRX3liUmV6QQ
SUPG: https://docs.google.com/open?id=0B9qsTlUZxTpdYUxvYktFc2pUVy1INkdUSllkYzk1UQ

If there are any corrections to be made, they are more than welcome!
Thanks again for all your help!

Marcos