# Third derivative

In the above equation system, I have two questions. I am not very clear about the correct expression of the p-equation, and I am also not clear about how to express the third-order derivative term

I wrote their forms based on my own understanding, as shown below
p equation

third-order derivative term

At present, the code running display is in a divergent state，The error is as follows：

My complete code is as follows，please provide me with some guidance and assistance.
vdw1.edp (1.3 KB)

You are doing a lot wrong here. I suggest you start with a simple Stokes equation solver and gradually increase complexity until you get to your intended problem.

Note that you can derive a weak form for this system that does not involve any 3rd order derivatives. I would suggest you take that approach.

Thank you for taking the time to reply to my questions amidst your busy schedule.I have implemented a finite element method based on feature lines for the example of a cavity.That is a relatively simple and easy code.
I have also made preliminary modifications to the above issues, but its convergence is not good and still diverges.
I feel it is necessary to reiterate my question here

For continuity equations， the weak form is as follows

problem cons(rho,rhot)=
int2d(Th)(rho*rhot/dt)-int2d(Th)((convect([uold,vold],-dt,rhoold)*rhot/dt))+int2d(Th)(rho*dx(uold)*rhot+rho*dy(vold)*rhot)
+int1d(Th,1,2,3,4)(dx(rho)*rhot*N.x+dy(rho)*rhot*N.y);


For the p-equation, I did not handle it in a weak form, but directly calculated it using an expression, as follows

fespace Ph(Th,P1);Ph p;
p=8*rho*sigma/(3-rho)-3*rho^2;


The velocity equation is as follows

problem NSK([u,v],[uu,vv])=
int2d(Th)(rho*(u*uu/dt+v*vv/dt))-int2d(Th)(rho*((convect([uold,vold],-dt,uold)*uu+convect([uold,vold],-dt,vold)*vv)/dt))
+int2d(Th)((1/Re)*(((4/3)*dx(u)-(2/3)*dy(v))*dx(uu)+(dy(u)+dx(v))*dy(uu)+(dx(v)+dy(u))*dx(vv)+((4/3)*dy(v)-(2/3)*dx(u))*dy(vv)))
+int2d(Th)(dx(p)*uu+dy(p)*vv)
+int2d(Th)((1/We)*(dx(rho)*(dxx(rho)+dyy(rho))*uu+dy(rho)*(dxx(rho)+dyy(rho))*vv))
+int2d(Th)((1/We)*(rho*(dxx(rho)+dyy(rho))*dx(uu)+rho*(dxx(rho)+dyy(rho))*dy(vv)))
-int1d(Th,1,2,3,4)((1/We)*(rho*(dxx(rho)+dyy(rho))*uu*N.x+rho*(dxx(rho)+dyy(rho))*vv*N.y))
-int1d(Th,1,2,3,4)((1/Re)*(((4/3)*dx(u)-(2/3)*dy(v))*uu*N.x+(dy(u)+dx(v))*uu*N.y+(dx(v)+dy(u))*vv*N.x+((4/3)*dy(v)-(2/3)*dx(u))*vv*N.y))
+on(1,2,3,4,u=0,v=0);


Especially for the third-order derivative term, I adopted the following weak form

+int2d(Th)((1/We)*(dx(rho)*(dxx(rho)+dyy(rho))*uu+dy(rho)*(dxx(rho)+dyy(rho))*vv))
+int2d(Th)((1/We)*(rho*(dxx(rho)+dyy(rho))*dx(uu)+rho*(dxx(rho)+dyy(rho))*dy(vv)))
-int1d(Th,1,2,3,4)((1/We)*(rho*(dxx(rho)+dyy(rho))*uu*N.x+rho*(dxx(rho)+dyy(rho))*vv*N.y))


The speed always increases until a large number, and I don’t know where my weak form problem lies.Please help me, I will be extremely grateful
vdw1.edp (1.6 KB)

You cannot handle the dispersive term (third order derivative) explicitly as you do, it leads definitely to instability.
To make it implicit you can linearize rho around the value at the old time step. This gives
vdw1.edp (1.8 KB)
In particular you solve in a hidden way a bilaplace type equation in rho.
Your problem is stiff because of the dispersive term, and overall because of the temperature T that you take less than 1. T less than 1 implies that the underlying pressure/velocity wave system is not hyperbolic, hence it is very much unclear if the whole system is well-posed. Then numerical results are not guaranteed at all!
For T greater than 1 I am more confident.
Anyway you get numerical oscillations because of the wave system (to avoid them you would need upwinding, which is not easily done in FreeFem).

