Hello,

I would like to simulate a collapsing water wolumn by using my own Eulerian FEM code (incopressible viscous).

Problem: I have problems calculating the global pressure of the system. Since I need the inverse of the nodal mass, which is "zero" at the nodes of void elements.

Does somebody know how to deal with the mass of void elements?

I tried to used 1.0d-6 as void mass but it leads to huge numerial problems.

Thank you in advance for any help!

Boni

Hi Boni

for two phase flow problems with high density ratio between the phases you must have a jump in the properties (density and viscosity) of the fluids near the interface, but the mass of your void fluid will not be completely "zero".

Are you computing the density and viscosity for both phases correctly?

In the Navier Stokes equation you'll have to compute these quantities taking the following relationship

rho = phi*rho_A + (1-phi)*rho_B mu = phi* mu_A + (1-phi)* mu_B

considering phi the volume/concentration "marker" = {1 for fluid A and 0 for fluid B)

Cheers

Renato.

Hallo Renato

>for two phase flow problems with high density ratio between >the phases you must have a jump in the properties (density >and viscosity) of the fluids near the interface, but the mass >of your void fluid will not be completely "zero".

at interface you are right but beyond the interface the density is zero. My Eulerian model looks like this:

#################################### #################################### #################################### #################################### @@@@@@@@@################### @@@@@@@@@################### @@@@@@@@@################### @@@@@@@@@################### @@@@@@@@@################### @@@@@@@@@################### @@@@@@@@@################### @@@@@@@@@###################

@ for water(elements) # for void(elements)

As for incompressible fluids, I keed the density constant. The same for the viscosity. I used a fractional step method for decopling velocity and pressure. The pressure is then solved by the pressure-Poisson equation: [C^T*M^-1*C]p = C^T*a

I get numerical problems since M^-1 (inverse of lumped mass matrix) has the huge numbers (1/m) for nodes of void element...

Thanks for helping

How does your air could have a zero density? I can't see.

I've just simulated this problem in my FEM code too, for this purpose I've used the following fluid properties:

air: rho_# = 1 Kg/m^3; mu_# = 0.0001 Kg/ms water: rho_@ = 1000 Kg/m^3; mu_@ = 0.01 Kg/ms

now, suppose you're computing the density in the air region, so:

rho = phi*rho_@ + (1-phi)*rho_# rho = 0.0 * 1000 + (1.0-0.0) * 0.0001 = 0.0001 Kg/ms

Regards

Renato.

-------------------------------------------- NOTE: NICE PICTURE YOU DID )

Thanks Renato for constructive postings!

I see, you gaved the void elements a material (air). I did not. For explicit dynamic (solid) problems I just avoid the calculation for element without material (filling rate < tol)...

OK, let's assume a 2-phase simulation. Do you use the same (incompressible) solver for both water and air? means, do you consider air as incompressible?

Thank you once more. Regards Boni

After you said you were working with VOF (Eulerian approach) to solve a dam break problem I understood a two phase immiscible flow.

> Do you use the same (incompressible) solver for both water and air? means, do you consider air as incompressible? <

Yes, I consider the air phase incompressible and use the same solver to compute both fluids flow at the same time. in fact, the air has little influence in this case and you can even avoid the air calculation doing some kind of pressure and velocity projection on elements near the interface. I don't do it because I'd rather keep my solver working for problems with and without high property ratios.

Renato.

I see. It sounds pretty reasonable the small effect of air, so I will try with the same code too.

>air: rho_# = 1 Kg/m^3; mu_# = 0.0001 Kg/ms water: rho_@ = 1000 Kg/m^3; mu_@ = 0.01 Kg/ms

Is there any reason why you use above (dynamic) viscosities? They are on factor 10 bigger. From Wikipedia I have that mu_# is around 0.00001 Kg/ms and mu_@ is around 0.001 Kg/ms (Pa.s) Or was it a typing mistake?

Thanks for great hints.

The properties I cited were got from the first article I found in a Google searching (A. Kaceniauskas, "Dam break flow simulation by the pseudo-concentration method", Mechanika - 2005 (151)1:39-43), but after reading some other papers, I realized that since the fluid flow is mainly driven by the water column weight, the viscous forces in the air region will also have little influence in the problem.

Renato.

I created a new material "Air" with properties posted here.

