Get the weak form help

Dear forum,
Once again I ask for your kindness to guide me through a problem. In the figure attached (I didn’t find a better way to show the equations) there is a set of 4 equations that I am trying to solve, more precisely, the stability equations written in the cylindrical coordinates for an axisymmetric flow. As I don’t get the expected results I guess there is something wrong with my weak formulation or some transcript problem.

As a short nomenclature, what you see there with u _ R,theta,Z represents the radial, azimuthal and axial coordinates of the eigenmodes/velocities. h _ R,theta,Z are the test functions for the three directions and q is the test function for mass continuity equation.

So here are my questions:

  1. I highlighted with red some terms which are the subject of a reduction in the weak formulation, is this the correct way to bring them to the weak form when we refer to cylindrical coordinates?

  2. Do I have to multiply all the terms by R, which in FreeFEM++ transcript become x?

Thank you in advance,

Robert

I think in the StabFem library, you might be able to find the solution:

Your equations look fine from what I can tell. Yes, the integral in cylindrical coordinates introduces a factor of r due to the definition of the differential volume. Please feel free to follow up with any more questions!

Thank you @aszaboa for the answer, I am aware of StabFEM librari but unfortunately I didn’t find an example to help me in this case/or I didn’t understand all the script.

@cmd thank you, I added the “r” term which I missed from the DV=rdxdr. I attach the new variant of the equations, I am still not sure if this is the way to apply the Green/Gauss method to reduce the order of the laplacian highlighted with red and for the pressure gradient in cylindrical coordinate.

Anyway, I wrote the following code:

load "UMFPACK64"
include "macro_ddm.idp"		
load "gmsh"					
load "iovtk"

mesh Th=gmshload("Test.msh");
// plot(Th, wait=true);

fespace XFES(Th, P1);
fespace XMFES(Th, P2);
fespace XXXXXFES(Th,[P1,P1,P1,P1]); 	

		XFES Axial, Radial, Tangential, Pressure;
		XXXXXFES[u,v,w,p];
		XXXXXFES<complex>[ue,ve,we,pe];  //Eigenmodes: uaxial, uradial, utheta
		XXXXXFES<complex>[hy,hx,hz,q];	 //Test functions, axial, radial, tangential directions

		XFES EMur, EMui;
func E = [P1, P1, P1, P1];				

// u = Axial Velocity;
// v = Radial velocity;
// w = Tangential Velocity;
// p = Static Pressure;

{
ifstream file("Axial.txt");
file >> Axial[];
}
{
ifstream file("Radial.txt");
file >> Radial[];
}
{
ifstream file("Tangential.txt");
file >> Tangential[];
}
{
ifstream file("Pressure.txt");
file >> Pressure[];
}

// Combine the files above in one finite element space - XXXXXFES instead of using XFES for each of them
[u,v,w,p] = [Axial, Radial, Tangential, Pressure]; //Time and azimuthal averaged flow components
// plot(u, wait=true);
// plot(v, wait=true);
// plot(w, wait=true);

// plot(u, wait=true);
// plot(v, wait=true);
// plot(w, wait=true);
// plot(p, wait=true);

real Re = 500;                      	// Reynolds number                       
int alpha = 1;						    // Wave number
real shiftr = 0;                    	// Real part of the shift value
real shifti = 0.2;					    // Imaginary part of the shift value

complex shift = shiftr+shifti*1i;      // Shift value (eigenvalue)     
real eps=1e-6;                  	   // Convergence criteria 
int nev=3;                             // Number of wanted eigenvalues
int ncv=10*nev;                        // Number of eigenvectors
						   
int[int] Order = [1,1];
string DataName = "DEMu DEMv ";

varf LNSE ([ue,ve,we,pe],[hy,hx,hz,q])								
=int2d(Th)
		  (
		  -(v*dx(ve)+ve*dx(v)+u*dy(ve)+ue*dy(v)-2*w*we/x+w*ve*1i*alpha/x)*hx*x
		  +pe*dx(hx)*x
		  +(1/Re)*((1/x)*dx(ve)-ve/(x^2)-2*we*1i*alpha/(x^2)-ve*(alpha^2)/(x^2))*hx*x
		  -(1/Re)*(dx(ve)*dx(hx)+dy(ve)*dy(hx))*x
		  -(v*dx(we)+ve*dx(w)+u*dy(we)+ue*dy(w)+v*we/x+ve*w/x+w*we*1i*alpha/x)*hz*x
		  -pe/x*1i*alpha*hz*x
		  +(1/Re)*((1/x)*dx(we)-we/(x^2)+2*ve*1i*alpha/(x^2)-we*(alpha^2)/(x^2))*hz*x
		  -(1/Re)*(dx(we)*dx(hz)+dy(we)*dy(hz))*x
		  -(v*dx(ue)+ve*dx(u)+u*dy(ue)+ue*dy(u)+w*ue*1i*alpha/x)*hy*x
		  +pe*dx(hy)*x
		  +(1/Re)*((1/x)*dx(ue)-ue*(alpha^2)/(x^2))*hy*x
		  -(1/Re)*(dx(ue)*dx(hy)+dy(ue)*dy(hy))*x
		  +(dx(ve)+we*1i*alpha/x+ve/x+dy(ue))*q*x
		  -shift*(ue*hy+ve*hx+we*hz)*x
		  	  		  )
 	  	  
		+on(7,12,11, ue=0., ve=0., we=0.)  //7, 12, 11 -- inlet and walls
		+on(10, ue=0.);			           //symmetry axis, pressure needed?? pe=0.
			