Sorry there was a mistake in my code correction: we cannot compute the Laplacian of rho directly because dx(rho) and dy(rho) have jumps between adjacent cells.
We need to compute a projection dxrho of dx(rho) and a projection dyrho of dy(rho) in P1 first. This gives the following fix
vdw1.edp (1.9 KB)

Dear François Bouchut,
Thank you very much for your help with my question. I will carefully consider the issue of not being able to directly calculate the Laplace term of density, but rather calculating the first derivative term of density first.
In addition, I also tried to distinguish between temperatures below and above one, and it is indeed difficult to make the speed converge to a very small value when the temperature is below 1, leading to numerical oscillations. Regarding whether the windward format you mentioned is a velocity based windward format, can you use jump and mean to retrieve its windward format, and where to retrieve its windward format. Is it impossible to implement what you said about being difficult to implement in freeFEM++?
Thank you again for your help. I would greatly appreciate it

My fix of laplacian of rho improves a lot the instabilities for T<1.

Concerning upwinding, this is a principle to avoid oscillations in hyperbolic systems or wave problems. It is implemented in finite difference methods or finite volume methods, but not in finite element methods. A replacement for finite element methods is to put stabilization terms in discontinuous Galerkin methods. An exemple is given in the FreeFem documentation

p. 62 “solution by discontinuous Galerkin FEM”.

Your problem is too complicate, I do not know how to formulate such stabilization.

Dear François Bouchut，
Thank you very much for your reply again. I will consider adding a stable term to the discontinuous finite element method.
I have carefully studied the code you modified last time and have a few questions
Thank you very much for your reply again. I will consider adding a stable term to the discontinuous finite element method.
I have carefully studied the code you modified last time and have a few questions
Firstly, regarding the processing of the pressure equation, is the method you have adopted universal，What kind of thinking is this based on

-int2d(Th)((pression(rhoold)-derpression(rhoold)*rhoold)*(dx(uu)+dy(vv)))-int2d(Th)(derpression(rhoold)*rho*(dx(uu)+dy(vv)))


Secondly, you did not consider boundary integrals when integrating parts, as if you set all boundary integrals to zero. Is this correct.For example, these three items


-int2d(Th)((pression(rhoold)-derpression(rhoold)*rhoold)*(dx(uu)+dy(vv)))-int2d(Th)(derpression(rhoold)*rho*(dx(uu)+dy(vv)))

+int2d(Th)((1/Re)*(((4./3.)*dx(u)-(2./3.)*dy(v))*dx(uu)+(dx(v)+dy(u))*(dx(vv)+dy(uu))+((4./3.)*dy(v)-(2./3.)*dx(u))*dy(vv)))

+int2d(Th)((1/We)*(dx(rhoold)*uu+rhoold*dx(uu)+dy(rhoold)*vv+rhoold*dy(vv))*(dx(dxrho)+dy(dyrho)))
+int2d(Th)((dx(rho)-dxrho)*dxrhot+(dy(rho)-dyrho)*dyrhot)


I would like to thank you again for your help over the past few days. Thank you very much

Dear Guangyuan wei,

Regarding the pressure, I did not follow a general method.
I was guided by the sound wave model which is
d_t rho +div u=0,
The most simple classical treatment of this equation is to solve implicit in both equations.
Your model has this system in it, so that you need similarly to make
implicit the velocity in the continuity equation, and make implicit the pressure in the momentum equation, giving p(rho^{n+1}).
But since we cannot resolve in a simple way such nonlinear term, I put a linearization around rho^n, which gives the formula.

Regarding boundary conditions, since you impose (u,v)=0 on the boundary, there is no way to set any other conditions on the boundary. Indeed the condition (u,v)=0 on the boundary implies that meaningfull test functions (uu,vv) also vanish on the boundary. Looking at any boundary integral you see that it involves either (u,v) or (uu,vv), thus such integral will do nothing.
At the continuous level it means similarly that the problem should be well-posed with the only condition (u,v)=0 and that you cannot impose anything more.
Best wishes,
Francois.

