Conservation of Mass

The finite element formulation of the Lagrangian form of the mass conservation equation is given by:

d ρ d t | x = − ( ρ / V ) d V d t | x MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aadaWcaaqaaiaadsgacqaHbpGCaeaacaWGKbGaamiDaaaadaabbaqa aiaadIhaaiaawEa7aiabg2da9iabgkHiTmaabmaabaGaeqyWdiNaai 4laiaadAfaaiaawIcacaGLPaaadaWcaaqaaiaadsgacaWGwbaabaGa amizaiaadshaaaWaaqqaaeaacaWG4baacaGLhWoaaaa@4D14@

When transformed into the ALE formulation it gives:

Applying a Galerkin variation form for the solution of Equation 2:

Where, ψ is the Weighting function.

Using a finite volume formulation

Where, ψ =1

ρ = constant density over control volume V MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamOvaaaa@36D1@ .

Therefore:

∫ V ∂ ρ ∂ t d V + ∫ V ρ ∂ ν K ∂ x K d V = 0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aadaWdrbqaamaalaaabaGaeyOaIyRaeqyWdihabaGaeyOaIyRaamiD aaaacaWGKbGaamOvaiabgUcaRmaapefabaGaeqyWdihaleaacaWGwb aabeqdcqGHRiI8aaWcbaGaamOvaaqab0Gaey4kIipakmaalaaabaGa eyOaIyRaeqyVd42aaSbaaSqaaiaadUeaaeqaaaGcbaGaeyOaIyRaam iEamaaBaaaleaacaWGlbaabeaaaaGccaWGKbGaamOvaiabg2da9iaa icdaaaa@5447@

Using the divergence theorem leads to:

∫ V ∂ ρ ∂ t d V + ∫ S ρ ( ν ⋅ j     n j ) d S = 0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aadaWdrbqaamaalaaabaGaeyOaIyRaeqyWdihabaGaeyOaIyRaamiD aaaacaWGKbGaamOvaiabgUcaRmaapefabaGaeqyWdihaleaacaWGtb aabeqdcqGHRiI8aaWcbaGaamOvaaqab0Gaey4kIipakmaabmaabaGa eqyVd42aaSraaSqaaiaadQgaaeqaaOGaeyyXICTaaGPaVlaayIW7ca WGUbWaaSbaaSqaaiaadQgaaeqaaaGccaGLOaGaayzkaaGaamizaiaa dofacqGH9aqpcaaIWaaaaa@5889@

Further expansion gives:

d d t ∫ v ρ d V = ∫ s ρ ( w j − v j ) n j d S ︸ m a s s   f l u x MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aadaWcaaqaaiaadsgaaeaacaWGKbGaamiDaaaadaWdrbqaaiabeg8a YjaadsgacaWGwbGaeyypa0Zaa8quaeaadaagaaqaaiabeg8aYnaabm aabaGaam4DamaaBaaaleaacaWGQbaabeaakiabgkHiTiaadAhadaWg aaWcbaGaamOAaaqabaaakiaawIcacaGLPaaacaWGUbWaaSbaaSqaai aadQgaaeqaaOGaamizaiaadofaaSqaaiaad2gacaWGHbGaam4Caiaa dohacaaMc8UaamOzaiaadYgacaWG1bGaamiEaaGccaGL44paaSqaai aadohaaeqaniabgUIiYdaaleaacaWG2baabeqdcqGHRiI8aaaa@5E39@

This formula is still valid if density ρ is not assumed uniform over volume V MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamOvaaaa@36D1@ .
Figure 1. Mass Flux Across a Surface


The density, ρ i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaeqyWdi3aaS baaSqaaiaadMgaaeqaaaaa@38D0@ , is given computed:

ρ i = 1 2 ρ I { 1 + η   s i g n ( ϕ i ) } + 1 2 ρ J { 1 − η   s i g n ( ϕ i ) }

Where, 0 ≤ η ≤ 1 is the upwind coefficient given on the input card.

If η =0, there is no upwind.

Therefore, ρ i = ρ I + ρ J 2 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaeqyWdi3aaS baaSqaaiaadMgaaeqaaOGaeyypa0ZaaSaaaeaacqaHbpGCdaWgaaWc baGaamysaaqabaGccqGHRaWkcqaHbpGCdaWgaaWcbaGaamOsaaqaba aakeaacaaIYaaaaaaa@4117@ .

If η =1, there is full upwind.

The smaller the upwind factor, the faster the solution; however, the solution is more stable with a large upwind factor. This upwind coefficient can be tuned with parameter from /UPWIND keyword (not recommended, this keyword has been obsolete as of version 2018).

For a free surface: ρ J MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaeqyWdi3aaS baaSqaaiaadMgaaeqaaaaa@38D0@ = ρ I MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaeqyWdi3aaS baaSqaaiaadMgaaeqaaaaa@38D0@