First guess for the Navier-Stokes and heat equations

I want to solve the incompressible Navier-Stokes equations (NVSE) coupled with a heat equation (convection and diffusion term). The variables are the velocity vector (u, v), the pressure p and the temperature T. The temperature appears in the NVSE with a buoyancy term. I want to solve for the stationary solution by using a Newton’s algorithm. Hence, I want to have a pretty good first guess before using the Newton’s algorithm. The procedure for making the first guess is:

• step 1: starting from zero initialisation of the variables
• step 2: solve for u, v and p by solving the NVSE (without the convection term so that the NVSE is linear)
• step 3: solve for T with the heat equation
• repeat steps 2-3 until convergence.

Do not hesitate if you need the set of equations.

But I get the following error:

``````   64 : 		// buoyancy term
65 : 		+ bulk*T*(gx*uh + gy*vh)
66 : 		// penalization of the continuity equation
67 : 		- error operator +  <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E>, <10LinearCombI6MDroit4C_F0E>
...
Error line number 67, in file heated_cavity.edp, before  token -

current line = 65
Compile error :
line number :67, -
error Compile error :
line number :67, -
code = 1 mpirank: 0
close  MUMPS_SEQ: MPI_Finalize
``````

Here are the definitions of the NVSE and heat equation:

``````problem NVSGuess(u, v, p, uh, vh, q, solver=CG, tgv=-1) =
int2d(Oh)(
//1/rho * (dx(p) * uh + dy(p) * vh)
- 1/rho * p * (dx(uh) + dy(vh))
// laplacian
+ nu * (dx(u)*dx(uh) + dy(u)*dy(uh) + dx(v)*dx(vh) + dy(v)*dy(vh))
// buoyancy term
+ bulk*T*(gx*uh + gy*vh)
// penalization of the continuity equation
- q * (dx(u)+dy(v))
)
- int2d(Oh)(
(gx*uh+gy*vh)
)
// velocity/no slip (u,v)
+ on(adiab, hot, cold, u = 0, v = 0);
//+ on(hot, T=Thot)
//+ on(cold, T=cold);

problem HeatGuess(T, w, solver=sparsesolver) =
// Heat equation
int2d(Oh)(
cp * w * (u*dx(T) + v*dy(T))
+ cond/rho*(dx(T)*dx(w) + dy(T)*dy(w))
)
// dirichlet conditions on T (hot/cold)
+ on(hot, T = Thot)
+ on(cold, T = Tcold);
``````

Here is the full code if needed:

``````include "getARGV.idp"

// Mesh
real xmin = -1;
real xmax =  1;
real ymin = -1;
real ymax =  1;

int hot = 1;
int cold = 2;
border left(t=ymin,ymax){ x=xmin; y=ymax+ymin-t ; label= hot;}
border right(t=ymin,ymax){ x=xmax; y=t; label= cold;}
border bot(t=xmin,xmax){ x=t; y=ymin; label= adiab;}
border top(t=xmin,xmax){ x=xmax+xmin-t; y=ymax; label= adiab;}

int n = getARGV("-n",100);
mesh Oh = buildmesh(
left(n) + bot(n) + right(n) + top(n)
);

//if (mpirank==0){
//	savevtk("initial_mesh.vtu", Oh);
//}
savevtk("initial_mesh.vtu", Oh);

// parameters
real cp = 1.0;
real rho = 1.0;
real cond = getARGV("-cond", 1.0);
real nu = getARGV("-nu", 1.0);
real Thot  = 1.0;
real Tcold = -1.0;
real bulk = 0.1;
real gx = 0;
real gy = -10.0;

// Approximated solution space
fespace Vh2(Oh, P2);
Vh2 u, v;
Vh2 du, dv;
Vh2 uh, vh;

fespace Vh1(Oh, P1);
Vh1 dT, T, w;
Vh1 dp, p, q;

u[] = 0;
v[] = 0;
p[] = 0;
T[] = 0;

int nGuess = getARGV("-nguess", 1);
problem NVSGuess(u, v, p, uh, vh, q, solver=CG, tgv=-1) =
int2d(Oh)(
//1/rho * (dx(p) * uh + dy(p) * vh)
- 1/rho * p * (dx(uh) + dy(vh))
// laplacian
+ nu * (dx(u)*dx(uh) + dy(u)*dy(uh) + dx(v)*dx(vh) + dy(v)*dy(vh))
// buoyancy term
+ bulk*T*(gx*uh + gy*vh)
// penalization of the continuity equation
- q * (dx(u)+dy(v))
)
- int2d(Oh)(
(gx*uh+gy*vh)
)
// velocity/no slip (u,v)
+ on(adiab, hot, cold, u = 0, v = 0);
//+ on(hot, T=Thot)
//+ on(cold, T=cold);

problem HeatGuess(T, w, solver=sparsesolver) =
// Heat equation
int2d(Oh)(
cp * w * (u*dx(T) + v*dy(T))
+ cond/rho*(dx(T)*dx(w) + dy(T)*dy(w))
)
// dirichlet conditions on T (hot/cold)
+ on(hot, T = Thot)
+ on(cold, T = Tcold);

for(int i=0; i<nGuess; i++){
cout << "Guess iteration "<< i+1 <<endl;
/*
solve NVSGuess(u, v, p, uh, vh, q, solver=CG, tgv=-1) =
int2d(Oh)(
//1/rho * (dx(p) * uh + dy(p) * vh)
- 1/rho * p * (dx(uh) + dy(vh))
// laplacian
+ nu * (dx(u)*dx(uh) + dy(u)*dy(uh) + dx(v)*dx(vh) + dy(v)*dy(vh))
// buoyancy term
+ bulk*Told*(gx*uh + gy*vh)
// penalization of the continuity equation
- q * (dx(u)+dy(v))
)
- int2d(Oh)(
(gx*uh+gy*vh)
)
// velocity/no slip (u,v)
+ on(adiab, hot, cold, u = 0, v = 0);
//+ on(hot, T=Thot)
//+ on(cold, T=cold);

solve HeatGuess(T, w, solver=sparsesolver) =
// Heat equation
int2d(Oh)(
cp * w * (u*dx(T) + v*dy(T))
+ cond/rho*(dx(T)*dx(w) + dy(T)*dy(w))
)
// dirichlet conditions on T (hot/cold)
+ on(hot, T = Thot)
+ on(cold, T = Tcold);
*/
NVSGuess;
HeatGuess;
Told[] = T[];
}
int[int] Order = [1,1,1];
savevtk("initial_guess.vtu", Oh, [u,v,0], p, T, dataname="velocity p T", order=Order);
``````

You must write the bilinear terms (which are linear in the unknown) and the linear terms (which do not contain the unknown) in separate int2d().
The buoyancy term is of the second kind, it must be put in another int2d().