Dear François Bouchut,
Thank you very much these days.I will continue to study further based on your suggestions.
Best wishes,
Guangyuan Wei.

Dear Guangyuan wei,
Another remark: the third-order term induces in your system a Schrodinger type behaviour

Consequently, as in the Schrodinger equation, oscillations are expected in the solution. The observed oscillations may not (only) be caused by the numerical method.
Francois.

Dear François Bouchut,
Thank you for your continued response about the shock you said might be expected, not just a numerical format issue. But some researchers use finite difference to solve this problem, in T<1, the speed converges to a small value, as described below
PRE23_Zaleskin.pdf (2.8 MB)
Best wishes,
Guangyuan Wei

Dear François Bouchut,
I used the spectral method to solve this problem and obtained good results, which may not be a problem with the equation itself.

Best wishes，
Guangyuan wei

Hello,

I am sorry to post on this topic but I am trying to solve the compressible Navier-Stokes equations with Freefem and the continuity and momentum equations are quite the same to the one proposed by M. Guangyuan Wei (except the triple derivative of rho). I have used the following equations :

I have added the equation for energy and the dependance of pressure with respect to temperature and density (ideal gas law), the viscosity as a function of temperature and diffusivity based on a constant Prandtl number. I have tried to write the corresponding weak formulation to solve these equations. In order to try the implementation, I have made a test on the flow around a cylinder at low Reynolds number (212.6).
I would like to know if the finite element space I used are the correct one to solve such kind of equations or if other implementation should be used.

I thank you for your help.

real xmin  = -100;
real xmax  = 200;
real rmax  = 25;

// 1 : axis , 2 : outlet, 3  : cylinder, 4 : inlet, 5 : refinement zone, 6 : top
border a00(t=xmin,-2.5) {x=t; y=0; label=1;}
border a01(t=-2.5,-0.5) {x=t; y=0; label=1;}
border a02(t=pi,0) {x=0.5*cos(t); y=0.5*sin(t); label=3;}
border a03(t=0.5,5.25) {x=t; y=0; label=1;}
border a04(t=0,2) {x=5.25; y=t; label=5;}
border a05(t=5.25,-2.5) {x=t; y=2; label=5;}
border a06(t=2,0) {x=-2.5; y=t; label=5;}
border a07(t=5.25,xmax) {x=t; y=0; label=1;}
border a08(t=0,rmax) {x=xmax; y=t; label=2;}
border a09(t=xmax,xmin) {x=t; y=rmax; label=6;}
border a10(t=rmax,0) {x=xmin; y=t; label=4;}

// Define an initial mesh
int n=3;
mesh Th = buildmesh(a00(20*n)+a01(20*n)+a02(50*n)+a03(20*n)+a04(-20*n)+a05(-40*n)+a06(-20*n)+a07(90*n)+a08(10*n)+a09(60*n)+a10(10*n));

// Save mesh features
{ ofstream file("SOLUT4/coordinates.dat");
for (int j=0;j<Th.nv; j++) {
file << Th(j).x << " " << Th(j).y<< endl;}
}
{ ofstream file("SOLUT4/connectivity.dat");
int nbtriangle = Th.nt;
for (int i=0;i<Th.nt; i++){
file << Th[i][0]+1 << " " << Th[i][1]+1 << " " << Th[i][2]+1 << endl;}
}

real T=1.,dt=0.000005; // total time and time step
real Re=212.6; // Reynolds number (characeristic size =1 and characteristic velocity = 1)
real visco=1.0/Re;
real gamma=1.4;
real rgaz=287.0;
real cp=(gamma*rgaz)/(gamma-1.0); // perfect gas constant
real prandtl=0.7;

//pressure based on the ideal gas law
func real pression(real argrho, real argtemp){
return argrho*rgaz*argtemp;
}

//viscosity as a function of temperature, Sutherland law
func real viscofunc(real argtemp){
return 1.715e-5*(argtemp/273.15)^(3.0/2.0)*((273.15+110.4)/(argtemp+110.4));
}

//conductivity as a function of viscosity considering a constant Prandtl number
func real kappafunc(real argvisco){
return argvisco*cp/prandtl;
}

