Boundary conditions in eigenvalue problem

[maximefio (Maxime)]

Hello,

I am currently doing a 1 dimensional linear stabiity solver (for fluid mechanics). To do so I solve an eigenvalue problem as AX = lambda BX where A and B are terms from the linearized Navier-Stokes equations and of the mean flow and X is the small perturbation that is unknown (u, v, w, p).

In these matrices A and B must be added boundary conditions, for some case only Neumann or Dirichlet boundary conditions need to be imposed (case m = 0 and m>1 in page 3). For these two cases, Freefem works perfectly well.However, for the case m=±1, the first perturbation u depends on the second perturbation v as 2u’(0)+mv’(0)=0.

Would there be a way to impose such kind of boundary condition if Freefem ? I have put the case that work for m = 0. In the code F stands for u, G for v, H for w and Pr for p.

real ymin = 0.;
real ymax = 1.;
border yline(t=ymin,ymax){x=0.;y=t;}
int np = 1000;

meshL Th=buildmeshL(yline(np));

fespace Xh(Th, P2);
fespace Mh(Th, P1);

Xh W=1-y^2, V=0, U=0, Re=100.0, Visco=1/Re, PrMean = 0.;
Xh alpha=10.0, m=0.0;

fespace Vh(Th, [P2, P2, P2, P1]);
Vh [F, G, H, Pr], [X1, X2, X3, X4];

varf opLHS([F, G, H, Pr], [X1, X2, X3, X4]) =
int1d(Th)( -(1iy^2/Re)dy(F)dy(X1)-((1i/Re)((m^2+1)+y^2alpha^2))FX1+y^2alphaWFX1-(1i)(y/Re)Fdy(X1)-(21im/Re)GX1+y^2Prdy(X1)
+(1i)(-(y^2/Re)dy(G)dy(X2))-((1i)/Re)((m^2+1)+y^2alpha^2)GX2+y^2alphaWGX2-(1i)(y/Re)Gdy(X2)-(2m(1i)/Re)FX2+myPrX2
-(1i)(y^2/Re)dy(H)dy(X3)-(1i)(y/Re)Hdy(X3)+y^2alphaWHX3-(1i/Re)(m^2+alpha^2y^2)HX3+y^2dy(W)FX3+1iy^2alphaPrX3
-X4*(F+mG+alphayH)+yF*dy(X4)
)
+on(1, F=0, G=0)
+on(2, F=0, G=0, H=0);

// get the RHS function
varf opRHS([F, G, H, Pr], [X1, X2, X3, X4]) = int1d(Th)(y^2FX1 + y^2GX2 + y^2HX3)
+on(1, F=0, G=0)
+on(2, F=0, G=0, H=0);

[maximefio (Maxime)]

Maybe I can provide in better way the system of equations :

[maximefio (Maxime)]

I thank you for your help.

Maxime

[cmd (Chris)]

While you can achieve what you want using Lagrange multipliers, that is not necessary.

I assume these BC came from the older works of Khorrami from the late 80’s and 90’s. Note that u+mv=0 and 2u'+mv'=0 necessarily implies u'=v'=0, which is a much simpler requirement. You can impose this naturally as a Neumann BC, and, as far as I am aware, this constraint is more popular nowadays for stability analysis in cylindrical coordinates.

[maximefio (Maxime)]

Hello Chris,

I thank you very much for your answer. This study is indeed based on the works from Khorrami as you stated. I have tried what you suggested by imposing u’=v’=0 for the case Re=2200, m=1 and alpha = 1 and I get the following most unstable mode : re(omega)=0.897, im(omega)=-0.078. The values that obtained Khorrami for these same parameters are re(omega)=0.89663 and im(omega)=−0.048114. So the real part of the most unstable mode seems to be matching but the amplification seems different. Would you know if this is what is expected ? (The modes match well both for the real and imaginary part with Khorrami for m=0 and m>1.) Also, more generally, I would be curious to know how to impose such kind of boundary conditions that are linear conbinations of the perturbations or spatial derivatives, would you have an example of the use of the Lagrange multipliers ? I thank you very much for you help.

[cmd (Chris)]

I would certainly expect your results to agree with Khorrami’s. Would you mind sharing your code as an attachment?