varf b([ue,ve,we,pe],[hy,hx,hz,q])
=int2d(Th)((ue*hy+ve*hx+we*hz)*x);

matrix<complex> OP = LNSE(XXXXXFES,XXXXXFES, solver=sparsesolver);
matrix<complex> B = b(XXXXXFES,XXXXXFES, solver=CG);
complex[int] EigenVal(nev); 		
XXXXXFES<complex>[int] [EigenVecu,EigenVecv,EigenVecw,EigenVecp](nev);

int k;
int i;
k = EigenValue(OP,B,sigma=shift,value=EigenVal,nev=nev,vector=EigenVecu,ncv=ncv,tol=eps);
//PLOTS//
for (i=0; i<nev; i++)
{
 EMur=real(EigenVecu[i]);		//Real part of the streamwise component direct eigenmode
 EMui=imag(EigenVecu[i]);		//Imaginary part of the streamwise component direct eigenmode
savevtk("EigenModes"+i+".vtu", Th, [EMur,EMui], dataname = DataName, order = Order);
}

I also attach the mesh, pressure and velocity files. As you can see the real part of the direct eigenmodes (x component of EMur if you open it in paraview) is quite strange/blurry and I was expecting to get something like in the following picture:
image

Could you tell me if there is a typo error in the code or just a mathematical issue? As for the boundary conditions I should impose the following:
image
where cZ and p are the ue and pe terms in my code.

Thank you again,

Regards,
Robert

Test.rar (2.2 MB)

One thing I can definitely tell you off the bat is that your mixed finite element space is not inf-sup stable. You have to use an inf-sup stable FEspace – something like Taylor–Hood ([P2,P2,P2,P1]) or MINI ([P1b,P1b,P1b,P1]) elements. You may also have to algebraically simplify your weak form (i.e. enforce that x*(1/x)=1) so that you dont divide everything by zero along the coordinate axis.

Also, it does not look like you integrated the viscous term by parts correctly. Note that you are integrating by parts using vectorial test and trial functions – you cannot treat the x, r, and theta momentum equations separately when integrating by parts.

Thank you, I changed to P2 elements type and enforced the term x*(1/x)=1, huge differences indeed. As with the mathematical part, I surely miss something. Do you refer to the viscous terms of second order partial derivative, highlighted in red?


Regards,
Robert

Not just the parts highlighted in red – the entire viscous term. For example, see slide 9 here. Note that in your case, the gradient and divergence operators must be the correct ones for a cylindrical coordinate system.

Thank you, Chris for your patience. So I took separately only the viscous term and this is what I’ve got.


I guess there is still something misunderstood in the equations or a typo in the code

Regards,
Robert

I think the main issue is that you are still integrating by parts incorrectly. You have to integrate the VECTOR Laplacian by parts. What you are doing doesn’t make sense since you are treating the viscous term as three separate instances of a scalar Laplacian.

Wouldn’t that be simply
image
where the term involving the partial derivative in respect to azimuthal direction is null?

No, you are integrating a vector-valued function by parts with respect to a vector-valued test function.

Thank you for your guidance, but I still miss something to get the correct weak formulation. Is it possible to show me the solution?

Regards,
Robert

Here is one way to write it for the axisymnmetric case. Here, x and y are along the axial and radial directions, respectively. [vax,vra,vth] is the velocity test function and [uax, ura, uth] is the velocity trial function.

1/Re*(y*(dx(vax)*dx(uax) + dy(vax)*dy(uax) + dx(vra)*dx(ura) + dy(vra)*dy(ura) + dx(vth)*dx(uth) + dy(vth)*dy(uth)) + (vth*uth + vra*ura)/y)

Thank you for your time!

So, basically the problem is the following:

And I wrote the following script:

load "UMFPACK64"
include "macro_ddm.idp"		
load "gmsh"					
load "iovtk"

mesh Th=gmshload("Test.msh");
// plot(Th, wait=true);

fespace XFES(Th, P1);
fespace XMFES(Th, P2);
fespace XXXXXFES(Th,[P2,P2,P2,P1]); 	

		XFES Axial, Radial, Tangential, Pressure;
		XXXXXFES[u,v,w,p];
		XXXXXFES<complex>[ue,ve,we,pe];  //Eigenmodes: uaxial, uradial, utheta
		XXXXXFES<complex>[hy,hx,hz,q];	 //Test functions, axial, radial, tangential directions

		XFES EMur, EMui;
