Third derivative

General Discussion
1

Hello developer, I am a beginner about this software.I encountered difficulties in solving the following problem
image
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
image
I wrote their forms based on my own understanding, as shown below
p equation
image
third-order derivative term

image
image1154×61 2.69 KB

At present, the code running display is in a divergent state，The error is as follows：
image
image805×31 1.79 KB

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

2

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.

3

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
image

image
image904×646 49.3 KB

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)

4

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).

5

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)

6

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

7

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.

8

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
image

-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

9

Dear Guangyuan wei,
It was my pleasure to help you on this interesting problem.

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,
d_t u+c^2 grad rho=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.

10

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

11

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

eulerwaves
eulerwaves1512×1172 141 KB

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.

12

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

13

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.
image
Best wishes，
Guangyuan wei