Get the weak form help

General Discussion
21

Yes, this was only a test. A better solution would be to just overwrite the singular terms using the tgv=-1 trick. You will also need to copy the BC from LNSE into the varf for b to get this to work properly.

matrix<complex> OP = LNSE(XXXXXFES,XXXXXFES, tgv = -1, solver=sparsesolver);
matrix<complex> B = b(XXXXXFES,XXXXXFES, tgv = -10, solver=CG);

I also wonder what the alpha in your equations represents. If it is the azimuthal wavenumber it should be set to zero since your formulation is written only for the axisymmetric case.

22

If you need to do non-axisymmetric calculations and the coordinate singularity turns out to be an issue, you might consider using the approach taken by Meliga et al. (Phys. Fluids, 2011), which uses a premultiplication by the radial coordinate to eliminate the singularity. (Note that another factor of r is introduced through the volume integral in cylindrical coordinates). You can find a variant of the resulting weak form here. For simplicity, you can just neglect the boundary integrals in your case.

23

Hello, sorry for the late response,

I do have some papers, here Paper_1 and Paper_2.
So firstly I did a 3D simulation in another program, I performed a time and azimuthally average on the flow components and interpolated the results on this 2D axisymmetric geometry. Now I have to do the linear global stability analysis in 2D.
Indeed alpha is the azimuthal wave number, I set it to 1 as it seems to be the same in these papers.

Regards,
Robert

24

In that case, you will have to include all of the non-axisymmetric terms in the governing equations.

25

Thank you, Chris

After studying more the problem I reached to this form of the 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,[P2,P2,P2,P1]); 	

		XFES Axial, Radial, Tangential, Pressure;
		XXXXXFES[u,v,w,p];
		XXXXXFES<complex>[uy,ux,uz,up];  //Eigenmodes: uaxial, uradial, utheta
		XXXXXFES<complex>[hy,hx,hz,hp];	 //Test functions, axial, radial, tangential directions

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

// uy = Axial Velocity;
// ux = Radial velocity;
// uz = Tangential Velocity;
// up = Static Pressure;

{
ifstream file("DAxial.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(uy, wait=true);
// plot(ux, wait=true);
// plot(uz, wait=true);
// plot(up, wait=true);


real Re = 500;                      	// Reynolds number                       
int  m = 1;						    // Wave number
real shiftr = 0;                    	// Real part of the shift value
real shifti = 0.2;					    // Imaginary part of the shift value
real nu=1/Re;
complex shift = shiftr+shifti*1i;      // Shift value (eigenvalue)     
real eps=1e-12;                  	   // Convergence criteria 
int nev=5;                             // Number of wanted eigenvalues
int ncv=10*nev;                        // Number of eigenvectors
// real epss = 1e-12;					   
int[int] Order = [1,1];
string DataName = "DEMu DEMv ";
macro JJJ x // multiply by "r" which is "x"

macro div(im,u) (dy(u#y)+dx(u#x)+u#x/x+im/x*u#z)// macro for divergence 
macro Grad(im,u) [
					[dy(u#y), dx(u#y),  im/x*u#y ], 
					[dy(u#x), dx(u#x),  im/x*u#x-u#z/x],
					[dy(u#z), dx(u#z),  im/x*u#z+u#x/x ]
				 ] // macro for vecocity gradient tensor
macro D(im,u) [	
				[dy(u#y), 				.5*(dy(u#x)+dx(u#y)),  .5*(im/x*u#y+dy(u#z)) ], 
				[.5*(dy(u#x)+dx(u#y)), 	dx(u#x),				.5*(im/x*u#x-u#z/x+dy(u#z))],
				[.5*(im/x*u#y+dy(u#z)),  .5*(im/x*u#x-u#z/x+dx(u#z)), im/x*u#z+u#x/x]
				] // macro for rate-of-deformation tensor

varf LNSE ([uy,ux,uz,up],[hy,hx,hz,hp])								
=int2d(Th)
		  (
		  (
		  -(v*dx(ux)+ux*dx(v)+u*dy(ux)+uy*dy(v)-2*w*uz/(x+(x==0))+w*ux*1i*m/(x+(x==0)))*hx
		  -(v*dx(uz)+ux*dx(w)+u*dy(uz)+uy*dy(w)+v*uz/(x+(x==0))+ux*w/(x+(x==0))+w*uz*1i*m/(x+(x==0)))*hz	  
		  -(v*dx(uy)+ux*dx(u)+u*dy(uy)+uy*dy(u)+w*uy*1i*m/(x+(x==0)))*hy	
		  +up*div(-1i*m,h)
		  +div(1i*m,u)*hp
		  -2*nu*(D(1i*m,u):D(-1i*m,h))
		  -shift*(uy*hy+ux*hx+uz*hz)
		  +eps*up*hp
		  )*JJJ
		  )	  
		+on(7,12,11, uy=0., ux=0., uz=0.)  //7, 12, 11 -- inlet and walls
		+on(10, uy=0., pe=0.);			           
			
varf b([uy,ux,uz,up],[hy,hx,hz,hp])
=int2d(Th)((uy*hy+ux*hx+uz*hz)*JJJ);

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);
}

However, I see in this example here a boundary condition written as an integral on line on the symmetry axis, which is labeled with “6” in the code (see lines 186-187).

*I changed the wave number name from alpha to m.
Normally, in such cases I should impose on:
Inlet and walls: velocity components of the eigenmodes equal to 0;
Symmetry axis: axial velocity component and pressure equal to 0;

Is there some boundary condition that I am missing?

Best regards,
Rober

26

I’m not sure that I understand your question. The integral in that code is basically equivalent to on() with tgv=1e30. Also, you should definitely not be imposing any boundary conditions on the pressure.