Hi, i am trying to solve NS equation with Newton method in a 3D cubic domain with non zero value of the y-component of the velocity at the inflow boundary.
The problem I have is that the velocity is almost not updating at each iteration and each component of the final velocty has exactly the same values.
which could be the reason of this behaviour?
here is the code I ame running:
real L = 1; // Characteristic length of the cell (m)
real mu = 1e-1; // Air viscosity at 300K
real U = 0.6; // Characteristic velocity and inlet boundary velocity (Dirichlet condition NS)
real absu = 0.; // Module of the velocity
real Re = rhoaUL/mu;
load “msh3”;
load “iovtk”;
load “medit”; // Per visualizzazione 3D con medit
//Mesh size
int meshSize = 20;
// Mesh generation
mesh3 Th = cube(meshSize, meshSize, meshSize, [Lx,Ly,L*z]);
// Visualize the mesh using the plot command
plot(Th,cmm = “3D Mesh Visualization”);
fespace Vh(Th, [P2, P2, P2]);
Vh [ux,uy,uz] , [vx,vy,vz]; //Fluid velocity, test function
Vh [ux1,uy1,uz1] , [dux,duy,duz]; //Fluid Velocity at the next iteration, variation between one iteration and another
fespace Qh(Th, P1);
Qh p, q; // pressure e test function
Qh p1; // pressure at the next iteration
Qh Q;
// Stabilization term definitions
Qh hK = hTriangle; // stabilization parameter based on the element size (hTriangle is usually a built-in function in FreeFem++ that returns the size of the element).
Qh tauK; //will be used to store stabilization coefficients.
//Macro NB: ‘* rappresenta il prodotto scalare
macro grad(v) [dx(v), dy(v), dz(v)] //
macro div(vx,vy,vz) ( dx(vx) + dy(vy) + dz(vz) ) //
macro ugradv(ux,uy,uz,vx,vy,vz) [ [ux,uy,uz]’*grad(vx) , [ux,uy,uz]‘*grad(vy) , [ux,uy,uz]’*grad(vz) ] //
func uIn = U; // Dirichlet BC
// Stokes (linear) problem, used to initialized Newton’s method:
problem Stokes( [ux,uy,uz,p] , [vx,vy,vz,q], solver=sparsesolver )
= int3d(Th)( 1/Re*(grad(ux)‘*grad(vx) + grad(uy)’*grad(vy) + grad(uz)'grad(vz) )
- div(vx,vy,vz) * p
+ div(ux, uy, uz) * q ) / l.h.s. della weak formulation /
+ on(2,4,5,6, ux = 0 , uy = 0 , uz = 0 ) / Dirichlet BCs on the walls /
+ on(1, ux = 0 , uy = uIn , uz = 0 ); / Dirichlet BCs at the inflow /
/ Neumann BCs at the outflow (automatically set by FreeFem) */
// Linearized Navier Stokes:
problem linNS( [ux1,uy1,uz1,p1] , [vx,vy,vz,q], solver=sparsesolver )
= int3d(Th)(ugradv(ux,uy,uz,ux1,uy1,uz1)' * [vx,vy,vz] +
ugradv(ux1,uy1,uz1,uz,uy,uz)' * [vx,vy,vz] +
1/Re * ( grad(ux1)'*grad(vx) + grad(uy1)'*grad(vy) + grad(uz1)'*grad(vz) ) -
div(vx,vy,vz) * p1 +
div(ux1, uy1, uz1) * q)
- int3d(Th)(ugradv(ux,uy,uz,ux,uy,uz)' * [vx,vy,vz])
+ on(2,4,5,6, ux1 = 0 , uy1 = 0 , uz1 = 0 )
+ on(1, ux1 = 0 , uy1 = uIn , uz1 = 0 );
real toll = 1e-06; // Error variation tolerance (of fluid velocity u)
int kmax = 30; // Maximum number of iteration for Newton to diverge
// Macro containing Newton’s iterations:
macro NavierStokes(){
Stokes;
for(int k = 0; k < kmax; k++){
linNS;
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) ) );
ux = ux1;
uy = uy1;
uz = uz1;
p = p1;
if( err < toll ){
cout << "Value of the error after " << k+1 << " iterations: " << err << endl ;
break;
}
if( k == kmax - 1 )
cout << “WARNING: Newton method does not converge in " << kmax << " iterations!” << endl;
}
} //
// Non-linear problem solution with Newton’s method:
NavierStokes;
absu = sqrt(int3d(Th)([ux,uy,uz]'[ux,uy,uz]))+ 1e-8; // norma L2 del campo di velocità su tutto il dominio, che è una misura dell’energia cinetica totale del flusso.
tauK = hK[].max/(2absu);
cout << "|u|= " << absu << endl;
ofstream fileOut(“Pressure_velocity.txt”); // Create a text file for Temperature
// Loop through the mesh points and save the scalar values
for (int i = 0; i < Th.nv; i++) {
// Access the x, y, z coordinates and corresponding scalar value
real x = Th(i).x;
real y = Th(i).y;
real z = Th(i).z;
real valueP = p[i]; // Access the value of the scalar field at each node
real valueux = ux[i]; // Access the value of the scalar field at each node
real valueuy = uy[i]; // Access the value of the scalar field at each node
real valueuz = uz[i]; // Access the value of the scalar field at each node
// Write coordinates and value to the file
fileOut << x << " " << y << " " << z << " " << valueP << " " << valueux << " " <<valueuy<< " " <<valueuz << endl;
};
plot(p, cmm = “pressure Solution”, value = true, fill=true, nbiso = 30);
Thank you in advance!