You can see an example with Lagrange multipliers here or, if using the PETSc interface, here.

[maximefio (Maxime)]

Yes, for sure, the code is the following (I an sorry it seems that I cannot deposit attachements)

real ymin = 0.;
real ymax = 1.;
border yline(t=ymin,ymax){x=0.;y=t;}
int np = 2000;

meshL Th=buildmeshL(yline(np));

fespace Xh(Th, P2);
fespace Mh(Th, P1);
// interpolate implicitely on P2 space
Xh W=1-y^2, V=0, U=0, Re=2200.0, Visco=1/Re, PrMean = 0.;
Xh alpha=1.0, m=1.0;

fespace Vh(Th, [P2, P2, P2, P1]);
Vh<complex> [F, G, H, Pr], [X1, X2, X3, X4];

varf opLHS([F, G, H, Pr], [X1, X2, X3, X4]) =
int1d(Th)( -(1i*y^2/Re)*dy(F)*dy(X1)-((1i/Re)*((m^2+1)+y^2*alpha^2))*F*X1+y^2*alpha*W*F*X1-(1i)*(y/Re)*F*dy(X1)-(2*1i*m/Re)*G*X1+y^2*Pr*dy(X1)
           +(1i)*(-(y^2/Re)*dy(G)*dy(X2))-((1i)/Re)*((m^2+1)+y^2*alpha^2)*G*X2+y^2*alpha*W*G*X2-(1i)*(y/Re)*G*dy(X2)-(2*m*(1i)/Re)*F*X2+m*y*Pr*X2
           -(1i)*(y^2/Re)*dy(H)*dy(X3)-(1i)*(y/Re)*H*dy(X3)+y^2*alpha*W*H*X3-(1i/Re)*(m^2+alpha^2*y^2)*H*X3+y^2*dy(W)*F*X3+1i*y^2*alpha*Pr*X3
           -X4*(F+m*G+alpha*y*H)+y*F*dy(X4)
           )
           +on(1, H=0, Pr=0)
           +on(2, F=0, G=0, H=0);


varf opRHS([F, G, H, Pr], [X1, X2, X3, X4]) = int1d(Th)(y^2*F*X1 + y^2*G*X2 + y^2*H*X3)
                                              +on(1, H=0, Pr=0)
                                              +on(2, F=0, G=0, H=0);



// assemble matrices
cout << "Assemble matrices" << endl;
cout << "Matrix L"<< endl;
matrix<complex> L = opLHS(Vh, Vh, solver=GMRES, tgv=1.0e30);
cout << "Matrix B"<< endl;
matrix<complex> B = opRHS(Vh, Vh, solver=GMRES, eps=1.0e-20, tgv=1.0e30);

// save all the information
load "bfstream"
cout << "Write files "<< endl;
complex[int] C(1);
int[int] I(1),J(1);
[I,J,C]=L;
cout << "Matrix L"<< endl;
{
	ofstream fid("LIN/LNS-I.dat");
	fid.write(I);
}
{
	ofstream fid("LIN/LNS-J.dat");
	fid.write(J);
}
{
	ofstream fid("LIN/LNS-C.dat");
	fid.write(C);
}
[I,J,C]=B;
cout << "Matrix B"<< endl;
{
	ofstream fid("LIN/B-I.dat");
	fid.write(I);
}
{
	ofstream fid("LIN/B-J.dat");
	fid.write(J);
}
{
	ofstream fid("LIN/B-C.dat");
	fid.write(C);
}

// get some information about the state vector
Vh [iu1,iu2,iu3,iu4] = [0,1,2,3];
Vh [u1x,u2x,u3x,u4x] = [x,x,x,x];
Vh [u1y,u2y,u3y,u4y] = [y,y,y,y];
Vh [u1b,u2b,u3b,u4b] = [U,V,W,PrMean];

{ofstream fid("LIN/dofs.dat");
	for (int j=0; j<iu1[].n; j++)
		fid << iu1[][j] << " "<< u1x[][j] << " " << u1y[][j] << endl;
}

{ofstream fid("LIN/base.dat");
	for (int j=0; j<iu1[].n; j++)
		fid << u1b[][j] << endl;
}

