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:
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:
where cZ and p are the ue and pe terms in my code.
Thank you again,
Regards,
Robert
Test.rar (2.2 MB)
