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.