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