fespace Rh(Th,P1); Rh rho=1.0,rhot,rhoold,press,temp=300.0,tempt,tempold;
fespace Uh(Th,P1); Uh u=0.0,v=0.0,uu,vv,uold,vold;

problem NSK([rho,u,v,temp],[rhot,uu,vv,tempt])=
int2d(Th)(rho*rhot/dt)-int2d(Th)((convect([uold,vold],-dt,rhoold)*rhot/dt)) +int2d(Th)(rhoold*(dx(u)+dy(v))*rhot) //continuity equation

+int2d(Th)(rhoold*(u*uu+v*vv)/dt)-int2d(Th)(rhoold*(convect([uold,vold],-dt,uold)*uu/dt))-int2d(Th)(convect([uold,vold],-dt,vold)*vv/dt)
-int2d(Th)(pression(rhoold,tempold)*(dx(uu)+dy(vv)))
+int2d(Th)(viscofunc(tempold)*(((4./3.)*dx(u)-(2./3.)*dy(v))*dx(uu)+(dx(v)+dy(u))*(dx(vv)+dy(uu))+((4./3.)*dy(v)-(2./3.)*dx(u))*dy(vv))) // momentum equations

+int2d(Th)(rhoold*cp*(temp*tempt/dt))-int2d(Th)(rhoold*cp*(convect([uold,vold],-dt,tempold)*tempt/dt))
+int2d(Th)(kappafunc(viscofunc(tempold))*(dx(temp)*dx(tempt)+dy(temp)*dy(tempt)))
+int2d(Th)(pression(rhoold,tempold)*dx(u)*tempt+pression(rhoold,tempold)*dy(v)*tempt)
-int2d(Th)(viscofunc(tempold)*(uold*((4./3.)*dx(u)-(2./3.)*dy(v)))*dx(tempt)+viscofunc(tempold)*(vold*(dx(v)+dy(u)))*dx(tempt)+viscofunc(tempold)*(uold*(dx(v)+dy(u)))*dy(tempt)+viscofunc(tempold)*(vold*((4./3.)*dy(v)-(2./3.)*dx(u)))*dy(tempt)) // Temperature equation
+ on(2,4,6,u=1.0,v=0.0,rho=1.0,temp=300.0) + on(3,u=0.0,v=0.0) + on(1,v=0.0); // boundary conditions

int i=0;

for(real t=0.;t<T;t+=dt)
{
uold=u;
vold=v;
rhoold=rho;
tempold=temp;
press=rho*rgaz*temp;
NSK;
if (i%10==0){
{
{ofstream fid("SOLUT4/solut_ux_"+i+".dat");
for (int j=0; j<u[].n; j++)
fid << u[][j] << endl;
}
}

{ofstream fid("SOLUT4/solut_uy_"+i+".dat");
for (int j=0; j<v[].n; j++)
fid << v[][j] << endl;
}
}
{
{ofstream fid("SOLUT4/solut_rho_"+i+".dat");
for (int j=0; j<rho[].n; j++)
fid << rho[][j] << endl;
}
}
{
{ofstream fid("SOLUT4/solut_p_"+i+".dat");
for (int j=0; j<press[].n; j++)
fid << press[][j] << endl;
}
}
cout << "time " << t+dt << " min rho = " << rho[].min << " max rho = " << rho[].max << endl;
i=i+1;
}



Maxime,
You choice of P1 for all unknowns is probably good. But the issue of stability of the resolution is a hard problem because you consider a large Reynolds number (hence small viscosity). In this regime the stability strongly relies of the stability of the hyperbolic part (material waves and sound waves). And it depends on the order of magnitude of the Mach number Ma=|u|/c (c the sound speed) that you expect, and of the timestep you can afford (small timestep is more stable but costs a lot). A reasonable timestep is dt |u|\sim h. If you can afford c dt\sim h then an explicit scheme can be stable. Otherwise if c dt>>h you have definitely to implicit the pressure.
You should first look at a simpler problem where you replace the temperature equation by a polytropic law p=C\rho^\gamma. Then the issue is the pressure term in the momentum equation: explicit (as you write in your code) or implicit. An exemple of implicit scheme is the one I wrote above (linearisation of p around rhoold).
Validated schemes for the compressible inviscid problem are more complicate (finite volume methods), for example
A low cost semi-implicit low-Mach relaxation scheme for the full Euler equations. J. Scientific Comput., 83:24, 2020. DOI: 10.1007/s10915-020-01206-z
thus making your finite element scheme work with only little viscosity is difficult, at best it would be stable but with lot of oscillations.

