Newton 3 D

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

Page last modified on November 26, 2014, at 03:48 PM