FreeFEM Documentation on GitHub

stars - forks

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{split}\begin{array}{rcl} \partial_t\theta+u\nabla\theta-\nabla\cdot(\kappa_T^m\nabla\theta) &=& 0\\ \partial_t u +u\nabla u -\nabla\cdot(\mu_T\nabla u) +\nabla p+ e(\theta-\theta_0)\vec e_2 &=&0\\ \nabla\cdot u &=& 0\\ \mu_T &=& c_\mu\frac{k^2}\epsilon\\ \kappa_T &=& \kappa\mu_T\\ \partial_t k + u\nabla k + \epsilon -\nabla\cdot(\mu_T\nabla k) &=& \frac{\mu_T}2|\nabla u+\nabla u^T|^2\\ \partial_t\epsilon+u\nabla\epsilon + c_2\frac{\epsilon^2} k -\frac{c_\epsilon}{c_\mu}\nabla\cdot (\mu_T\nabla\epsilon) &=& \frac{c_1}2 k|\nabla u+\nabla u^T|^2\\ \end{array}\end{split}\]

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{split}\begin{array}{rcl} \frac 1{\delta t}(\theta^{m+1}-\theta^m \circ X^m)-\nabla\cdot(\kappa_T^m\nabla\theta^{m+1}) &=& 0\\ \frac1{\delta t}(u^{m+1}-u^m \circ X^m) -\nabla\cdot(\mu_T^m\nabla u^{m+1}) +\nabla p^{m+1}+ e(\theta^{m+1}-\theta_0)\vec e_2 &=& 0\\ \nabla\cdot u^{m+1} &=& 0\\ \frac1{\delta t}(k^{m+1}-k^m \circ X^m) + k^{m+1}\frac{\epsilon^m}{k^m} -\nabla\cdot(\mu_T^m\nabla k^{m+1}) &=& \frac{\mu_T^m}2|\nabla u^m+{\nabla 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}\nabla\dot(\mu_T^m\nabla\epsilon^{m+1}) &=& \frac{c_1}2 k^m|\nabla u^m+{\nabla 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{array}\end{split}\]

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;
ALargeFluidProblem_T_100

Fig. 27 Temperature at time step 100

ALargeFluidProblem_U_100

Fig. 28 Velocity at time step 100

ALargeFluidProblem_T_200

Fig. 29 Temperature at time step 200

ALargeFluidProblem_U_200

Fig. 30 Velocity at time step 200

ALargeFluidProblem_T_300

Fig. 31 Temperature at time step 300

ALargeFluidProblem_U_300

Fig. 32 Velocity at time step 300

ALargeFluidProblem_T_400

Fig. 33 Temperature at time step 400

ALargeFluidProblem_U_400

Fig. 34 Velocity at time step 400

ALargeFluidProblem_T_500

Fig. 35 Temperature at time step 500

ALargeFluidProblem_U_500

Fig. 36 Velocity at time step 500

A large fluid problem

Table of content