# 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 :

\begin{eqnarray} \p_t\theta+u\n\theta-\n\cdot(\kappa_T^m\n\theta) &=& 0\\ \p_t u +u\n u -\n\cdot(\mu_T\n u) +\n p+ e(\theta-\theta_0)\vec e_2 &=&0\\ \n\cdot u &=& 0\\ \mu_T &=& c_\mu\frac{k^2}\epsilon\\ \kappa_T &=& \kappa\mu_T\\ \p_t k + u\n k + \epsilon -\n\cdot(\mu_T\n k) &=& \frac{\mu_T}2|\n u+\n u^T|^2\\ \p_t\epsilon+u\n\epsilon + c_2\frac{\epsilon^2} k -\frac{c_\epsilon}{c_\mu}\n\cdot (\mu_T\n\epsilon) &=& \frac{c_1}2 k|\n u+\n u^T|^2\\ \end{eqnarray}

We use a time discretization which preserves positivity and uses the method of characteristics ($X^m(x)\approx x-u^m(x)\delta t$)

\begin{eqnarray} \frac 1{\delta t}(\theta^{m+1}-\theta^m \circ X^m)-\n\cdot(\kappa_T^m\n\theta^{m+1}) &=& 0\\ \frac1{\delta t}(u^{m+1}-u^m \circ X^m) -\n\cdot(\mu_T^m\n u^{m+1}) +\n p^{m+1}+ e(\theta^{m+1}-\theta_0)\vec e_2 &=& 0\\ \n\cdot u^{m+1} &=& 0\\ \frac1{\delta t}(k^{m+1}-k^m \circ X^m) + k^{m+1}\frac{\epsilon^m}{k^m} -\n\cdot(\mu_T^m\n k^{m+1}) &=& \frac{\mu_T^m}2|\n u^m+{\n u^m}^T|^2\\ \frac1{\delta t}(\epsilon^{m+1}-\epsilon^m \circ X^m) + c_2\epsilon^{m+1}\frac{\epsilon^m} {k^m} -\frac{c_\epsilon}{c_\mu}\n\dot(\mu_T^m\n\epsilon^{m+1}) &=& \frac{c_1}2 k^m|\n u^m+{\n u^m}^T|^2\\ \mu_T ^{m+1} &=& c_\mu\frac{{k^{m+1}}^2}{\epsilon^{m+1}}\\ \kappa_T^{m+1} &=& \kappa\mu_T^{m+1} \end{eqnarray}

In variational form and with appropriated boundary conditions the problem is :

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 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 reylnods =500; 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.0001, kp=k; Vh ep=0.0001, 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 Grad(U) [grad(U#x), grad(U#y)] // macro Div(U) (dx(U#x) + dy(U#y)) // // Functions func g = (x) * (1-x) * 4; // Problem real alpha = 0.; problem Temperature(T, q) = int2d(Th)( alpha * T * q + kappa* grad(T)' * grad(q) ) + int2d(Th)( - alpha*convect([Upx, Upy], -dt, Tp)*q ) + on(b6, T=25) + on(b1, b2, T=30) ; problem KineticTurbulence(k, q) = int2d(Th)( (epp/kp + alpha) * k * q + muT* grad(k)' * grad(q) ) + int2d(Th)( prodk * q - alpha*convect([Upx, Upy], -dt, kp)*q ) + on(b5, b6, k=0.00001) + on(b1, b2, k=beta*numu*stress) ; problem ViscosityTurbulence(ep, q) = int2d(Th)( (1.92*epp/kp + alpha) * ep * q + muT * grad(ep)' * grad(q) ) + int1d(Th, b1, b2)( T * q * 0.001 ) + int2d(Th)( prode * q - alpha*convect([Upx, Upy], -dt, epp)*q ) + on(b5, b6, ep=0.00001) + on(b1, b2, ep=beta*nuep*pow(stress,1.5)) ; // Initialization with stationary solution solve NavierStokes ([Ux, Uy, p], [Vx, Vy, q]) = int2d(Th)( alpha * [Ux, Uy]' * [Vx, Vy] + muT * (Grad(U) : Grad(V)) + p * q * Penalty - p * Div(V) - Div(U) * q ) + int1d(Th, b1, b2, b4)( Ux * Vx * 0.1 ) + int2d(Th)( eps * (T-35) * Vx - alpha*convect([Upx, Upy], -dt, Upx)*Vx - alpha*convect([Upx, Upy], -dt, Upy)*Vy ) + on(b6, Ux=3, Uy=0) + on(b5, Ux=0, Uy=0) + on(b1, b4, Uy=0) + on(b2, Uy=-Upx*N.x/N.y) + on(b3, Uy=0) ; 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=true, cmm="Ux x=0.5 cup"); plot([yy, pp], wait=true, cmm="p y=0.999 cup"); } // Initialization dt = 0.1; //probably too big int nbiter = 3; 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); nu = 1./reylnods; 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.0001); epp = max(ep, 0.0001); // to be secure: should not be active muT = 0.09*kp*kp/epp; // Solve NS NavierStokes; // Update prode = -0.126*kp*(pow(2*dx(Ux),2)+pow(2*dy(Uy),2)+2*pow(dx(Uy)+dy(Ux),2))/2; prodk = -prode*kp/epp*0.09/0.126; kappa = muT/0.41; stress = abs(dy(Ux)); // 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; // Adaptmesh 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); plot(Th); // Update dt = dt * coefdt; tol = tol * coeftol; cut = cut * coefcut; } cout << "Total Time = " << clock()-T0 << endl; 
Fig. 1: Temperature at time steps 100, 200, 300, 400, 500. Fig. 2: Velocity at time steps 100, 200, 300, 400, 500.