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)(
// gradient pressure
//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"
load "iovtk"
load "MUMPS"
//load "PETSc"
// Mesh
real xmin = -1;
real xmax = 1;
real ymin = -1;
real ymax = 1;
int hot = 1;
int cold = 2;
int adiab = 3;
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)(
// gradient pressure
//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)(
// gradient pressure
//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);