Chart.js Example #1

\(\psi_{l} \times 10^7\):
\(\psi_{v} \times 10^7\):

Boundary Layer Mesh Calculator

I will implement some calculators to estimate the proper settings for the boundary prism layer meshing.

Simple mathematical calculations, such as the four arithmetic operations and power, can be done using HTML and JavaScript.

Boundary Layer Mesh

l_{n} = l_{1} r^{n-1}
l_{tot} = \sum_{k=1}^{n} l_{1} r^{k-1} = \frac{l_{1} \left( r^n – 1 \right)}{r – 1}    \left( r \neq 1 \right)

First layer thickness \(l_{1}\):
# Prism layers \(n\):
Growth rate \(r\):
Write precision:  6


Last layer thickness \(l_{n}\):
Total thickness \(l_{tot}\):

Boundary Layer Thickness – Internal Flow

Under construction – please check again later

Boundary Layer Thickness – External Flow

Under construction – please check again later

Matplotlib – How to Plot Mathematical expressions (2D)

What We Want to Do

We want to plot the following three mathematical expressions in 2D that describe the speed of sound \(c\) as a function of a vapor fraction \(\gamma\).

c = \frac{1}{\sqrt{\psi}}. \tag{1} \label{eq:psi-sos}

    • linear

\psi_{linear} = \gamma \psi_v + \left(1 – \gamma \right)\psi_l \tag{2} \label{eq:linear}

    • Wallis

\psi_{Wallis} = \left( \gamma \rho_{v, Sat} + \left( 1 – \gamma \right)\rho_{l, Sat} \right)
\left( \gamma \frac{\psi_v}{\rho_{v, Sat}} + \left( 1 – \gamma \right)\frac{\psi_l}{\rho_{l, Sat}} \right) \tag{3} \label{eq:wallis}

    • Chung

\psi_{Chung} = \left( \left( \frac{1 – \gamma}{\sqrt{\psi_v}} + \gamma \frac{s_{fa}}{\sqrt{\psi_l}} \right)\frac{\sqrt{\psi_v \psi_l}}{s_{fa}} \right)^2 \tag{4} \label{eq:chung}
s_{fa} = \sqrt{ \frac{\frac{\rho_{v, Sat}}{\psi_v}}{\left( 1- \gamma \right)\frac{\rho_{v, Sat}}{\psi_v} + \gamma \frac{\rho_{l, Sat}}{\psi_l}} } \tag{5} \label{eq:chung-sfa}

For more details of the above expressions, you might want to read the following post.

Sample Code

Free Access to A Course of Theoretical Physics by L.D. Landau & E.M. Lifshitz

The famous physics textbook written by the theoretical physicist L. D. Landau and E. M. Lifshitz are freely accessible from the following links.

Mechanics 3rd Edition (Course of Theoretical Physics Vol. 1)

The Classical Theory of Fields (Course of Theoretical Physics Vol. 2)

Quantum Mechanics (Course of Theoretical Physics Vol. 3)

Relativistic Quantum Physics Part1 (Course of Theoretical Physics Vol. 4)

Statistical Physics (Course of Theoretical Physics Vol. 5)

Fluid Mechanics (Course of Theoretical Physics Vol. 6)

Theory of Elasticity (Course of Theoretical Physics Vol. 7)

Electrodynamics of Continuous Media (Course of Theoretical Physics Vol. 8)

cavitatingFoam – barotropicCompressibilityModel (v1812)

Last Updated: June 9, 2019

The cavitatingFoam is a transient cavitation solver based on the homogeneous equilibrium model (HEM) from which the compressibility of the liquid/vapour ‘mixture’ is obtained, whose density varies from liquid density to vapor one according to the chosen barotropic equation of state.

Compressibility and speed of sound

\frac{D \rho}{Dt} = \psi \frac{Dp}{Dt} \tag{1} \label{eq:eos}

\psi = \rho \beta_s \tag{2} \label{eq:psi}

  • \(\rho\): density
  • \(p\): pressure
  • \(\beta_s\): isentropic compressibility

