load "iovtk"
verbosity=0;
// Parameters
int nn = 15;
int nnPlus = 5;
real l = 1.;
real L = 15.;
real hSlope = 0.1;
real H = 6.;
real h = 0.5;
real beta = 0.01;
real eps = 9.81/303;
real nu = 1.;
real numu = nu/sqrt(0.09);
real nuep = pow(nu,1.5)/4.1;
real dt = 0.;
real Penalty = 1.e-6;
// Mesh
border b1(t=0, l){x=t; y=0;}
border b2(t=0, L-l){x=1.+t; y=-hSlope*t;}
border b3(t=-hSlope*(L-l), H){x=L; y=t;}
border b4(t=L, 0){x=t; y=H;}
border b5(t=H, h){x=0; y=t;}
border b6(t=h, 0){x=0; y=t;}
mesh Th=buildmesh(b1(nnPlus*nn*l) + b2(nn*sqrt((L-l)^2+(hSlope*(L-l))^2)) + b3(nn*(H + hSlope*(L-l))) + b4(nn*L) + b5(nn*(H-h)) + b6(nnPlus*nn*h));
plot(Th);
// Fespaces
fespace Vh2(Th, P1b);
Vh2 Ux, Uy;
Vh2 Vx, Vy;
Vh2 Upx, Upy;
fespace Vh(Th,P1);
Vh p=0, q;
Vh Tp, T=35;
Vh k=0.27, kp=k;
Vh ep=0.002104, epp=ep;
fespace V0h(Th,P0);
V0h muT=1;
V0h prodk, prode;
Vh kappa=0.25e-4, stress;
// Macro
macro grad(u) [dx(u), dy(u)] //
macro norm(u1,u2,p,t,k,ep) (sqrt(u1^2+u2^2+p^2+t^2+k^2+ep^2))//
macro Grad(u1,u2) [dx(u1),dy(u1), dx(u2),dy(u2)] //
macro Gradt(u1,u2) [dx(u1),dx(u2), dy(u1),dy(u2)] //
macro Div(u1,u2) (dx(u1) + dy(u2)) //
macro UgradV(u1,u2,v1,v2) [[u1,u2]'*[dx(v1),dy(v1)],[u1,u2]'*[dx(v2),dy(v2)]]//
macro ugradv(u1,u2,v) [u1,u2]'*grad(v)//
// Problem
real alpha = 0.;
problem Temperature(T, q)
= int2d(Th)(
alpha*T*q
+ kappa*(grad(T)'*grad(q))
)
+ int2d(Th)(
ugradv(Upx, Upy, Tp)
)-int2d(Th)(alpha*Tp*q)
+ on(b6, T=25)
+ on(b1, b2, T=30)
;
problem KineticTurbulence(k, q)
= int2d(Th)(alpha*q*(k-kp)+q*ugradv(Upx,Upy,k)+epp*q)
+int2d(Th)(muT*grad(k)'*grad(q)+prodk * q)
+ on(b5, b6, k=0.27)
;
problem ViscosityTurbulence(ep, q)
= int2d(Th)(
(1.92*epp/kp + alpha) * ep * q-alpha*epp
+ (1./1.3)*muT * grad(ep)' * grad(q)
)
+ int1d(Th, b1, b2)(
-q * (1./1.3)*2.4952*muT*(ep*epp)/((kp)^1.5)
)
+ int2d(Th)(
prode * q
+ugradv(Upx,Upy,ep)
)
+ on(b6, ep=0.002104)
;
// Initialization with stationary solution
solve NavierStokes ([Ux, Uy, p], [Vx, Vy, q])
= int2d(Th)(
alpha * ([Ux, Uy]' * [Vx, Vy]-[Upx, Upy]' * [Vx, Vy])
+ muT * ((Grad(Ux,Uy)+Gradt(Ux,Uy))'*Grad(Vx,Vy))
+ p * q * Penalty
- p * Div(Vx,Vy)
- Div(Ux,Uy) * q
)
+ int1d(Th, b1, b2)(
Ux * Vx * 0.2245662488*kp^0.5
)
+ int2d(Th)(
eps * (T-35) * Vy
+UgradV(Ux,Uy,Upx,Upy)
)
+ on(b6, Ux=3, Uy=0)
+ on(b5, Ux=0, Uy=0)
+ on(b2, Uy=-Upx*N.x/N.y)
;
plot([Ux, Uy], p, value=true, coef=0.2, cmm="[Ux, Uy] - p");
{
real[int] xx(21), yy(21), pp(21);
for (int i = 0 ; i < 21; i++){
yy[i] = i/20.;
xx[i] = Ux(0.5,i/20.);
pp[i] = p(i/20.,0.999);
}
cout << " " << yy << endl;
plot([xx, yy], wait=false, cmm="Ux x=0.5 cup");
plot([yy, pp], wait=false, cmm="p y=0.999 cup");
}
// Initialization
dt = 0.05; //probably too big
int nbiter = 2000;
real coefdt = 0.25^(1./nbiter);
real coefcut = 0.25^(1./nbiter);
real cut = 0.01;
real tol = 0.5;
real coeftol = 0.5^(1./nbiter);
T = T - 10*((x<1)*(y<0.5) + (x>=1)*(y+0.1*(x-1)<0.5));
// Convergence loop
real T0 = clock();
for (int iter = 1; iter <= nbiter; iter++){
cout << "Iteration " << iter << " - dt = " << dt << endl;
alpha = 1/dt;
// Time loop
real t = 0.;
for (int i = 0; i <= 500; i++){
t += dt;
cout << "Time step " << i << " - t = " << t << endl;
// Update
Upx = Ux;
Upy = Uy;
kp = k;
epp = ep;
Tp = max(T, 25); //for beauty only should be removed
Tp = min(Tp, 35); //for security only should be removed
kp = max(k, 0.27); epp = max(ep, 0.002104); // to be secure: should not be active
muT = 0.09*kp*kp/epp;
// Solve NS
NavierStokes;
// Update
prode = -0.0648*kp*(pow(2*dx(Ux),2)+pow(2*dy(Uy),2)+2*pow(dx(Uy)+dy(Ux),2));
prodk = prode*(kp/epp)*0.045/0.0648;
kappa = muT/0.9;
// Solve k-eps-T
KineticTurbulence;
ViscosityTurbulence;
Temperature;
// Plot
plot(T, value=true, fill=true);
plot([Ux, Uy], p, coef=0.2, cmm=" [Ux, Uy] - p", WindowIndex=1);
// Time
cout << "\tTime = " << clock()-T0 << endl;
}
// Check
if (iter >= nbiter) break;
// Update
dt = dt * coefdt;
tol = tol * coeftol;
cut = cut * coefcut;
}
cout << "Total Time = " << clock()-T0 << endl;
when I try to compile this code, it says that
65 : ugradv(Upx, Upy, Tp) [Upx, Upy]'*grad( Tp) [dx( Tp), dy( Tp)]
66 : )- error operator + <12FormBilinear>, <d>
the error occurs at temperature part