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 :

  1load "iovtk"
2
3verbosity=0;
4
5// Parameters
6int nn = 15;
7int nnPlus = 5;
8real l = 1.;
9real L = 15.;
10real hSlope = 0.1;
11real H = 6.;
12real h = 0.5;
13
14real reylnods =500;
15real beta = 0.01;
16
17real eps = 9.81/303.;
18real nu = 1;
19real numu = nu/sqrt(0.09);
20real nuep = pow(nu,1.5)/4.1;
21real dt = 0.;
22
23real Penalty = 1.e-6;
24
25// Mesh
26border b1(t=0, l){x=t; y=0;}
27border b2(t=0, L-l){x=1.+t; y=-hSlope*t;}
28border b3(t=-hSlope*(L-l), H){x=L; y=t;}
29border b4(t=L, 0){x=t; y=H;}
30border b5(t=H, h){x=0; y=t;}
31border b6(t=h, 0){x=0; y=t;}
32
33mesh 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));
34plot(Th);
35
36// Fespaces
37fespace Vh2(Th, P1b);
38Vh2 Ux, Uy;
39Vh2 Vx, Vy;
40Vh2 Upx, Upy;
41
42fespace Vh(Th,P1);
43Vh p=0, q;
44Vh Tp, T=35;
45Vh k=0.0001, kp=k;
46Vh ep=0.0001, epp=ep;
47
48fespace V0h(Th,P0);
49V0h muT=1;
50V0h prodk, prode;
51Vh kappa=0.25e-4, stress;
52
53// Macro
56macro Div(U) (dx(U#x) + dy(U#y)) //
57
58// Functions
59func g = (x) * (1-x) * 4;
60
61// Problem
62real alpha = 0.;
63
64problem Temperature(T, q)
65   = int2d(Th)(
66        alpha * T * 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
76problem KineticTurbulence(k, q)
77   = int2d(Th)(
78        (epp/kp + alpha) * k * 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
89problem ViscosityTurbulence(ep, q)
90   = int2d(Th)(
91        (1.92*epp/kp + alpha) * ep * 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
106solve NavierStokes ([Ux, Uy, p], [Vx, Vy, q])
107   = int2d(Th)(
108        alpha * [Ux, Uy]' * [Vx, Vy]
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
129plot([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
144dt = 0.1; //probably too big
145int nbiter = 3;
146real coefdt = 0.25^(1./nbiter);
147real coefcut = 0.25^(1./nbiter);
148real cut = 0.01;
149real tol = 0.5;
150real coeftol = 0.5^(1./nbiter);
151nu = 1./reylnods;
152
153T = T - 10*((x<1)*(y<0.5) + (x>=1)*(y+0.1*(x-1)<0.5));
154
155// Convergence loop
156real T0 = clock();
157for (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
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}
211cout << "Total Time = " << clock()-T0 << endl;


A large fluid problem

Table of content