Hello,
I am trying to write Freefem++ code for the steady Navier Stokes equations in 3D using Newton’s method taking my cue from the 2D model in the FreeFem documentation. The code I have written is as follows:
load “msh3”
load “medit”
//Data
real mu=1.825*10^(-5); // viscosity
real g=9.81; //gravitational acceleration
real uMax = 10. ;// Inlet boundary velocity (for Dirichlet condition)
//Mesh Generation
real L = 1;
mesh3 Th = cube(20, 20, 20, [Lx,Ly,L*z]);
medit(“Mesh”, Th);
// {2,4,5,6} = wall; 1 = inflow boundary ; 3 = outflow boundary
//Finite Elements Spaces:
fespace Vh(Th, [P2, P2, P2]);
Vh [ux,uy,uz] , [vx,vy,vz] ; // Fluid velocity, test function
Vh [ux1,uy1,uz1] , [dux,duy,duz] ;
fespace Qh(Th, P1);
Qh p , q ; // pressure e test function
Qh p1 ;
Qh rhoa = 1.293; /* Just to define the macro, in this case I am assuming rhoa constant but in further implementations I shall need to define more complex functions for rhoa, so I want to keep the code as general as possible */
//Macro
macro grad(v) [dx(v), dy(v), dz(v)] //
macro div(vx,vy,vz) ( dx(vx) + dy(vy) + dz(vz) ) //
macro DivRho(vx,vy,vz,rho) ( grad(rho)’ * [vx,vy,vz] + rho * div(vx,vy,vz) ) //
macro ugradv(ux,uy,uz,vx,vy,vz) [ [ux,uy,uz]‘*grad(vx) , [ux,uy,uz]’*grad(vy) , [ux,uy,uz]'*grad(vz) ] //
// Macro for solving the steady nonlinear problem using Newton method:
real toll = 1e-09; // Tolerance (for fluid velocity u)
int kmax = 15; // max number of iterations
macro NavierStokes(){
Stokes; // I start solving Stokes
for(int k = 0 ; k < kmax ; k++ ){
linNS; // I compute u_(k+1) using linNS, with u_k taken from the previous iteration
// I compute the error:
dux[] = ux1[] - ux[];
duy[] = uy1[] - uy[];
duz[] = uz1[] - uz[];
real err = sqrt( int3d(Th) ( grad(dux)'*grad(dux) + grad(duy)'*grad(duy) + grad(duz)'*grad(duz) ) ) /
sqrt( int3d(Th) ( grad(ux)'*grad(ux) + grad(uy)'*grad(uy) + grad(uz)'*grad(uz) ) );
// I update the fluid velocity and the pressure:
ux[] = ux1[];
uy[] = uy1[];
uz[] = uz1[];
p[] = p1[];
if( err < toll ){
// I print the error at the last iteration:
cout << "Error after " << k+1 << " iterations: " << err << endl ;
break;
}
if( k == kmax - 1 )
cout << "WARNING: Newton method does not converge in " << kmax << " iterations!" << endl;
}
} //
func h = 0; // Neumann boundary datum
func uIn = uMax; // Dirichlet bounday datum (both are constant functions)
// Stokes problem (cut out the nonlinear term):
problem Stokes( [ux,uy,uz,p] , [vx,vy,vz,q] )
= int3d(Th)( mu * ( grad(ux)'*grad(vx) + grad(uy)'*grad(vy) + grad(uz)'*grad(vz) ) +
mu * div(ux,uy,uz) * div(vx,vy,vz) -
div(vx,vy,vz) * p +
DivRho(ux, uy, uz, rhoa) * q ) // l.h.s. della weak formulation
- int3d(Th)( rhoa * (-g) * vz ) // r.h.s. della weak formulation
+ on(2,4,5,6, ux = 0 , uy = 0 , uz = 0 ) // Dirichlet BCs at the wall
+ on(1, ux = 0 , uy = uIn , uz = 0 ) //Dirichlet BCs at the inflow
+ int2d(Th,3)( h*vy ); // Neumann BCs at the outflow
// Navier Stokes linarized with semi-implicit expression of the nonlinear term:
problem linNS( [ux1,uy1,uz1,p1] , [vx,vy,vz,q] )
= int3d(Th)(
rhoa * ugradv(ux,uy,uz,ux1,uy1,uz1)' * [vx,vy,vz] +
mu * ( grad(ux1)'*grad(vx) + grad(uy1)'*grad(vy) + grad(uz1)'*grad(vz) ) +
mu * div(ux1,uy1,uz1) * div(vx,vy,vz) -
div(vx,vy,vz) * p1 +
DivRho(ux1, uy1, uz1, rhoa) * q )
- int3d(Th)( rhoa * (-g) * vz )
+ on(2,4,5,6, ux1 = 0 , uy1 = 0 , uz1 = 0 )
+ on(1, ux1 = 0 , uy1 = uIn , uz1 = 0 )
+ int2d(Th,3)( h*vy );
// Finally I solve the problem:
NavierStokes;
// Plot:
plot( p , ps=“Power.ps”);
plot( [ux,uy,uz] , ps=“FluidVelocity.ps”);
Nevertheless, the program does not compile and I do not understand what mistake I am making since I followed step by step the model found in the documentation (obviously adapting it to my case).
Thank you in advance for your help.
Best regards,
Marianna