The definition of the isentropic compressibility can be deformed in the following way:
\beta_s &\equiv -\frac{1}{V} \left( \frac{\partial V}{\partial p} \right)_s \tag{3} \label{eq:betas} \\
&= -\frac{1}{V} \left(\frac{\partial V}{\partial \rho} \frac{\partial \rho}{\partial p} \right)_s \\
&= \rho \frac{1}{\rho^2} \left(\frac{\partial \rho}{\partial p} \right)_s \\
&= \frac{1}{\rho c^2}
where \(c\) is the speed of sound. We can substitute this expression for \(\beta_s\) in \eqref{eq:psi}, obtaining the following expression:
\psi = \frac{1}{c^2}. \tag{4} \label{eq:psi-sos}
The compressibility of the mixture equals to the reciprocal of the square of the speed of sound.

We can specify a compressibility function \(\psi\) that defines the coupling between the density and pressure in the thermophysicalProperties file.


The following three models for \(\psi\) are available in OpenFOAM v1812.

  • linear
  • Chung
  • Wallis

  • psiv: compressibility \(\psi\) of the vapor phase
  • psil: compressibility \(\psi\) of the liquid phase
  • pSat: saturation pressure
  • rhovSat: density of vapor at saturation point
  • rholSat: density of liquid at saturation point

The expressions of the speed of sound \(c\) in the liquid/vapour mixture are different according to the selected function.


\psi = \gamma \psi_v + \left(1 – \gamma \right)\psi_l \tag{5} \label{eq:linear}


s_{fa} = \sqrt{ \frac{\frac{\rho_{v, Sat}}{\psi_v}}{\left( 1- \gamma \right)\frac{\rho_{v, Sat}}{\psi_v} + \gamma \frac{\rho_{l, Sat}}{\psi_l}} } \tag{6} \label{eq:chung-sfa}
\psi = \left( \left( \frac{1 – \gamma}{\sqrt{\psi_v}} + \gamma \frac{s_{fa}}{\sqrt{\psi_l}} \right)\frac{\sqrt{\psi_v \psi_l}}{s_{fa}} \right)^2 \tag{7} \label{eq:chung}


\psi = \left( \gamma \rho_{v, Sat} + \left( 1 – \gamma \right)\rho_{l, Sat} \right)
\left( \gamma \frac{\psi_v}{\rho_{v, Sat}} + \left( 1 – \gamma \right)\frac{\psi_l}{\rho_{l, Sat}} \right) \tag{8} \label{eq:wallis}
\frac{\psi}{\gamma \rho_{v, Sat} + \left( 1 – \gamma \right)\rho_{l, Sat}} = \gamma \frac{\psi_v}{\rho_{v, Sat}} + \left( 1 – \gamma \right)\frac{\psi_l}{\rho_{l, Sat}} \tag{8′} \label{eq:wallis2}


[1] Christopher Earls Brennen, Fundamentals of Multiphase Flow. Cambridge University Press (2005)

twoPhaseEulerFoam – Diameter Model (v1812)


Constant dispersed-phase particle diameter model.


Isothermal dispersed-phase particle diameter model.
p \times \frac{4}{3} \pi \left( \frac{d}{2} \right)^3 &= p_0 \times \frac{4}{3} \pi \left( \frac{d_0}{2} \right)^3 \\
d &= d_0 \left( \frac{p_0}{p} \right)^{1/3} \tag{1} \label{eq:levmdiffusion}


IATE (Interfacial Area Transport Equation) bubble diameter model.
Solves for the interfacial curvature per unit volume of the phase rather than interfacial area per unit volume to avoid stability issues relating to the consistency requirements between the phase fraction and interfacial area per unit volume.

The following bubble interaction mechanisms are considered in this model.

    • wakeEntrainmentCoalescence
      Bubble coalescence by wake entrainment
    • randomCoalescence
      Bubble coalescence by random collision
    • turbulentBreakUp
      Bubble break-up due to the impact of turbulent eddies
      \phi_{TI} = \frac{C_{ti}}{3} \frac{U_t}{d} \sqrt{1 – \frac{We_{Cr}}{We}} e^{-\frac{We_{Cr}}{We}}

      • \(We\): Weber number [-]
      • \(U_t\): Turbulent fluctuation velocity [m/s]

twoPhaseEulerFoam – phaseProperties

In solving two-phase flow problems, we can use a variety of simulation techniques, such as Eulerian-Eulerian, Eulerian-Lagrangian methods, Volume of Fluid (VOF) method and homogeneous equilibrium model (HEM).

