Hello,

I am developing a 3D navier-stokes code. Alhtough is works fine in uniform spacing grids, it seems to give strange results when using uneven spacing grids.

I think it is due to the factor between the volume of the cell and the sum of surface area of the cell (Gauss theorem). As the volumes increases, it's surface does not increase in the same degree. Please check this image to see what I am talking about:

http://xflow.50megs.com/ngrid.GIF

Any help is appreciated.

Thank you in advance,

x-flow.

I am putting a part of the code that calculates the volume by gauss theorem x1, y1, z1 etc are vertices of the triangle of one of the face by the volume in question

gutil.getCrossProduct(A,B,tAr);

gets the cross product of two vectors A(from node n1 to n2) and B(from node n2 to n3) into array tAr ..so tAr[1] , tAr[2] and tAr[3] are x y z components of twice the area of triangle,

then the contribution from this face is

+ (txc * arx1 + tyc * ary1 + tzc * arz1 ) * (1.0 / 3.0) ;

where all the area components are assumed to be pointing out,

add for all the triangles faces those can make the volume, further u can have a look at the code its clear, (cr is volume in right and cl is cell in left of face by n1 n2 n3 with right hand rule)

x1 = *( Xori + n1 -1) ;

y1 = *( Yori + n1 -1) ;

z1 = *( Zori + n1 -1) ;

x2 = *( Xori + n2 -1) ;

y2 = *( Yori + n2 -1) ;

z2 = *( Zori + n2 -1) ;

x3 = *( Xori + n3 -1) ;

y3 = *( Yori + n3 -1) ;

z3 = *( Zori + n3 -1) ;

// first triangle

A.setX(x2 - x1);

A.setY(y2 - y1);

A.setZ(z2 - z1);

B.setX(x3 - x2);

B.setY(y3 - y2);

B.setZ(z3 - z2);

gutil.getCrossProduct(A,B,tAr);

arx1 = tAr[1] * 0.5 ;

ary1 = tAr[2] * 0.5 ;

arz1 = tAr[3] * 0.5 ;

// for volume

txc = (1.0 / 3.0) * (x1 + x2 + x3);

tyc = (1.0 / 3.0) * (y1 + y2 + y3);

tzc = (1.0 / 3.0) * (z1 + z2 + z3);

if (cl > 0) {

*(volume + cl) = *(volume + cl) + (txc * arx1 + tyc * ary1 + tzc * arz1 ) * (1.0 / 3.0) ;

}

if (cr > 0) {

*(volume + cr) = *(volume + cr) - (txc * arx1 + tyc * ary1 + tzc * arz1 ) * (1.0 / 3.0) ;

}

// now store the area

*(arX + i) = arx1;

*(arY + i) = ary1 ;

*(arZ + i) = arz1;

Hi zxaar,

Thanks for your answer but I think the problem is not the value of the volume or area. This code uses simple orthogonal hexahedra (bricks). The method I used to calculate the volume is correct because I have checked it. I decompose every quad face into 2 triangles to calculate the area, and the hexahedra into 12 tetrahedra sharing the centroid of the hexahedra. This makes it easy to calculate both surface area and volume.

I know you can calculate the volume using the Gauss theorem but I did not use this. I believe this is an approximation and introduces some error in non-orthogonal grids.

I was referring to the gauss theorem used in the evaluation of convection and dissipation terms which are the first and second terms of the equation shown in the figure. I think the problem should be in the sum of the area over all the faces of each cell. I decompose the momentum in 3D and use the sum of all faces not only the ones in each direction. For example, for u-velocity (x direction), I use all the faces of the cell (in all directions). Is this correct?

x-flow.

I forgot to say but the vector plotted in the figure is the velocity (u, v, w). You can see that in the cells with larger volume (twice as high) the velocity is 1.34 times larger (in u or x direction) than in smaller cells. Also the factor of (1/V1 * sum Af1) / (1/V2 * sum Af2) is 1.375. If I use evenly spaced grids, this factor is 1 (because all cells have the same volume and surface area) and the velocity seems to be correct.

I meant second and third terms. Sorry.

http://xflow.50megs.com/ngrid.GIF

actually i am not able to see this figure , and perhaps i mistook your message ..can u sen dme the image to my email ..if possible i will be glad to help you out with this

The image is posted at:

http://xflow.50megs.com/

x_flow,

I am not sure your FV equation has all the terms correctly, but might be wrong. If you wish me to clarify, drop me an e-mail. BTW, are you assuming incompressible flow? Laminar? And how is the continuity eq. treated?

Rami