The Courant condition however gives very small time steps (1.d-8).

dt = CFL * minlen / (maxvel+minC) - CFL = Courant number = 0.01 - minlen is the smallest element edge - maxvel is the maximal velocity - minC is the smallest speed of sound (here Air has 340m/s)

Do you have same time steps? The model size is (I hope the same with yours) 0.6m x 0.4m. The water column is 0.15m x 0.3m

Thanks and sorry for taking your time. Boni

I followed the dimensions described in the article of the previous post which followed the dimension of the famous Koshizuka's experiment (Koshizuka S., Tamako H., Oka Y, "A particle method for incompressible viscous flow with fluid fragmentation", Comput. Fluid Dynamics, 1995 (4) 1:29-46)

A = 0.146 m (characteristic length) Box = 2.4*A x 4*A x A/2 water column = 2*A x A x A/2

Regarding the time step, at the first time I run this problem I used a fixed time step of 0.01 and the results were good up to the leading edge reached the opposite wall. After the water stricken on the opposite wall the Reynolds number (in the water phase) tends to increase quickly. It was confirmed after running this problem with adapted time step. However, my shortest time steps were never less than 10^-3. Off course, it'll depend on your mesh size.

I saw you're dealing with CFL number in the air region as it was a compressible flow, why? You don't need to worry about your air region, it doesn't matter in this problem. This problem is mainly driven by the water region weight. Take in mind you're working with a water flow where the air region was inserted in your problem only for completeness and is used only to help when marking the interface between the fluids.

Regards

Renato.

Hi Renato

0.01 seconds? Do you solved the NSE (momentum eq.) by using explicit or implicit time integration? And FEM or FDM?

Actually I am trying to extend my explicit FEM code for solid dynamics to fluid dynamics. As you know explicit codes need smaller time steps than implicit ones. Therefore I use the CFL condition each time step.

BTW, exactly Koshizuka's experiment I trying to simulate

Regards

Boni

As I said, 0.01 was a first attempt, but the code was able to solve the problem and the results were good.

Basically, my code runs an edge based SUPG/PSPG FEM implicit solver in meshes with linear tetrahedra where the NS equation is fully coupled solved (momentum and mass conservations are solved in only one system of equations).

There's a short description of my steady flow solver in the following articles:

R. N. Elias, M. A. D. Martins, A. L. G. A. Coutinho, Parallel Edge-Based Inexact Newton Solution of Steady Incompressible 3D Navier-Stokes Equations, Euro-Par 2 sept 2005, Lisbon, Portugal

R. N. Elias, A. L. G. A. Coutinho, M. A. D. Martins Parallel Edge-Based Solution of Viscoplastic Flows with the SUPG/PSPG Formulation , FEF05 IACM Special Interest Conference supported by ECCOMAS, April 4-6, 2005, Swansea, Wales, UK - published in Computational Mechanics 2005

R. N. Elias, A. L. G. A. Coutinho, M. A. D. Martins, R. M. Sydenstricker, Edge-based Data Structures for the Parallel Solution of 3D Steady Incompressible Fluid Flow with the SUPG/PSPG Finite Element Formulation , XXVI CILAMCE - Iberian Latin American Congress on Computational Methods - Out/2005, Guarapari, ES, Brasil

these articles are available to download in http://www.nacad.ufrj.br/~rnelias/pubs.html

Hope to be helpful.

Renato.

wow, highly interesting publications! Specially the combination FEM + implicit + parallel computing.

Implicit SUPG/PSPG FEM allows for sure those small time steps. I use explicit time integration assuming a lumped mass matrix. To keep the calculation stable, smaller time steps are required.

At first I will use 4-node quadrilateral Q1P0 elements. Different from your code, I use a fractional step (or projection) method for decoupling velocity and pressure.

Later I plan to do solid-fluid interaction with one code by using a mixture theory for handling multi-material elements.

Renato, very much obligado for very helpful postings!!

I will come to this thread back if I still have the problem. Boni

[latdatngoc]

I am newbie CFD, I am trying to solve the unsteady NS-using SUPG/PSPG scheme employ the one-step method. I use the cavity lid driven as a test case. But I have a trouble in treating the Dirichlet conditions. Does anybody can tell me some experiences for this?

Thanks in advance
Jammy