Dear Michael,

thank you for your reply, but actually I have mainly  two doubts.

1) First doubt regards the assignment of Neumann boundary conditions. As an example, free slip conditions along an horizontal plane (x,y) require that tangential stresses are null sigma_{xz}=sigma_{yz}=0. Typically the other components of the stress are not assigned. Does this mean that they must be taken null, or that they must be computed from the actual solution?

2) If the non specified components of stress are to be computed from the actual solution, we have to calculate the stress tensor that involves derivatives of the velocity. In finite element formulation derivatives are obtained by interpolation of the nodal values of the dof (in this case velocity) with the derivatives of the shape functions. How can be defined the derivative of shape functions over a face of a 3D element? I'm not able to find a proper mapping between global 3D mesh coordinates (x,y,z) and local 2D representation of the element face (xi,eta).

Antonella

Dear S. M. H. Goushegir,

thank you for your answer. Unfortunately I still have doubts. Are you suggesting me to perform computations on the Hexahedron face of the reference element maintaining the 3D coordinates? As far as I know, when you map a the global coordinates of an element face into the reference coordinates of a face, the number of coordinates decreases of 1.

As an example, in 2D coordinates, global elements are general quadrangles in (x,y), and the reference element is a square in (ξ,η). Global element sides are lines in (x,y) coordinates and I map them into a reference 1D line along the ξ coordinate. For the line, shape functions are N(1)=0.5*(1-ξ), N(1)=0.5*(1+ξ). Their derivatives with respect to ξ are -0.5 and +0.5. With respect to global coordinates, if dxdξ > 0, I get dN(1)/dx = -1.0/(x(2) - x(1)), dN(2)/dx = 1.0/(x(2) - x(1)), dN(1)/dy = dN(2)/dy = 0.0. Similar thing if dxdη > 0.

But what about the Hexhaedron?

My approach is to transform the 3D face element to a 2D space. For instance, take a 3D triangle of 3 nodes Xi = (xi, yi, zi), i = 1, 2, 3. Then we could define a new base ex, ey, ez as follows:

ex = (X2 - X1)/norm(X2-X1)

ez = (X2-X1) x (X3-X1) / (norm(X2-X1) * norm(X3-X1))

ey = ez x ex

Then, eliminating the z-component, the nodes in that new space become for instance:

X1 -> (0,0,0)

X2 -> ((X2-X1) . ex, 0, 0)

Z3 -> ((X3-X1) . ex, (X3-X2) . ey, 0)

Now you can use the standard FEM procedure for integration on the 3D triangle.

Hope this helps.

X2