// Discussion: // // Vortex caused by free convection // // Solves the Laplace Equation on a parallel plate capacitor boundary. // P1 Elements are used. // // Export the data, using ffmatlib, for graphic processing by MATLAB or Octave. // // Licensing: // // Copyright (C) 2018 Chloros2 // // This program is free software: you can redistribute it and/or modify it // under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // This program is distributed in the hope that it will be useful, but // WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program. If not, see // . // // Modified: // // 2018-05-19 // // Author: // // Chloros2 // // Local: // // real dt, the time step. // // real g, gravitational coefficient. // // real Tcold, the temperature of the cold wall, in centigrade. // // real tend, the final time. // // real Thot, the temperature of the hot wall, in centigrade. // include "ffmatlib.idp" load "UMFPACK64"; cout << endl; cout << "convective_rolls:" << endl; cout << " FreeFem++ version." << endl; cout << " Solve the Laplace equation on a parallel plate capacitor boundary." << endl; cout << " Record data about the resulting convective rolls." << endl; cout << " Using the ffmatlib interface, export the data for graphic" << endl; cout << " processing by MATLAB/Octave." << endl; // // Geometry. // real width = 0.30; real height = 0.05; int m = 300; real tend = 5.0; real dt = 0.04; // // Set physical parameters. // real g = 9.81; real Thot = 50.0; real Tcold = 0.0; real beta = (1.0/(273.0+0.5*(Thot+Tcold))); // Expansion coefficient //properties air @25°C / 1 bar real rho = 1.168; // density kg/m^3 real nu = 15.82e-6; // kinematic viscosity m^2/s real cp = 1.007e+3; // specific isobar heat capacity J/kg*K real lambda = 0.0261; // heat conductivity W/mK real mu = nu * rho; real a = lambda / ( rho * cp ); // // Define the boundaries of the region. // int C1=1; int C2=2; int C3=3; int C4=4; int C10=5; border floor(t=0,width){ x=t; y=0; label=C1;}; border ceiling(t=width,0){ x=t; y=height; label=C2;}; border right(t=0,height){ x=width; y=t; label=C3;}; border left(t=height,0){ x=0; y=t; label=C4;}; // // Create the mesh. // mesh Th = buildmesh ( floor ( width * m ) + right ( height * m ) + ceiling ( width * m ) + left ( height * m ) ); // // Plot the mesh. // plot ( Th, wait = false, ps = "convective_rolls_mesh.ps" ); // // Set finite element spaces for velocity and pressure. // fespace Vh ( Th, P2 ); fespace Qh ( Th, P1 ); // //Fluid stuff // Vh u1,u2,v2,v1,up=0,vp=0; Qh p=0, pp; Qh psi,phi; // //r is a temperature rise/difference // Qh rr, rp=0,r=0; Qh Tfluid; //Heatflux density Qh qx, qy; bool reuseMatrix=false; int n=0; // // Define Navier Stokes problem for velocity and pressure. // problem NS([u1,u2,p],[v1,v2,pp],solver=UMFPACK, init=reuseMatrix) = int2d(Th)( u1*v1 + u2*v2 //Remains from advective part in lagrange coordinates + nu*dt*(dx(u1)*dx(v1) + dy(u1)*dy(v1) //Viscous part + dx(u2)*dx(v2) + dy(u2)*dy(v2)) + p*pp*1.e-6 //Stabilization - dt*(p*(dx(v1) + dy(v2))/rho //Pressure part - dt*pp*(dx(u1) + dy(u2)))) //Mass conservation (continuity) - int2d(Th) (convect([up,vp],-dt,up)*v1 //Advection Lagrangian coordinates + convect([up,vp],-dt,vp)*v2 + dt*beta*g*r*v2) //Buyoancy + on(C1,C2,C3,C4,u2=0) + on(C1,C2,C3,C4,u1=0); //no-slip on box surface // // Define convection-diffusion problem for heat. // problem convectdiffusion(r,rr,solver=UMFPACK, init=reuseMatrix) = int2d(Th)( r*rr //Remains from advective part in lagrangian coordinates + a*dt*(dx(r)*dx(rr)+dy(r)*dy(rr))) //Diffusion part - int2d(Th)(convect([up,vp],-dt,rp)*rr) //Advection Lagrangian coordinates + on(C2,C3,C4,r=0) + on(C1,r=Thot-Tcold); // // in order to calculate the heat transfer coefficient // Vh o=1; real lenfloor = int1d(Th,C1,qfe=qf2pE)(o); // // Define problem for stream function. // problem streamlines(psi,phi) = int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi)) + int2d(Th)( -phi*(dy(u1) - dx(u2))) + on(C1,psi=0.); cout.precision(3); for ( real t = 0.0; t < tend; t = t + dt ) { NS; up=u1; vp=u2; //plot([u1,u2],p,value=1, cmm="t="+t); convectdiffusion; rp=r; plot(r, wait = false, value=false, nbiso=20, fill=true, cmm="Temperature rise r at t="+t); reuseMatrix=true; cout << "Timestamp:\t\t" << t << "sec / " << tend << "sec" << endl; streamlines; //plot(psi,nbiso=20, fill=false,value=true,cmm="Streamlines [(m^3/s)/m]; t="+t+"sec; alpha="+alpha+"W/m^2K",wait=0); n=n+1; } qx = -lambda*dx(r); qy = -lambda*dy(r); real heatflux = int1d(Th,C1,qfe=qf2pE)(qx*N.x+qy*N.y); real alpha = -heatflux/(lenfloor*(Thot-Tcold)); cout << endl; cout << "Heat transfer:\t\t" << heatflux << "W/m" << endl; cout << "Heat transfer coef:\t" << alpha << "W/m^2K" << endl; cout << "Prandtl Number:\t\t" << nu/a << endl; cout << "Beta:\t\t\t" << beta << "1/K" << endl; cout << "Thermal diffusivity:\t" << a << "m^2/s" << endl; cout << endl; cout << "NbBoundaryElements:\t" << Th.nbe << endl; cout << "NbTriangles:\t\t" << Th.nt << endl; cout << "NbVertices:\t\t" << Th.nv << endl; cout << "nDoF (T & Psi):\t\t" << Qh.ndof << endl; cout << "nDoF (u1,u2):\t\t" << Vh.ndof << endl; cout << endl; Tfluid = r + Tcold; // // Write data to files. // savemesh ( Th, "convective_rolls.msh" ); ffSaveVh ( Th, Vh, "convective_rolls_vh.txt" ); ffSaveVh ( Th, Qh, "convective_rolls_qh.txt" ); ffSaveData ( psi, "convective_rolls_psi.txt" ); ffSaveData ( Tfluid, "convective_rolls_temperature.txt" ); ffSaveData2 ( u1, u2, "convective_rolls_ux_uy.txt" ); // // Terminate. // cout << "\n"; cout << "convective_rolls:\n"; cout << " Normal end of execution.\n"; exit ( 0 );