func E = [P2, P2, P2, P1];				

// u = Axial Velocity;
// v = Radial velocity;
// w = Tangential Velocity;
// p = Static Pressure;

{
ifstream file("Axial.txt");
file >> Axial[];
}
{
ifstream file("Radial.txt");
file >> Radial[];
}
{
ifstream file("Tangential.txt");
file >> Tangential[];
}
{
ifstream file("Pressure.txt");
file >> Pressure[];
}

// Combine the files above in one finite element space - XXXXXFES instead of using XFES for each of them
[u,v,w,p] = [Axial, Radial, Tangential, Pressure]; //Time and azimuthal averaged flow components
// plot(u, wait=true);
// plot(v, wait=true);
// plot(w, wait=true);

// plot(u, wait=true);
// plot(v, wait=true);
// plot(w, wait=true);
// plot(p, wait=true);

real Re = 500;                      	// Reynolds number                       
int alpha = 1;						    // Wave number
real shiftr = 0;                    	// Real part of the shift value
real shifti = 0.2;					    // Imaginary part of the shift value

complex shift = shiftr+shifti*1i;      // Shift value (eigenvalue)     
real eps=1e-6;                  	   // Convergence criteria 
int nev=3;                             // Number of wanted eigenvalues
int ncv=10*nev;                        // Number of eigenvectors
						   
int[int] Order = [1,1];
string DataName = "DEMu DEMv ";

varf LNSE ([ue,ve,we,pe],[hy,hx,hz,q])								
=int2d(Th)
		  (
		  -((v*dx(ve)+ve*dx(v)+u*dy(ve)+ue*dy(v))*hx*x+(w*ve*1i*alpha-2*w*we)*hx)
		  -((v*dx(we)+ve*dx(w)+u*dy(we)+ue*dy(w))*hz*x+(v*we+ve*w+w*we*1i*alpha)*hz)
		  -((v*dx(ue)+ve*dx(u)+u*dy(ue)+ue*dy(u))*hy*x+(w*ue*1i*alpha)*hy)
		  +pe*dx(hx)*x+pe*dx(hy)*x-pe/x*1i*alpha*hz*x
		  -(1/Re)*(x*(dx(ve)*dx(hx)+dy(ve)*dy(hx)+dx(we)*dx(hz)+dy(we)*dy(hz)+dx(ue)*dx(hy)+dy(ue)*dy(hy)))
		  -(1/Re)*(we*hz+ve*hx)/x 
		  -shift*(ue*hy+ve*hx+we*hz)*x
		  )  	  
		+on(7,12,11, ue=0., ve=0., we=0.)  //7, 12, 11 -- inlet and walls
		+on(10, ue=0., pe=0.);			           //symmetry axis, pressure needed?? pe=0.
			
varf b([ue,ve,we,pe],[hy,hx,hz,q])
=int2d(Th)((ue*hy+ve*hx+we*hz)*x);

matrix<complex> OP = LNSE(XXXXXFES,XXXXXFES, solver=sparsesolver);
matrix<complex> B = b(XXXXXFES,XXXXXFES, solver=CG);
complex[int] EigenVal(nev); 		
XXXXXFES<complex>[int] [EigenVecu,EigenVecv,EigenVecw,EigenVecp](nev);

int k;
int i;
k = EigenValue(OP,B,sigma=shift,value=EigenVal,nev=nev,vector=EigenVecu,ncv=ncv,tol=eps);
//PLOTS//
for (i=0; i<nev; i++)
{
 EMur=real(EigenVecu[i]);		//Real part of the streamwise component direct eigenmode
 EMui=imag(EigenVecu[i]);		//Imaginary part of the streamwise component direct eigenmode
savevtk("EigenModes"+i+".vtu", Th, [EMur,EMui], dataname = DataName, order = Order);
}

And I get the following error:

I also attach the files that I used. Is it possible to tell me where the error might come from?

Regards,
Robert
Test.rar (488.5 KB)

Did you put the edp in the rar file? The copy from screen seems to drop asterisks
and make other problems on my Chrome.

Hello,

No, I attach a new .rar with the .edp file too

Regards,
Robert
Test.rar (488.5 KB)

I wonder if there is an issue with your varf and the singularity at r=0. Try replacing the 1/y with 1/(y+(y==0)).

Note that one could probably find a better way to avoid dividing by zero.

Thank you very much, Chris. It seems that did the trick, it is working now. The results are however different than what I expected, I will keep looking to the equations/boundary conditions/perform some mesh sensitivity analysis to find the problem.

Regards,
Robert

What results did you expect "? Do you have a paper to cite?
AFAICT that change was just intended to determine if the divide by
zero was causing the immediate problem. It may require some reformulation
to get a useful result. On a quick look of the literature this seems
to be a bit of an open topic. Apparently there are some frequency
domain approaches. These may help with some singularities
but are still difficult with non-linear terms.