Hello François,

Thank you very much for your message. I was thinking to stay explicit in time, but as you said the timestep ensuring stability will be decreased. Do CFL number is relevant as in other discretization methods. For example setting a CFL number of 1 lead in incompressible to set
dt=dx/u and for compressible equations as in the present case to dt=dx/(u+c). Maybe by staying explicit in time and using ffddm approach of Freefem to parallelize the solving of the problem would alleviate the problem of having a small dt (but for the moment I would be happy if it would work sequentially)

So, at first you would not consider the energy equation but would use the polytropic law to express pressure in terms of the density if I understood well. I can definitely do that, thank you for the advise. In fact, I was thinking it would be possible to have the full compressible Navier-Stokes equations as Meliga et al. used this set of equations with Freefem to do the study in the paper attached. But there may be a lot of work of validation for the stability of the scheme so I can start with the solutions you proposed me.

Thank you very much once again.
meliga2010rp-1pp.pdf (5.9 MB)

In the paper you mention they consider a steady base flow, that they compute with an iterative method. Then they compute the time dependent linearized flow around this base flow. This is possible to do that, but an am not at all an expert on this. What I was commenting is the question of solving the time dependent (transient) nonlinear problem, which is more difficult for the stability because you need to define a scheme that is at the same time implicit and nonlinear. This is what I understood you wanted to do.
An additional comment on the scales is that if you want the viscosity to stabilize the hyperbolic part you need that
(|u|+c)h/\mu<1
or in other words
(1+\frac{1}{Ma})Re\frac{h}{L}<1.
Therefore if you can affort such a fine mesh, you could use simple explicit \nabla p in the momentum equation and \nabla\cdot (up) in the temperature equation (as you tried) without care. And in this case I think you do not have to satisfy any CFL condition. For the dissipation term in the temperature equation, you can avoid to discretize it by writing the equation on the total energy \rho|u|^2/2+\rho c_p T ie you discretize this equation instead of the one on T.

Hello François,

thanks for the explanations. I tought Meliga et al. were considering and solving the full unsteady compressible Navier-Stokes equations (non linear) as shown in the set of eq. 2.1a 2.1d since there is the non linear convective term on the left hand side of momentum equations \rho \nabla u .u if I am not mistaking (but I agree that it should be the term \nabla.(up) in the energy equation and this has been transformed into p \nabla.u).

But also, I am not too sure about what I am saying and you are surely right since the case is at low Reynolds number, the solution is laminar and steady, and so an iterative process you presented should be enough. It is just that in the paper it seems that they solved the unsteady compressible equations with freefem to converge towards a steady solution. Then, this solution was used as a baseflow for the linearized compressible Navier-Stokes equations to obtain the natural instabilities associated to this base flow.

What I was thinking for solving the unsteady compressible Navier-Stokes equations is to be explicit in time (maybe with a simple Euler scheme for the temporal terms) and for the non linear terms use the method of characteristics (convect function of Freefem). Then, because the condition of stablity leads to low values of \Delta t, use the ffddm with petsc to parallelize it in order to alleviate this low time step. But If I am not wrong, this is not possible as ffddm needs to invert a global problem (matrix system) which leads to be implicit in time, or maybe could you tell me if I am wrong or not and whether you would have reather used an implicit approach in time ?

I thank you for the advice about how to handle the \nabla p and \nabla. (up) terms as well as the way to write the trasnport equation on the total energy rather than temperature, I will do it this way.

Thank you.

It you use an explict scheme as you describe for the hyperbolic part (but implicit for the viscous part) with small timestep satisfying a CFL condition, you could hopefully parallelize both the linear system to solve (viscous part) and the explicit part. But you would nevertheless get oscillations. To avoid oscillations you need either to use specific schemes for the hyperbolic part, involving approximate Riemann solvers, or use a fine enough mesh (as described in my previous message) so that at the mesh level the viscosity smoothes out the oscillations. If Re is moderate and Ma is not too small this last option is excellent because in this case hopefully you do not need to satisfy any CFL condition.