A Large Fluid Problem
A friend of one of us in Auroville-India was building a ramp to access an air conditioned room. As I was visiting the construction site he told me that he expected to cool air escaping by the door to the room to slide down the ramp and refrigerate the feet of the coming visitors. I told him “no way” and decided to check numerically.
The fluid velocity and pressure are solution of the Navier-Stokes equations with varying density function of the temperature.
The geometry is trapezoidal with prescribed inflow made of cool air at the bottom and warm air above and so are the initial conditions; there is free outflow, slip velocity at the top (artificial) boundary and no-slip at the bottom. However the Navier-Stokes cum temperature equations have a RANS \(k-\epsilon\) model and a Boussinesq approximation for the buoyancy. This comes to :
We use a time discretization which preserves positivity and uses the method of characteristics (\(X^m(x)\approx x-u^m(x)\delta t\))
In variational form and with appropriated boundary conditions the problem is :
1 load "iovtk"
2
3 verbosity=0;
4
5 // Parameters
6 int nn = 15;
7 int nnPlus = 5;
8 real l = 1.;
9 real L = 15.;
10 real hSlope = 0.1;
11 real H = 6.;
12 real h = 0.5;
13
14 real reylnods =500;
15 real beta = 0.01;
16
17 real eps = 9.81/303.;
18 real nu = 1;
19 real numu = nu/sqrt(0.09);
20 real nuep = pow(nu,1.5)/4.1;
21 real dt = 0.;
22
23 real Penalty = 1.e-6;
24
25 // Mesh
26 border b1(t=0, l){x=t; y=0;}
27 border b2(t=0, L-l){x=1.+t; y=-hSlope*t;}
28 border b3(t=-hSlope*(L-l), H){x=L; y=t;}
29 border b4(t=L, 0){x=t; y=H;}
30 border b5(t=H, h){x=0; y=t;}
31 border b6(t=h, 0){x=0; y=t;}
32
33 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));
34 plot(Th);
35
36 // Fespaces
37 fespace Vh2(Th, P1b);
38 Vh2 Ux, Uy;
39 Vh2 Vx, Vy;
40 Vh2 Upx, Upy;
41
42 fespace Vh(Th,P1);
43 Vh p=0, q;
44 Vh Tp, T=35;
45 Vh k=0.0001, kp=k;
46 Vh ep=0.0001, epp=ep;
47
48 fespace V0h(Th,P0);
49 V0h muT=1;
50 V0h prodk, prode;
51 Vh kappa=0.25e-4, stress;
52
53 // Macro
54 macro grad(u) [dx(u), dy(u)] //
55 macro Grad(U) [grad(U#x), grad(U#y)] //
56 macro Div(U) (dx(U#x) + dy(U#y)) //
57
58 // Functions
59 func g = (x) * (1-x) * 4;
60
61 // Problem
62 real alpha = 0.;
63
64 problem Temperature(T, q)
65 = int2d(Th)(
66 alpha * T * q
67 + kappa* grad(T)' * grad(q)
68 )
69 + int2d(Th)(
70 - alpha*convect([Upx, Upy], -dt, Tp)*q
71 )
72 + on(b6, T=25)
73 + on(b1, b2, T=30)
74 ;
75
76 problem KineticTurbulence(k, q)
77 = int2d(Th)(
78 (epp/kp + alpha) * k * q
79 + muT* grad(k)' * grad(q)
80 )
81 + int2d(Th)(
82 prodk * q
83 - alpha*convect([Upx, Upy], -dt, kp)*q
84 )
85 + on(b5, b6, k=0.00001)
86 + on(b1, b2, k=beta*numu*stress)
87 ;
88
89 problem ViscosityTurbulence(ep, q)
90 = int2d(Th)(
91 (1.92*epp/kp + alpha) * ep * q
92 + muT * grad(ep)' * grad(q)
93 )
94 + int1d(Th, b1, b2)(
95 T * q * 0.001
96 )
97 + int2d(Th)(
98 prode * q
99 - alpha*convect([Upx, Upy], -dt, epp)*q
100 )
101 + on(b5, b6, ep=0.00001)
102 + on(b1, b2, ep=beta*nuep*pow(stress,1.5))
103 ;
104
105 // Initialization with stationary solution
106 solve NavierStokes ([Ux, Uy, p], [Vx, Vy, q])
107 = int2d(Th)(
108 alpha * [Ux, Uy]' * [Vx, Vy]
109 + muT * (Grad(U) : Grad(V))
110 + p * q * Penalty
111 - p * Div(V)
112 - Div(U) * q
113 )
114 + int1d(Th, b1, b2, b4)(
115 Ux * Vx * 0.1
116 )
117 + int2d(Th)(
118 eps * (T-35) * Vx
119 - alpha*convect([Upx, Upy], -dt, Upx)*Vx
120 - alpha*convect([Upx, Upy], -dt, Upy)*Vy
121 )
122 + on(b6, Ux=3, Uy=0)
123 + on(b5, Ux=0, Uy=0)
124 + on(b1, b4, Uy=0)
125 + on(b2, Uy=-Upx*N.x/N.y)
126 + on(b3, Uy=0)
127 ;
128
129 plot([Ux, Uy], p, value=true, coef=0.2, cmm="[Ux, Uy] - p");
130
131 {
132 real[int] xx(21), yy(21), pp(21);
133 for (int i = 0 ; i < 21; i++){
134 yy[i] = i/20.;
135 xx[i] = Ux(0.5,i/20.);
136 pp[i] = p(i/20.,0.999);
137 }
138 cout << " " << yy << endl;
139 plot([xx, yy], wait=true, cmm="Ux x=0.5 cup");
140 plot([yy, pp], wait=true, cmm="p y=0.999 cup");
141 }
142
143 // Initialization
144 dt = 0.1; //probably too big
145 int nbiter = 3;
146 real coefdt = 0.25^(1./nbiter);
147 real coefcut = 0.25^(1./nbiter);
148 real cut = 0.01;
149 real tol = 0.5;
150 real coeftol = 0.5^(1./nbiter);
151 nu = 1./reylnods;
152
153 T = T - 10*((x<1)*(y<0.5) + (x>=1)*(y+0.1*(x-1)<0.5));
154
155 // Convergence loop
156 real T0 = clock();
157 for (int iter = 1; iter <= nbiter; iter++){
158 cout << "Iteration " << iter << " - dt = " << dt << endl;
159 alpha = 1/dt;
160
161 // Time loop
162 real t = 0.;
163 for (int i = 0; i <= 500; i++){
164 t += dt;
165 cout << "Time step " << i << " - t = " << t << endl;
166
167 // Update
168 Upx = Ux;
169 Upy = Uy;
170 kp = k;
171 epp = ep;
172 Tp = max(T, 25); //for beauty only should be removed
173 Tp = min(Tp, 35); //for security only should be removed
174 kp = max(k, 0.0001); epp = max(ep, 0.0001); // to be secure: should not be active
175 muT = 0.09*kp*kp/epp;
176
177 // Solve NS
178 NavierStokes;
179
180 // Update
181 prode = -0.126*kp*(pow(2*dx(Ux),2)+pow(2*dy(Uy),2)+2*pow(dx(Uy)+dy(Ux),2))/2;
182 prodk = -prode*kp/epp*0.09/0.126;
183 kappa = muT/0.41;
184 stress = abs(dy(Ux));
185
186 // Solve k-eps-T
187 KineticTurbulence;
188 ViscosityTurbulence;
189 Temperature;
190
191 // Plot
192 plot(T, value=true, fill=true);
193 plot([Ux, Uy], p, coef=0.2, cmm=" [Ux, Uy] - p", WindowIndex=1);
194
195 // Time
196 cout << "\tTime = " << clock()-T0 << endl;
197 }
198
199 // Check
200 if (iter >= nbiter) break;
201
202 // Adaptmesh
203 Th = adaptmesh(Th, [dx(Ux), dy(Ux), dx(Ux), dy(Uy)], splitpbedge=1, abserror=0, cutoff=cut, err=tol, inquire=0, ratio=1.5, hmin=1./1000);
204 plot(Th);
205
206 // Update
207 dt = dt * coefdt;
208 tol = tol * coeftol;
209 cut = cut * coefcut;
210 }
211 cout << "Total Time = " << clock()-T0 << endl;