In this page, I will organize information on the physical models available in twoPhaseEulerFoam, an Eulerian-Eulerian solver available in OpenFOAM.

I will add descriptions of each model.

Selected OpenFOAM version:

drag model
  • Ergun
  • Gibilaro
  • GidaspowErgunWenYu
  • GidaspowSchillerNaumann
  • IshiiZuber
  • Lain
  • SchillerNaumann
  • SyamlalOBrien
  • TomiyamaAnalytic
  • TomiyamaCorrelated
  • WenYu
  • segregated

Read more

lift model
  • LegendreMagnaudet
  • Moraga
  • Tomiyama
  • constantCoefficient
  • none
heatTransfer model
  • RanzMarshall
  • spherical
turbulentDispersion model
  • Burns
  • Gosman
  • LopezDeBertodano
  • constantCoefficient
  • none
virtualMass model
  • Lamb
  • constantCoefficient
  • none
aspectRatio model
  • Tomiyama
  • VakhrushevEfremov
  • Wellek
  • constant
wallLubrication model
  • Antal
  • Frank
  • Tomiyama
  • none
diameter model
  • constant
  • isothermal
  • IATE

Read more

Non-Reflecting Boundary Conditions in OpenFOAM

What we want to achieve

When we simulate fluid flow, we have to cut a finite computational domain out of an entire flow region. For accurate simulation, we need to let fluid and sound wave flow smoothly out of the domain through the boundary.

The reflection at the boundary has a larger effect on the solution especially when we perform a compressible flow simulation as can be seen in the following movie (Upper: With reflection, Lower: Without reflection of sound wave).

In OpenFOAM, we can use two approximate non-reflecting boundary conditions:

They determine the boundary value by solving the following equation

\frac{D \phi}{D t} = \frac{\partial \phi}{\partial t} + \boldsymbol{U} \cdot \nabla \phi = 0, \tag{1} \label{eq:advection}

where \(D/Dt\) is the material derivative and \(\boldsymbol{U}(\boldsymbol{x}, t)\) is the advection velocity.

We assume that the advection velocity \(\boldsymbol{U}\) is parallel to the boundary (face) normal direction and rewrite the eqn. \eqref{eq:advection} as

\frac{D \phi}{Dt} \approx \frac{\partial \phi}{\partial t} + U_{n} \cdot \frac{\partial \phi}{\partial \boldsymbol{n}}= 0, \tag{2} \label{eq:advection2}

where \(\boldsymbol{n}\) is the outward-pointing unit normal vector.

These boundary conditions are different in how the advection speed (scalar quantity) \(U_{n}\) is calculated and it is calculated in advectionSpeed() member function.

advective B.C.

The advection speed is the component of the velocity normal to the boundary

U_n = u_n. \tag{3} \label{eq:advectiveUn}

waveTransmissive B.C.

The advection speed is the sum of the component of the velocity normal to the boundary and the speed of sound \(c\)

U_n = u_n + c = u_n + \sqrt{\gamma/\psi}, \tag{4} \label{eq:waveTransmissiveUn}

where \(\gamma\) is the ratio of specific heats \(C_p/C_v\) and \(\psi\) is compressibility.

What do lInf and fieldInf mean?

Coming soon.

Numerical Libraries – Linear Algebra

This page will be sequentially updated.

The linear algebraic operations, such as solving systems of linear equations and matrix multiplication, often appear in the field of machine learning, computational fluid dynamics, computer graphics, and so on. As the dimensions of systems or matrices will be bigger, these operations will become more time consuming and complicated.

There are many optimized open-source numerical libraries for linear algebra, such as LAPACK and Eigen, so we can use them to get better performance even when our matrices are very large.

For example, an open source CFD software SU2 uses BLAS and LAPACK, and moreover the performance of NumPy can be improved with the use of these optimized linear algebra libraries.

In this page, I try to succinctly summarize what we need to know at a minimum about open-source numerical libraries.

BLAS (Basic Linear Algebra Subprograms)
  • [URL]
  • BLAS libraries provide standard building blocks for performing basic vector and matrix operations
    • Level 1: scalar, vector and vector-vector operations
    • Level 2: matrix-vector operations
    • Level 3: matrix-matrix operations
BLACS (Basic Linear Algebra Communication Subprograms)
LAPACK (Linear Algebra PACKage)
ScaLAPACK (Scalable Linear Algebra PACKage)