I use Freefem to generate the matrices A and B and then I read the matrices with python and solve the eigenvalue problem with slepc5py that calls slepc. Would you want also this python script ?
I thank you once again for your help.

[cmd (Chris)]

Please edit that and place the code inside

a `Preformatted text` block

so that I can copy/paste.

[maximefio (Maxime)]

Yes, my apologies, I didn’t know that. This is done normally.

[cmd (Chris)]

It looks like you are imposing boundary conditions incorrectly. Your Dirichlet BC’s should depend on m and you should not enforce any BC on the pressure since you use a Taylor-Hood type FE space.

[maximefio (Maxime)]

Ok, thank you then maybe I understood badly what you meant. I Thought that the two conditions u+mv=0 and 2u'+mv'=0 necessarily implied u'=v'=0, which in the finite element approach are imposed by default. But maybe I also understood something wrong for the pressure, do I need to add something for the pressure disturbance in terms of boundary condition ? (Sorry if these are obvious questions).

[cmd (Chris)]

No you do not need to impose anything special for the pressure, simply eliminate the pressure from the on(). I would suggest something like the following:

Xh W=1-y^2, V=0, U=0;
real Re = 2200.0, alpha=1.0; // these dont need to be defined on an FE space
int m = 1; //m is an integer (makes logical statements easier)

then

varf opLHS([F, G, H, Pr], [X1, X2, X3, X4])
  = int1d(Th)(
    ...
  )
  + on((abs(m) != 1)*1, F = 0, G = 0)
  + on((abs(m) >  0)*1, H = 0)
  + on(2, F=0, G=0, H=0);

and, of course,

varf opRHS([F, G, H, Pr], [X1, X2, X3, X4])
  = int1d(Th)(
    ...
  )
  + on((abs(m) != 1)*1, F = 0, G = 0)
  + on((abs(m) >  0)*1, H = 0)
  + on(2, F=0, G=0, H=0);

Also, the implemented system of equations you’ve shared is not the same as equations 7-10 above, as you must have derived the corresponding weak form. I would make sure you are certain that you didn’t mix anything up during that derivation, as it is a bit nuanced in cylindrical coordinates. I haven’t caught or checked for any errors, just thought I’d mention that. You should derive the weak form for the general Navier Stokes equations, and then apply the parallel flow assumptions to the resulting weak form.

[maximefio (Maxime)]

hello Chris,

I thank you very much for the code, this is indeed working now, I recover the proper real and imaginary part for the most unstable mode as the one from Khorrami.

I have just one question, for the mode m=±1, normally along the axis we should impose w=0,p=0, u+mv=0 and 2u’+mv’=0, or in the code that you provided, the only condition imposed is on u=0 and v=0, could you explain me why you are doing so ?

Also, regarding the equation (7)-(10), to make them in the weak form I have transform all laplacian(perturbation) x func_test by -grad(perturbation) x grad(func_test) and grad(perturbation) x func_test by -perturbation x grad(func_test), is it what should be done ? And finally does it makes a difference to do this process of setting it in weak form the final form of the equations (i.e after linearization of the Navier-Stokes equations + perturbation decomposed as normal modes), or should this process be done right after the linearization ?

I thank you once again for your answers, it is of great help for me.

Maxime

[cmd (Chris)]

Hi Maxime,

The pressure BC in Khorrami’s earlier work is artificial. You can read more about that here.

It sounds like you have done it correctly if your results match. One thing to be wary about when integrating by parts is that the complex inner product involves a conjugation, so the Fourier derivatives of the test function use the conjugate of the wavenumber. In other words, the sign of m and alpha in the test function must be flipped. This should have no effect on flows that are in the O(2) symmetry group because the spectrum is symmetric with +/- m, but if you break the O(2) symmetry e.g. via swirl this is no longer the case.

Best,
Chris

[Amitkaush]

@cmd

and, of course,

hi Maxime,
I am also working on hydrodynamic stability in cartesian coordinates. Can you share the python code to calculate eigenvector and eigenvalues.

Thanks

[Amitkaush]

hi Maxime,
I am also working on hydrodynamic stability in cartesian coordinates. Can you share the python code to calculate eigenvector and eigenvalues.

Thanks