Using Boussinesq hypothesis, the equations in \(\Omega=(0,1)^3\) for velocity and temperature are (\(\vec{e}_3\) OZ-unit vector)
\begin{array}{rcl} (\vec{u}\cdot \vec{\nabla})\vec{u} -\nu \Delta \vec{u}+\vec{\nabla p}& =& c\,T\vec{e}_3,\quad \Omega\\ div(\vec{u})& =& 0,\quad \quad \quad\quad \Omega\\ (\vec{u}\cdot \vec{\nabla})T -\kappa \Delta T& =& 0,\quad \quad \quad \quad \Omega\\ \vec{u} & =& 0,\quad \quad\partial\Omega\\ T & =& 0,\quad\quad \Gamma=\mbox{two opposite walls}\\ \frac{\partial{T}}{\partial\vec{n}} & =& 0,\quad\partial\Omega\setminus \Gamma. \end{array}
include Dirichlet boundary conditions on \(U\) space for the triple \((\vec{u},p,T)\) and maintain the restriction \(div( \vec{u})=0\). Then take
$$F(\vec{u},p,T):=(-\nu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u}+\vec{\nabla p} +c\,T\vec{e}_3,\; div(\vec{u}),(\vec{u}\cdot \vec{\nabla})T -\kappa \Delta T)$$
again, a simple computation gives
$$<dF((\vec{u},p,T)),(\vec{h},q,H)>=(-\nu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u}+\vec{\nabla q},\; div(\vec{h}),\; -\kappa \Delta H+ (\vec{u}\cdot \vec{\nabla})H +(\vec{h}\cdot \vec{\nabla})T) $$
and Newton method is
Want \((\vec{u}_\star,p_\star,T_\star)\in U\) such that \(F(\vec{u}_\star,p_\star,T_\star)=0\). Given \((\vec{u},p,T) \approx (\vec{u}_\star,p_\star,T_\star)\) improve \((\vec{u},p,T)\) by taking \(\vec{u}:=\vec{u}+\vec{h}\), \(p:=p+q\) and \(T:=T+H\) where \((\vec{h},q,H)\) satisfies the linear problem \begin{array}{rcl} (-\nu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u}+\vec{\nabla q}+c\,H\vec{e}_3, \; div(\vec{h}),-\kappa \Delta H+ (\vec{u}\cdot \vec{\nabla})H +(\vec{h}\cdot \vec{\nabla})T )= & &\\ =(-\nu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u}+\vec{\nabla p} +c\,T\vec{e}_3,& div(\vec{u}),&-\kappa \Delta T+ (\vec{u}\cdot \vec{\nabla})T ) \end{array} Solved by a variational formulation using a triple of test function \((\vec{v} ,s,r)\).
Incompressible Navier Stokes and Temperature in a cube using Boussinesq approximation. Continuation on viscosity
Adapted on April 2014 from code by Author: F. Hecht, jan 2012.
![]() | ![]() | ![]() |
Temperature initial | Viscosity 0.01 | Viscosity 0.0025 |
download example: boussinesq3d.edp or return to main page
// // Adapted in April 2014 from a code by F. Hecht // /* Incompressible steady Navier Stokes + boussinesq coupled system with Taylor-Hood Finite element Non-linearity : Newton method continuation on Rayleigh Number Solution computed on a cube [0,1]^3 Problem: Navier-Sokes and Temperature on vertical direccion (Boussinesq hyp.) Equations: (u.grad)u-nu*Lap(u) +Grad p = -c*T*(0,0,1)' in the cube (u.grad)T-k*Lap(T) = 0 in the cube Boundary Conditions: u=0 on boundary of cube T=0 on opposite sidewalls 2 and 4 (d/dn)T=0 on the rest of walls */ // // 2D mesh // int nn=10; mesh th2d=square(nn,nn); // // labels for the cube 3d mesh // int[int] rup=[0,1]; //top of the cube labelled by cero int[int] rdown=[0,1];//bottom of the cube labelled by cero int[int] rmid=[1,1,2,2,3,3,4,4];// sidewalls labelled as the 2d edges they stand on // // 3D mesh // load "msh3" mesh3 Th=buildlayers(square(nn,nn),nn, zbound=[0.,1.], labelmid=rmid, labelup = rup, labeldown = rdown); // // MACROS....they mean // // for a scalar u grad(u) computes the gradient // macro grad(u) [dx(u),dy(u),dz(u)]// // // For a vector u of componentes u#1,u#2,u#3 // here # acts as wildcar character // Grad(u) compute the matrix of the gradients // macro Grad(u) [grad(u#1),grad(u#2),grad(u#3)]// // // For a vector u of componentes u#1,u#2,u#3 // div(u) compute the divergence // macro div(u) (dx(u#1)+dy(u#2)+dz(u#3)) // // // for a vector u of componentes u#1,u#2,u#3 a scalar v // ugrad(u,v) compute a scalar (u.Grad)v // macro ugrad(u,v) ([u#1,u#2,u#3]'*grad(v) ) // // // for a vector u of componentes u#1,u#2,u#3 // for a vector v of componentes v#1,v#2,v#3 // the temperature T // UgradV(u,v,T) gives a vector of componente // (u.grad)v#1, (u.grad)v#2, (u.grad)v#3, (u.grad)T // // This is a non linear term when u=v // macro UgradV(u,v,T) [ugrad(u,v#1),ugrad(u,v#2),ugrad(u,v#3),ugrad(u,T)]// // // Finite element spaces for computations // // Taylor-Hood P2-P1 pair for better precision. Migth be a bit expensive and // could be better to use P1b-P1 // fespace XXXMh(Th,[P2,P2,P2,P1,P1]);//vectorial FE space for velocity and pressure XXXMh [u1,u2,u3,p,T]; // XXXMh [uw1,uw2,uw3,pw,Tw]; //increments during the Newton iteration XXXMh [v1,v2,v3,q,TT]; // test functions XXXMh [uc1,uc2,uc3,pc,Tc]; // Current computed values out of the Newton iterate // // Physical constants // // Physical parameter // real nu= 1./100; // starting viscosity real nufinal=1/1000.; // wanted final viscosity // // cnu reducing factor for viscosity // We reduce viscosity by dividing by cnu // real cnu=2; // // Number of reductions on viscosity // int countermax=10; int counter=0; // real Pr = 56.2; // Prandtl number real k = 1. / Pr ; // Heat diffusion, inverse of Prandtl // // verbosity=0; real err=0; // stop test for Newton real eps=1e-6; real nuc=nu;// current viscosity cout << " ---------------------------- " << endl; cout << " --- START COMPUTATIONS --- " << endl; cout << " ---------------------------- " << endl; // bool adaptation=1; int one =adaptation; // to perform mesh adaptation // // initial data for Newton process // [u1,u2,u3,p,T]=[0,0,0,0,x]; // cero velocity and heating from above plot(T,wait=0,cmm="Initial temperature T=x", coef=0.3,fill=1,value=1,nbiso=20); // // Continuation on viscosity number // Compute first for a given viscosity and then after convergence of // Newton method reduce the viscosity number until nufinal // while((nuc>nufinal) && (counter<countermax)&&(cnu>1.05)) { // start loop on viscosity counter =counter+1; real rey=1./nuc; real c= - rey/Pr; // ratio Reynolds/Prandtl int n; // counter for Newton iteration for (n=0;n<=20;n++)// nonlinear Newton iteration { // initial data very first iteration is //[u1,u2,u3,p,T]=[0,0,0,0,1-x]; // // variables uw1,uw2,uw3,pw,Tw contain the increments computed by // using the Gateaux o Frechet derivative // // v represents---> [v1,v2,v3] (test functions) // TT test function for temperature // solve BoussinesqNL([uw1,uw2,uw3,pw,Tw],[v1,v2,v3,q,TT]) = // start now with the term // <dF(u,T),(uw,Tw)> // int3d(Th) ( // //result of the nonlinear parts on the differential mappping //nonlinear parts (u*grad)uw+(uw*grad)u+(u*grad)Tw+(uw*grad)T // // // u represents---> [u1,u2,u3] // uw represents---> [uw1,uw2,uw3] // // Solving <dF(u,T),(uw,Tw)>=-F(u,T) using variational formulation // and then // u:=u+uw // T:=T+Tw // UgradV(u,uw,Tw)' * [v1,v2,v3,TT] + UgradV(uw,u,T)' * [v1,v2,v3,TT] // // nu*Lap(uw)*v (here uw represents---> [uw1,uw2,uw3]) // + ( Grad(uw):Grad(v)) * nuc // // effect of Boussinesq hypothesis body force is -c*T*(0,0,1)' // + c*Tw*v3 // //-k*Lap(Tw) + grad(Tw)'*grad(TT)*k // // incompressibility for uw // - div(uw)*q // // grad(pw) // -div(v)*pw // // stabilization term // + 1e-8*pw*q ) // // Now the -F(u,T) term goes +F(u,T) on left hand side + int3d(Th)( // // non linear term on data (u,T) // (u*grad)u+(u*grad)T // UgradV(u,u,T)' * [v1,v2,v3,TT] // non linear term on data (u,T) // //the -nu*Lap(u) // + ( Grad(u):Grad(v) ) * nuc // // the body force c*T*(0,0,1)' // + c*T*v3 // // the -k*Lap(T) // + grad(T)'*grad(TT)*k // // the incompressibility condition // - div(u)*q // // the grad(p) // -div(v)*p // the stabilization term + 1e-8*p*q ) // Boundary Conditions for velocity field + on(1,2,3,4, uw1=0,uw2=0,uw3=0)// velocity cero top, bottom and sidewalls // Boundary Conditions for Temperature field + on(2,4,Tw=0);// Temperature cero opposite sidewalls 2 and 4 // on the rest...normal derivative of T is cero. // // update all the fields....velocities and temperature whithin Newton iteration // // // update values to compute again // u1[] += uw1[]; //u1[]=u1[]+uw1[] an update the rest in the same way // u2[] += uw2[]; //u2[]=u2[]+uw2[] // u3[] += uw3[]; //u3[]=u3[]+uw3[] // p[] += pw[]; //p[]=p[]+pw[] // T[] += Tw[]; //T[]=T[]+Tw[] // // computation of norms for the increments, takes five components at the same time // real solutionl2=u1[].l2; // l2 norm of all components [u1,u2,u3,p,T] real incrementl2=uw1[].l2; // l2 norm of all components [uw1,uw2,uw3,pw,Tw] // // normalization of the norms // err= incrementl2/solutionl2; cout << " iter = "<< n << " norm l2 of the solution is = " << solutionl2 << " nu = " << nuc << endl; cout << " iter = "<< n << " norm l2 of the increment is = " << incrementl2 << " nu = " << nuc << endl; cout << " iter = "<< n << " err " <<err << " nu = " << nuc << endl; // // plot iterate // plot(T,wait=0,cmm="Temperature using Newton Method for viscosity = " + nuc +", it = "+ n+ ", err = "+err, coef=0.3,fill=1,value=1,nbiso=20); // //--> In case of no convergence... // (n>3 AND err>2) is blowup // (n>12 AND err>eps) is lack of convergence or very slow // if( (n>3 && err > 2.) || (n>12 && err > 1) ) // (n>3 AND err>2) or (n>12 AND err>eps) {// start ---> if cout<<" NO convergence for nu = "+nuc+ ", it = "+ n+ ", err = "+err <<endl; nuc=nuc*cnu; // get back to previous nu cout << " LAST viscosity TO CONVERGE = " << nuc<< endl; // // take computed values for previous nuc to start again // u1[]=uc1[]; // // New reduction for cnu // cnu= cnu^(0.6); // no conv. => change lower (cnu^(0.6) makes cnu smaller) nuc = nuc/cnu; // new vicosity cout << " RESTART with NEW viscosity = " << nuc << endl; cout << " Using as Inital Data the solution for PREVIOUS viscosity " << endl; break; // Blowup ???? --->out of for loop go to while loop }// end ---> if( n>3 && err > 2.) // //-->In case of Convergence then... // if(err < eps) {// start ---> if(err < eps) cout << " CONVERGENCE for viscosity = " << nuc << " on iters = "<< n <<", err = "+err << endl; //uc1=u1; uc2=u2; uc3=u3; pc=p; Tc=T; // save computed values of this convergent solution uc1[]=u1[]; plot(T,wait=0,cmm="Solution for viscosity = " + nuc,coef=0.3,fill=1,value=1,nbiso=20); if( n < 4) cnu=cnu^1.5; // fast converge => change faster ( cnu^ 1.5 makes bigger cnu) nuc = max(nufinal, nuc/cnu); // reduce to a new vicosity cout << " Continuation to new viscosity = " << nuc <<" (dividing factor "<<cnu<<")"<< endl; cout << " Initial data IS the FINAL computation for previous viscosity " << endl; break; // converge ---> out of for loop. Start again with u1,u2,p and new nuc }// end ---> if(err < eps) }//end ---> for loop } // end while loop on viscosity plot(T,cmm="Temperature Last converge of Newton Method for nu = " + nuc,coef=0.3,wait=0,fill=1,value=1,nbiso=10); cout << " LAST viscosity TO CONVERGE nu = " << nuc << endl; plot(coef=0.2,cmm=" Temperature for viscosity "+nuc,T,fill=1,value=1,nbiso=20);