Boundary condition in eigenvalue problem

General Discussion
1

Hello,

I am currently trying to do the linear stability around a sphere flow at low Reynolds number (following the work of P. Meliga, J.-M.Chomaz and D.Sippa ( Unsteadiness in the wake of disks and spheres: Instability, receptivity and control using direct and adjoint global stability analyses). I have implemented the incompressible linearized Navier-Stokes equations in cylindrical coordinates and provided the base flow based on a simulations from Freefem++. My main question concerns the boundary conditions to be applied. These are the following one based on the domain notations in attached:

domain
domain807×332 18.9 KB

boundary_cond

The way to impose Dirichlet and Neumann boundary conditions seems relatively clear for me but I have more difficulites to understand how to impose the no-stress boundary condition on the outlet boundary.

Would it be possible to head me to some existing scripts or page that could help to impose such kind of boundary condition.

I Thank you very much,

Maxime Fiore

2

Hello Maxime,
Usually the variational formulation for the no-stress condition is obtained by taking
any test function without any boundary condition nor boundary integral.
For the free slip condition one has to impose
\vec u \cdot n = 0
and take any test function satisfying the same normality condition (still without any boundary integral).
I think that the formal writing of the free slip condition should include
\partial_r v - v/r instead of \partial_r v alone, because of the free rotating solution v=r that should be allowed.
François.

3

You can see some example implementations here. If you integrate the weak pressure gradient and viscous terms by parts, you simply neglect writing the boundary integral to enforce the free outflow condition. Note that this BC is NOT a no-stress condition, even though it is often incorrectly called that in the literature.

4

Hello François and Chris,

I thank you very much for your answers and the explanations, it seems clear to me now. I will have use the example proposed to implement it.

Thank you once again,

Maxime

5

Long Time ago I have write slide on this. Subject

The example are in the url directory

6

Hello,

I thank you for the document.

7

You may look at this post for the implementation of the slip condition on a curved boundary.

8

Hello François, I thank you for your answer, so if I understand well, the difficulty is not really on the outlet boundary condition but rather on inlet and sphere boundaries. Would it be possible to share with you my code because for inlet and sphere boudaries, I tought impose zero velosity on these conditions was enough but it doesn’t seem to be the case if I understand well:

include "getARGV.idp"

// read the parameters
real m = 1; // helical mode m=1
// data from a laminar simulation using Freefem++
string meshname = getARGV("-mesh", "BASE/ffem-mesh.msh");
string urname = getARGV("-ur", "BASE/ffem-ur.dat"); //ur
string utname = getARGV("-ut", "BASE/ffem-ut.dat"); // utheta
string uxname = getARGV("-ux", "BASE/ffem-ux.dat"); // ux
string nuname = getARGV("-nu", "BASE/ffem-nu.dat"); //viscosity

// import  mesh
mesh Th = readmesh(meshname);

// prepare the vector spaces
fespace Xh(Th, P2);
fespace Mh(Th, P1);

// import meanfields on P1 space
Mh ur, ut, ux, nu;
{
        ifstream fid(urname);
        fid >> ur[];
}
{
        ifstream fid(utname);
        fid >> ut[];
}
{
        ifstream fid(uxname);
        fid >> ux[];
}
{
        ifstream fid(nuname);
        fid >> nu[];
}
// interpolate implicitely on P2 space
Xh U=ur, V=ut, W=ux, Visco=nu;
Mh PrMean;

Xh dUdr=dr(U), dVdr=dr(V), dWdr=dr(W);
Xh dUdx=dx(U), dVdx=dx(V), dWdx=dx(W);


fespace Vh(Th, [P2, P2, P2, P1]);
Vh<complex> [F, G, H, Pr], [X1, X2, X3, X4]; // [ur_fluct,utheta_fluct,ux_fluct,p_fluct]
Xh Re = 1.0/Visco;

//  LHS matrix function (sorted by independent to m, linear to m and quadratic to m), y represents r, momentum equations are multiplied by i * y^2 to prevent issues on the axis and have directly the problem LHS * X= omega * RHS * X
varf opLHS([F, G, H, Pr], [X1, X2, X3, X4]) =
int2d(Th)(1i*(  - y^2*X1*(W*dx(F) + H*dUdx + U*dy(F) + F*dUdr) + y*X1*2*G*V
                            - y^2*X2*(W*dx(G) + H*dVdx + U*dy(G) + F*dVdr) - y*X2*U*G - y*X2*V*F
                            - y^2*X3*(W*dx(H) + H*dWdx + U*dy(H) + F*dWdr)
                            - 1/Re*(  y^2*dx(H)*dx(X3) + y^2*dy(H)*dy(X3)
                                            + y^2*dx(G)*dx(X2) + y^2*dy(G)*dy(X2) + G*X2
                                            + y^2*dx(F)*dx(X1) + y^2*dy(F)*dy(X1) + F*X1  )
                            + Pr*(y^2*dx(X3) + y*X1 + y^2*dy(X1))  )
                  -X4*(y*dx(H) + F + y*dy(F))
+1i*(-y*X1*1i*m*V*F
                          -y*X2*1i*m*V*G
                          -y*X3*1i*m*V*H
                          -1/Re*(+2i*m*X1*G - 2i*m*X2*F)
                          +Pr*y*(-1i*m*X2))
                  -X4*(+1i*m*G)
+1i*(-1/Re*(m^2*H*X3 + m^2*G*X2 + m^2*F*X1))
                 )

+on(1, F=0.0, G=0.0, H=0.0)
+on(2, H=0.0, Pr=0.0)
+on(3, F=0.0)
+on(4, F=0.0, G=0.0, H=0.0)
;

//  RHS matrix function
varf opRHS([F, G, H, Pr], [X1, X2, X3, X4]) =
int2d(Th)(y^2*(F*X1 + G*X2 + H*X3))
+on(1, F=0.0, G=0.0, H=0.0)
+on(2, H=0.0, Pr=0.0)
+on(3, F=0.0)
+on(4, F=0.0, G=0.0, H=0.0)
                ;

// The boundaries are defined this way:
//  o—————————3—————————  o
//  |                                                              |
//  4                                                            5
//  |            ╭—1—╮                                 |
//  o-----2-----o  0  o————2 —————    o
// 1: sphere, 2: axis, 3: top condition, 4: inlet, 5: outlet

I thank you for your kind help.

9

Did you use tgv=-1 or other for the matrix LHS and RHS ?
Maybe you should not put boundary condition in the RHS matrix if you use penalty to impose the boundary condition (ex. tgv =1e30) .
When you use EPSsolve to solve the eigenvalue problem, the shift-invert methods might be done by the EPSsolve and RHS^{-1} or something like that will be involved.

You could refer to the “note” in the document:

https://doc.freefem.org/models/eigen-value-problems.html

10

Besides, the outlet boundary condition might be called " free outflow" more properly. It is not the strict “no stress”. But for the incompressible NS equations, the results should be almost the same if the viscous terms (with “grad of velocity” or “rate of deformation tensor”) in your variational form is written correctly and computational domain is long enough.

11

Maxime,
My message was rather about the external boundary where you impose a slip condition. If I understand well this boundary is curved (in the sense that the normal is not constant on this boundary).
On the boundaries where you impose homogeneous Dirichlet I see no particular problem.
However you have to be sure (from you modeling situation) of the boundary conditions that you want.

12

Hello,

Thank you for your anwers. For the computational domain, I use the one I presented in my first message so the domain is simply a rectangular box except at the sphere location. So I was wondering whether the use of slip condition over a curved boundary implementation was better for the present case.

Also, for the solving of the eigenvalue problem, I output the left hand side and right hand side matrix and solve the eigenvalue problem with PETSc and SLEPc that use a shift-invert method. So I tought that apply the boundary conditions in Freefem++ on both left and right hand side matrices, output these matrices, and solve the problem with SLEPc would be correct, but it is not te case based on the previous messages ?

13

Ok it was not clear to me if the domain was 2d (rectangular) or 3d (cylinder).
For 2d rectangular the boundary is flat thus for the slip condition you can use any penalty method.

14

You should not use tgv=-1 on both the LHS and RHS. You should use -1 on the LHS, and -10 on the RHS.

15

Thank you for the answers, may I ask why different values of tgv are used for left and right hand side matrices. Also I am wondering what is the impact of setting such kind of values on the matrices themselves ?

Thank you.

My apologies, I have read this topic: Application of Dirichlet BCs
and now it seems clearer to me why different values of tgv are used.