The curious engineer is still thriving in me - checking the veracity of the continuity equation in computational fluid dynamics...

 It was early in the morning. Everywhere there was a pin-drop silence. After the regular morning exercise, I turned on my laptop. Today's plan was to verify the continuity equation of computational fluid dynamics. For noncompressible fluid, the continuity equation can be mathematically written as ∇⋅U=0. It's the conservation of mass in fluid simulation. Whatever goes into a volume has to come out - no creation, no destruction, just pure, seamless flow.

The paraFoam is good for visualizing the cavity example of the Tutorial that comes with openFoam. As I am new to using ParaView, i preferred raw data exported to a CSV file.

So I created the following Python script.

from paraview.simple import *

# Load the OpenFOAM case

foamCase = OpenFOAMReader(FileName="cavity.OpenFOAM")

foamCase.MeshRegions = ['internalMesh']

foamCase.CellArrays = ['U']  # Ensure velocity field is loaded


# Apply the new Gradient filter

gradFilter = Gradient(Input=foamCase)

gradFilter.ScalarArray = ['POINTS', 'U']  # Compute gradient of velocity U


# Show the gradient field in ParaView

gradDisplay = Show(gradFilter)

ColorBy(gradDisplay, ('POINTS', 'Gradient'))


# Apply a color map for better visualization

glyph = Glyph(Input=gradFilter, GlyphType='Arrow')

glyph.ScaleFactor = 0.01  # Adjust the size

Show(glyph)

Render()

# Optional: Save visualization as an image

SaveScreenshot("grad_U_visualization.png")

writer = CreateWriter("gradU_data.csv", gradFilter)

writer.UpdatePipeline()

Now the time is to test the veracity of the continuity equation.

So, another Python script as follows.

import pandas as pd
import matplotlib.pyplot as plt

# Load CSV file
df = pd.read_csv("/home/som/OpenFOAM/som-dev/run/cavity/cavity/gradU_data.csv")

# Compute divergence of U (sum of diagonal terms of gradient matrix)
df["divU"] = df["Gradient:0"] + df["Gradient:4"] + df["Gradient:8"]

# Plot divergence
plt.plot(df["Points:0"], df["divU"], marker="o", linestyle="-", label="∇·U")
plt.axhline(0, color='r', linestyle="--", label="Zero Line")
plt.xlabel("X Position")
plt.ylabel("Divergence of Velocity (∇·U)")
plt.title("Verification of Continuity Equation (∇·U)")
plt.legend()
plt.grid()
plt.show()

# Check if max deviation from zero is small
print("Max |∇·U| =", df["divU"].abs().max())

A little explanation for the code 

df["divU"] = df["Gradient:0"] + df["Gradient:4"] + df["Gradient:8"]

This is important because we must take the diagonal values of a 

3X3 matrix to verify the continuity equation.

Everything was looking fine except near the lid. The value was high. 

Max |∇·U| = 100.0.

The culprits were located near the moving lid and the corner regions of the cavity. This made sense—corners were prone to numerical singularities due to discontinuities in the velocity gradient.

Not to my surprise. But i wanted to verify using a histogram image.

So added the following python script.

plt.hist(df["divU"], bins=50)
plt.xlabel("∇·U values")
plt.ylabel("Frequency")
plt.title("Distribution of divergence (∇·U)")
plt.show()

And voila.

Most of the values are almost 0.

The scientist in me became happy to contribute to the learning community.

Happy learning...

Here's the simulation exported from ParaView...

That's it for today.

Reclaiming #WhoWeAre...