Hi, I am doing 3D simulation using SOLA-VOF and NASA-VOF. Does anyone know how to extend the 15 cases for 2D free surface boudary conditions in SOLA-VOF into 63 cases for 3D BC? I have do so but it seems that there are some problems in my extension of the BC. Thank you!

i paste my code (in C language) for implementation of free surface bc in 3D using 64 case method but i think that there are better methods (recently i don't use this).

/*

Free surface boundary condition using 64 case configuration

*/

#include "free_surface_bc.h

#define U(a,b,c) u[a+1][b][c] #define V(a,b,c) v[a][b+1][c] #define W(a,b,c) w[a][b][c+1]

/* macro definition of possible 64 case orientation of free surface cell in 3d */

#define NF_000000 0x000000 #define NF_000001 0x000001 #define NF_000010 0x000010 #define NF_000011 0x000011 #define NF_000100 0x000100 #define NF_000101 0x000101 #define NF_000110 0x000110 #define NF_000111 0x000111 #define NF_001000 0x001000 #define NF_001001 0x001001 #define NF_001010 0x001010 #define NF_001011 0x001011 #define NF_001100 0x001100 #define NF_001101 0x001101 #define NF_001110 0x001110 #define NF_001111 0x001111 #define NF_010000 0x010000 #define NF_010001 0x010001 #define NF_010010 0x010010 #define NF_010011 0x010011 #define NF_010100 0x010100 #define NF_010101 0x010101 #define NF_010110 0x010110 #define NF_010111 0x010111 #define NF_011000 0x011000 #define NF_011001 0x011001 #define NF_011010 0x011010 #define NF_011011 0x011011 #define NF_011100 0x011100 #define NF_011101 0x011101 #define NF_011110 0x011110 #define NF_011111 0x011111 #define NF_100000 0x100000 #define NF_100001 0x100001 #define NF_100010 0x100010 #define NF_100011 0x100011 #define NF_100100 0x100100 #define NF_100101 0x100101 #define NF_100110 0x100110 #define NF_100111 0x100111 #define NF_101000 0x101000 #define NF_101001 0x101001 #define NF_101010 0x101010 #define NF_101011 0x101011 #define NF_101100 0x101100 #define NF_101101 0x101101 #define NF_101110 0x101110 #define NF_101111 0x101111 #define NF_110000 0x110000 #define NF_110001 0x110001 #define NF_110010 0x110010 #define NF_110011 0x110011 #define NF_110100 0x110100 #define NF_110101 0x110101 #define NF_110110 0x110110 #define NF_110111 0x110111 #define NF_111000 0x111000 #define NF_111001 0x111001 #define NF_111010 0x111010 #define NF_111011 0x111011 #define NF_111100 0x111100 #define NF_111101 0x111101 #define NF_111110 0x111110 #define NF_111111 0x111111

LOOP_OVER_SURFACE_CELLS {

switch(configuration) {

// one empty neighbor

case NF_100000: U(i-1,j,k) = U(i,j,k) + dx*((V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_010000: U(i,j,k) = U(i-1,j,k) -dx*((V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_001000: V(i,j-1,k) = V(i,j,k) + dy*((U(i,j,k)-U(i-1,j,k))/dx + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_000100: V(i,j,k) = V(i,j-1,k) -dy*((U(i,j,k)-U(i-1,j,k))/dx + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_000010: W(i,j,k-1) = W(i,j,k) + dz*((U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_000001: W(i,j,k) = W(i,j,k-1) -dz*((U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy); break;

// two empty neighbors (opposite)

case NF_110000: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; div = (U(i,j,k)-U(i-1,j,k))/dx +(V(i,j,k)-V(i,j-1,k))/dy +(W(i,j,k)-W(i,j,k-1))/dz; U(i-1,j,k) += 0.5*dx*div; U(i,j,k) -= 0.5*dx*div; break;

case NF_001100: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; div = (U(i,j,k)-U(i-1,j,k))/dx+(V(i,j,k)-V(i,j-1,k))/dy +(W(i,j,k)-W(i,j,k-1))/dz; V(i,j-1,k) += 0.5*dy*div; V(i,j,k) -= 0.5*dy*div; break;

case NF_000011: W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; div = (U(i,j,k)-U(i-1,j,k))/dx +(V(i,j,k)-V(i,j-1,k))/dy +(W(i,j,k)-W(i,j,k-1))/dz; W(i,j,k-1) += 0.5*dz*div; W(i,j,k) -= 0.5*dz*div; break;

// two empty neighbors (adjacent)

case NF_101000: U(i-1,j,k)=U(i,j,k)+0.5*dx*((W(i,j,k)-W(i,j,k-1))/dz); V(i,j-1,k)=V(i,j,k)+0.5*dy*((W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_100100: U(i-1,j,k)=U(i,j,k)+0.5*dx*((W(i,j,k)-W(i,j,k-1))/dz); V(i,j,k)=V(i,j-1,k)-0.5*dy*((W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_011000: U(i,j,k)=U(i-1,j,k)-0.5*dx*((W(i,j,k)-W(i,j,k-1))/dz); V(i,j-1,k)=V(i,j,k)+0.5*dy*((W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_010100: U(i,j,k)=U(i-1,j,k)-0.5*dx*((W(i,j,k)-W(i,j,k-1))/dz); V(i,j,k)=V(i,j-1,k)-0.5*dy*((W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_001010: V(i,j-1,k)=V(i,j,k)+0.5*dy*((U(i,j,k)-U(i-1,j,k))/dx); W(i,j,k-1)=W(i,j,k)+0.5*dz*((U(i,j,k)-U(i-1,j,k))/dx); break;

case NF_001001: V(i,j-1,k)=V(i,j,k)+0.5*dy*((U(i,j,k)-U(i-1,j,k))/dx); W(i,j,k)=W(i,j,k-1)-0.5*dz*((U(i,j,k)-U(i-1,j,k))/dx); break;

case NF_000110: V(i,j,k)=V(i,j-1,k)-0.5*dy*((U(i,j,k)-U(i-1,j,k))/dx); W(i,j,k-1)=W(i,j,k)+0.5*dz*((U(i,j,k)-U(i-1,j,k))/dx); break;

case NF_000101: V(i,j,k)=V(i,j-1,k)-0.5*dy*((U(i,j,k)-U(i-1,j,k))/dx); W(i,j,k)=W(i,j,k-1)-0.5*dz*((U(i,j,k)-U(i-1,j,k))/dx); break;

case NF_100010: U(i-1,j,k)=U(i,j,k)+0.5*dx*((V(i,j,k)-V(i,j-1,k))/dy); W(i,j,k-1)=W(i,j,k)+0.5*dz*((V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_100001: U(i-1,j,k)=U(i,j,k)+0.5*dx*((V(i,j,k)-V(i,j-1,k))/dy); W(i,j,k)=W(i,j,k-1)-0.5*dz*((V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_010010: U(i,j,k)=U(i-1,j,k)-0.5*dx*((V(i,j,k)-V(i,j-1,k))/dy); W(i,j,k-1)=W(i,j,k)+0.5*dz*((V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_010001: U(i,j,k)=U(i-1,j,k)-0.5*dx*((V(i,j,k)-V(i,j-1,k))/dy); W(i,j,k)=W(i,j,k-1)-0.5*dz*((V(i,j,k)-V(i,j-1,k))/dy); break;

// three empty neighbors (adjacent/liquid corner)

case NF_010101: U(i,j,k) = U(i-1,j,k); V(i,j,k) = V(i,j-1,k); W(i,j,k) = W(i,j,k-1); break;

case NF_100101: U(i-1,j,k) = U(i,j,k); V(i,j,k) = V(i,j-1,k); W(i,j,k) = W(i,j,k-1); break;

case NF_011001: U(i,j,k) = U(i-1,j,k); V(i,j-1,k) = V(i,j,k); W(i,j,k) = W(i,j,k-1); break;

case NF_101001: U(i-1,j,k) = U(i,j,k); V(i,j-1,k) = V(i,j,k); W(i,j,k) = W(i,j,k-1); break;

case NF_010110: U(i,j,k) = U(i-1,j,k); V(i,j,k) = V(i,j-1,k); W(i,j,k-1) = W(i,j,k); break;

case NF_100110: U(i-1,j,k) = U(i,j,k); V(i,j,k) = V(i,j-1,k); W(i,j,k-1) = W(i,j,k); break;

case NF_011010: U(i,j,k) = U(i-1,j,k); V(i,j-1,k) = V(i,j,k); W(i,j,k-1) = W(i,j,k); break;

case NF_101010: U(i-1,j,k) = U(i,j,k); V(i,j-1,k) = V(i,j,k); W(i,j,k-1) = W(i,j,k); break;

// three empty neighbors (two opposite)

case NF_111000: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j-1,k) = V(i,j,k) + dy*((U(i,j,k)-U(i-1,j,k))/dx + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_110100: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j,k) = V(i,j-1,k) -dy*((U(i,j,k)-U(i-1,j,k))/dx + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_110010: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; W(i,j,k-1) = W(i,j,k) + dz*((U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_110001: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; W(i,j,k) = W(i,j,k-1) -dz*((U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_101100: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; U(i-1,j,k) = U(i,j,k) + dx*((V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_011100: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; U(i,j,k) = U(i-1,j,k) -dx*((V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_001110: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; W(i,j,k-1) = W(i,j,k) + dz*((U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_001101: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; W(i,j,k) = W(i,j,k-1) -dz*((U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_100011: W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; U(i-1,j,k) = U(i,j,k) + dx*((V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_010011: W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; U(i,j,k) = U(i-1,j,k) -dx*((V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_001011: W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; V(i,j-1,k) = V(i,j,k) + dy*((U(i,j,k)-U(i-1,j,k))/dx + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_000111: W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; V(i,j,k) = V(i,j-1,k) -dy*((U(i,j,k)-U(i-1,j,k))/dx + (W(i,j,k)-W(i,j,k-1))/dz); break;

// four empty neighbors (two opposite pairs)

case NF_111100: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; div = (U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy +(W(i,j,k)-W(i,j,k-1))/dz; U(i-1,j,k) += 0.25*dx*div; U(i,j,k) -= 0.25*dx*div; V(i,j-1,k) += 0.25*dy*div; V(i,j,k) -= 0.25*dy*div; break;

case NF_110011: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; div = (U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz; U(i-1,j,k) += 0.25*dx*div; U(i,j,k) -= 0.25*dx*div; W(i,j,k-1) += 0.25*dz*div; W(i,j,k) -= 0.25*dz*div; break;

case NF_001111: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; div = (U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz; V(i,j-1,k) += 0.25*dy*div; V(i,j,k) -= 0.25*dy*div; W(i,j,k-1) += 0.25*dz*div; W(i,j,k) -= 0.25*dz*div; break;

// four empty neighbors (one opposite pair)

case NF_110101: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j,k)=V(i,j-1,k)-0.5*dy*((U(i,j,k)-U(i-1,j,k))/dx); W(i,j,k)=W(i,j,k-1)-0.5*dz*((U(i,j,k)-U(i-1,j,k))/dx); break;

case NF_110110: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j,k)=V(i,j-1,k)-0.5*dy*((U(i,j,k)-U(i-1,j,k))/dx); W(i,j,k-1)=W(i,j,k)+0.5*dz*((U(i,j,k)-U(i-1,j,k))/dx); break;

case NF_111001: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j-1,k)=V(i,j,k)+0.5*dy*((U(i,j,k)-U(i-1,j,k))/dx); W(i,j,k)=W(i,j,k-1)-0.5*dz*((U(i,j,k)-U(i-1,j,k))/dx); break;

case NF_111010: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j-1,k)=V(i,j,k)+0.5*dy*((U(i,j,k)-U(i-1,j,k))/dx); W(i,j,k-1)=W(i,j,k)+0.5*dz*((U(i,j,k)-U(i-1,j,k))/dx); break;

case NF_011101: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; U(i,j,k)=U(i-1,j,k)-0.5*dx*((V(i,j,k)-V(i,j-1,k))/dy); W(i,j,k)=W(i,j,k-1)-0.5*dz*((V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_011110: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; U(i,j,k)=U(i-1,j,k)-0.5*dx*((V(i,j,k)-V(i,j-1,k))/dy); W(i,j,k-1)=W(i,j,k)+0.5*dz*((V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_101101: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; U(i-1,j,k)=U(i,j,k)+0.5*dx*((V(i,j,k)-V(i,j-1,k))/dy); W(i,j,k)=W(i,j,k-1)-0.5*dz*((V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_101110: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; U(i-1,j,k)=U(i,j,k)+0.5*dx*((V(i,j,k)-V(i,j-1,k))/dy); W(i,j,k-1)=W(i,j,k)+0.5*dz*((V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_010111: W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; U(i,j,k)=U(i-1,j,k)-0.5*dx*((W(i,j,k)-W(i,j,k-1))/dz); V(i,j,k)=V(i,j-1,k)-0.5*dy*((W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_011011: W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; U(i,j,k)=U(i-1,j,k)-0.5*dx*((W(i,j,k)-W(i,j,k-1))/dz); V(i,j-1,k)=V(i,j,k)+0.5*dy*((W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_100111: W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; U(i-1,j,k)=U(i,j,k)+0.5*dx*((W(i,j,k)-W(i,j,k-1))/dz); V(i,j,k)=V(i,j-1,k)-0.5*dy*((W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_101011: W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; U(i-1,j,k)=U(i,j,k)+0.5*dx*((W(i,j,k)-W(i,j,k-1))/dz); V(i,j-1,k)=V(i,j,k)+0.5*dy*((W(i,j,k)-W(i,j,k-1))/dz); break;

// five empty neighbors

case NF_011111: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; U(i,j,k) = U(i-1,j,k) -dx*((V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_101111: V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; U(i-1,j,k) = U(i,j,k) + dx*((V(i,j,k)-V(i,j-1,k))/dy + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_110111: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; V(i,j,k) = V(i,j-1,k) -dy*((U(i,j,k)-U(i-1,j,k))/dx + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_111011: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; V(i,j-1,k) = V(i,j,k) + dy*((U(i,j,k)-U(i-1,j,k))/dx + (W(i,j,k)-W(i,j,k-1))/dz); break;

case NF_111101: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; W(i,j,k) = W(i,j,k-1) -dz*((U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy); break;

case NF_111110: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; W(i,j,k-1) = W(i,j,k) + dz*((U(i,j,k)-U(i-1,j,k))/dx + (V(i,j,k)-V(i,j-1,k))/dy) ; break;

// six empty neighbors

case NF_111111: U(i-1,j,k) += gravity[0]*dt; U(i,j,k) += gravity[0]*dt; V(i,j-1,k) += gravity[1]*dt; V(i,j,k) += gravity[1]*dt; W(i,j,k-1) += gravity[2]*dt; W(i,j,k) += gravity[2]*dt; break;

}//end switch

}

There was a NASA-VOF/3D code developed.

The funding was stopped when NASA priorities were changed, and the LANL report was very brief. But I think the code was (and maybe IS) available under whatever restrictions were imposed after 09/11/2001.

Try Googling NASA-VOF/3D.

rt: Thank you very much for sharing your codes with me. What is the meaning of "gravity" in the codes? How to define this variable using other variables in SOLA-VOF? I remember there is no such a varible in SOLA-VOF.

Jim: I have NASA-VOF3D codes. But you know the treatment of free surface BC is different in NASA-VOF from that in SOLA-VOF. And NASA-VOF introduces variables: AC, AT,ABK,AR, etc.

Dear

first gravity means: gravitational accelleration vector gravity[0]=g_x, gravity[1]=g_y, gravity[2]=g_z,

it is different from sola code but i belive that it work better, sola is too old.

AC, AT,ABK,AR: these parameters is related to treating complex geometries as partial cell (in this manner you can have sloped boundary that is embedded in the structured grid), for more information read nasa report or refer to FLOW3D site (WWW) and read material related to FAVOR (partial cell) algorithm, in fact FLOW3D is 3D version of solavof code. after mastering this you can modify my free surface bc to incorporate this parameters.

note that nasavof3d is related to cylinderical coordinate and if probably it is better for you to extend solavof2d or nasavof2d to 3d (cartesian coordinate).

---------------

also i belive that 64 case mathod is not appropriate for example it miss symmetry in some test cases, for better method refer to litrature:

1. if you need high accurate method refer to paper due to S. Popinet, S. Zaleski JFM 2002 (2D) for 3D version see Hao and Prosperetti in JCP 2004. They contain velosity extrapolation with enforcing zero tangential stress free surface boundary condition.

2. if you need extrapolating velosity to narrow band near free surface (more than one row, suitable for semi-lagrange interface tracker algorithm) refer to paper due to Adalsteinsson and Sethian JCP 1999 and Sethian SIAM Review 1999.

3. less accurate and simplest method is waighted averaging (it gives reasonble accuracy and preserve symmetry), several litrature use this such as: Kim and Lee part I Int. J. Numer. Meth. Fluids 2003, or Babaei et al Comp. Meth. Appl. Mech. Eng. 2005 (it use solavof and it is better for you).

If you need any of these cited paper notify and give your email.

In the sola-vof documentation (LA8355, August, 1980, p. 36) the variables gx and gy are defined as "body accelerations", a generalization on "gravity".

In NASA-VOF, I think you can set the AC, AT, ABK, and AR to 1.0 (There's also some sort of 'cell volume' associated with the partial cell definitions). I don't know anything about differences in the VOF treatment in the two codes.

rt: Thank you for providing so much information! I have no CFD background. Ask a question: gravity[0]=g_x=0, gravity[1]=g_y = g(9.81), gravity[2]=g_z=0, right?

yes

rt: Thank you again!

Dear Fan, Jim_Park and rt

I have NASA technical report related to solavof2d, nasavof2d and nasavof3d (description of algorithm, variables and subrourine + contain code listing), but i need fortran source code of solavof2d. If you have it please send it to me, insteadly i can send to you PDF of reports (peterpendum@yahoo.com)

good luck.

Peter: I have hard copies of all these source codes. But don't have e-copies. If you cannot get them from others, I can mail you hard copies. I think you can get them from this forum.