Dear Prof. Abdul Aziz Jaafar,

thank you for asking.

Yes, I am referring to the Navier-Stokes equations (the example I mentioned was in 2D) of a Newtonian fluid, with equation of state say, of a perfect gas, or with a Boussinesq dependence of density upon temperature.

tau_xx is force per unit area in x direction exerted on the surface with x as normal. It gives not a net dilation as pressure does, but an elongation along the x axis. 

The fluid is isotropic, since it is with isotropic viscosity.

The code is based on the finite element space-time formulation by Prof. Hughes and Prof. Tezduyar and coworkers.

Antonella

Dear Prof. Aziz Jaafar,

I totally agree with your reasoning. You made me think that the dependence of stress components on velocity and pressure has to be accounted for by computing stresses from the actual value of the solution. Setting non-assigned value of stress to zero is incorrect.

Thank you.

Antonella

Dear Prof. Andriy Borisyuk,

thank you for your questions. The geometry may be completely arbitrary, from a perfect square gridded with square elements, to an irregular domain, discretized with general quadrangles. Take for example, the "flat plate", a standard benchmark in compressible fluid-dynamics (please see attached files p.6888 fig.6 (2001CMAME), and p.17 (1998CMAME)). In the flat plate case the domain is a rectangle. Some details follow on the boundary conditions which correspond to infinite, semiinfinite, and finite planes. Top boundary is far enough from the bottom wall to be considered at infinite distance from the perturbation induced by the bottom wall. The boundary conditions are of Dirichlet type, assuming, as I just said, that the flow is unperturbed. Left boundary is a free outflow, so it is to be considered as a far field condition, i.e. infinite distance. Symmetry conditions are set at the bottom boundary before the horizontal wall. Symmetry conditions are those of a semiinfinite plane (as far as I know): no penetrating conditions, i.e.  normal velocity equal to zero (Dirichlet condition) and tangential stress equal to zero (Neumann condition). Bottom wall is a no slip wall, i.e. finite plane. 

All boundaries are fixed (rigid). No elasticity is considered here.

The bottom wall is considered as thin.

Looking forward from your comments,

Antonella 

Dear Antonella,

1) no, you cannot set arbitrarily the other component of the stress tensor but you have to consider the suitable model. We worked on free-surface in this study that can help you:  https://www.researchgate.net/publication/223906138_A_finite_volume_dynamic_large-eddy_simulation_method_for_buoyancy_driven_turbulent_geophysical_flows?ev=prf_pub

2) In any case, the setting of the correct BC.s requires to set properly the derivatives of the shape function. For example, if the velocity (2D case) is prescribed by means of linear shape function on Lagrangian simplex, v(x,y) = N1(x,y)*v1+N2(x,y)*v2+N3(x,y)*v3, and you have to set n.Grad v = f, you simply apply Grad to all the shape functions.

Article A finite volume dynamic large-eddy simulation method for buo...

Dear Antonella, I am not quite sure what kind of boundary you are trying to model. I assume that it is some kind of mirroring condition along the x-axis. The first question is how many boundary conditions you need to get a well-posed problem. (I am referring to linear well-posedness and linear stability of the scheme as the non-linear counterparts are unknown. This does not mean that you have to solve a linear problem.) Since there is no flow through the boundary the correct number of bcs are 4 in a three dimensional setting or 3 in 2D.  (It is the compressible NS?)

By setting the normal velocity to 0 you bound the convective terms. In order to bound the viscous terms you need 3 conditions in 2D which implies that the boundary conditions must be mixed. (Sum of convective and diffusive). The problem with this is that you can not guarantee that the normal velocity remains 0. (The only way you can achieve that is by setting no-slip which is not what you wish to do.) Unless you have known data on some of the gradients a solution would be to set

A(u0-ud)-F=0.

u0 is the numerical solution at the boundary that you wish should have a 0 normal velocity. ud is the same as u0 but you flip the sign of the normal velocity. F is the viscous flux vector.  A is the Roe diffusion matrix.  (Flux Jacobian with eigenvalues set to their aboslute values.)

These conditions will not over specify your problem.

Hope it helps,

Magnus

Dear Filippo,

thank you for suggesting the paper on LES simulation of OML turbulent flows. It was an interesting reading. I read that you are setting free surface at the top of the see, periodic lateral boundary conditions, and solid walls at the bottom. First of all it comes up in my mind a question: how did you technically set the periodic boundary conditions? I faced with that in benchmarks of Kelvin-Helmholtz ans Rayleigh–Taylor instabilities but found no way to set them in an appropriate way in finite elements. Do you have any hint on how to manage the periodic repetition of degree of freedom and the scatter  procedure?

As regards the assignment of Neumann or mixed boundary conditions I agree with you: everything depends upon the particular case under consideration. However, if a Dirichlet condition is assigned on velocity along a direction (e.g. vz), stresses can be assigned only in the orthogonal directions (e.g. surface force components in the (x,y) plane but not in the z direction). In case of free supersonic outflow no boundary conditions has to be prescribed. So, in the end, everything depends on the model. Is it this that you mentioned?

As regards the computation of the velocity derivatives on the boundary, I totally agree with you that the Grad operator has to be applied to the shape functions following your formula of interpolation. I am using Lagrangian polynomials. But is impossible to apply the Grad operator directly in the global coordinates. It would be too a complex calculation. The mapping with the reference element allows to transfer computation in a viable system of coordinates. My problem is exactly in the mapping between the global (x,y,z) and the reference (xi,eta,zeta) coordinates when the domain is restricted to the face of a 3D Hexhaedron. The face in global coordinates is fully 3D, while the corresponding face of the reference element has one coordinate that is fixed (i.e. eta = 1) so it is in 2D. The Jacobian of the transformation is no more a 3x3 square matrix, but it is a 2x3 non-invertible matrix. This is the problem I'm tackling. Do you have any suggestion coming from the field of Finite Volume?

Thank you for your interest.

Best Regards,

Antonella   

Dear Antonella,

I'll try to give you some answers:

1) periodical BC.s are setting by considering that all functions repeat themselves with period given by a characteristic lenght. The implementation depends on the discrete method. For example, spectral methods implicitly imply periodicit, finite difference requires setting periodicity in value at some nodes, finite volume at some face. Finite elements, that are based on shape functions, have periodicity implemented in the construction of the shape function. In a simple 1D example, consider nodes for 1 to N with periodicity implied by 1=N. The method works resolving unknown functions from 2 to N and setting 1=N after that. You can easily see how the shape functions are prescribed to fulfill the periodicity condition.

2) the type and numbers of BC.s depends on the mathematical character of the equations. Euler equations are strictly hyperbolic, Navier-Stokes are incompletely parabolic. That changes the setting of the BC.s. However, for very small viscosity, NS can be seen as perturbed Euler equations.

3) If you cover your computational domain by a triangular tesselation, you can work solving the equations directly with FE in the physical space, without using any mapping into the computational space...For example, n.Grad can be (2D) written:

n = nx i + ny j  and Grad = d/dx i + d/dy j

n.Grad (N1(x,y))= nx d N1/dx + ny d N1/dy

and so on for the other shape functions and 3D case. Such derivatives can be calculated analytically.

If you need other details do not hesitate to contact me.