From fc02e4c5f8c8bbb014449dc27d7b69992ad6043f Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Tue, 16 Apr 2019 22:45:11 +0200
Subject: Add basic physical scaling and Knudsen quality criterion
The paper "Automatic grid refinement criterion for lattice Boltzmann method" [2015]
by Lagrava et al. describes a criterion for measuring the local simulation quality
using a comparison of the theoretical Knudsen number and the quotient of the cells's
non-equilibrium and equilibrium function.
While this criterion was developed to enable automatic selection of areas to be refined,
it also offers a interesting and unique perspective on the fluid structure.
As the criterion requires calculation of the modeled Reynolds-, Mach- and Knudsen-numbers
I took the time to set up the basics for scaling the simulation to actually model a physical
system. Or rather calculating which physical model is represented by the chosen resolution
and relaxation time.
[2015]: https://arxiv.org/abs/1507.06767
---
src/shader/code/collide.glsl | 50 ++++++++++++++++++++++++++++++++++++++------
src/shader/code/vertex.glsl | 17 +++++++++++----
2 files changed, 57 insertions(+), 10 deletions(-)
(limited to 'src')
diff --git a/src/shader/code/collide.glsl b/src/shader/code/collide.glsl
index 3bfd332..9ff6ab3 100644
--- a/src/shader/code/collide.glsl
+++ b/src/shader/code/collide.glsl
@@ -11,6 +11,12 @@ uniform uint nX;
uniform uint nY;
uniform uint iT;
+/// Fluid characteristics
+
+const float physCharLength = 1.0;
+const float physCharVelocity = 1.0;
+const float physViscosity = 0.1;
+
/// LBM constants
const uint q = 9;
@@ -19,9 +25,25 @@ const float weight[q] = float[](
1./9. , 4./9., 1./9. ,
1./36 , 1./9., 1./36.
);
+const float invCs2 = 1./3.;
+
+const float relaxationTime = 0.6;
+const float relaxationFrequency = 1 / relaxationTime;
+
+/// Unit conversion
+
+const float convLength = physCharLength / nX;
+const float convTime = (relaxationTime - 0.5) / invCs2 * convLength*convLength / physViscosity;
+const float convVelocity = convLength / convTime;
+const float convViscosity = convLength * convLength / convTime;
+
+const float latticeCharVelocity = physCharVelocity / convVelocity;
+
+/// Emergent fluid numbers
-const float tau = 0.65;
-const float omega = 1/tau;
+const float Re = physCharVelocity * physCharLength / physViscosity;
+const float Ma = latticeCharVelocity * sqrt(invCs2);
+const float Kn = Ma / Re;
/// Vector utilities
@@ -65,12 +87,18 @@ void set(uint x, uint y, int i, int j, float v) {
streamCells[indexOfLatticeCell(x,y) + indexOfDirection(i,j)] = v;
}
-void setFluid(uint x, uint y, vec2 v) {
+void setFluidVelocity(uint x, uint y, vec2 v) {
const uint idx = indexOfFluidVertex(x, y);
fluidCells[idx + 0] = v.x;
fluidCells[idx + 1] = v.y;
}
+void setFluidQuality(uint x, uint y, float knudsen, int quality) {
+ const uint idx = indexOfFluidVertex(x, y);
+ fluidCells[idx + 0] = knudsen;
+ fluidCells[idx + 1] = quality;
+}
+
int getMaterial(uint x, uint y) {
const uint idx = indexOfFluidVertex(x, y);
return int(fluidCells[idx + 2]);
@@ -149,19 +177,29 @@ void main() {
d = 1.0;
}
- setFluid(x,y,v);
+ //setFluidVelocity(x,y,v);
+
+ float knudsen = 0.0;
for ( int i = -1; i <= 1; ++i ) {
for ( int j = -1; j <= 1; ++j ) {
- set(x+i,y+j,i,j, get(x,y,i,j) + omega * (equilibrium(d,v,i,j) - get(x,y,i,j)));
+ const float feq = equilibrium(d,v,i,j);
+ const float fneq = get(x,y,i,j) - feq;
+ knudsen += abs(fneq / feq);
+
+ set(x+i,y+j,i,j, get(x,y,i,j) + relaxationFrequency * (feq - get(x,y,i,j)));
}
}
+
+ knudsen /= q;
+
+ setFluidQuality(x,y,knudsen,int(round(log2(knudsen / Kn))));
}
if ( isBounceBackCell(material) ) {
for ( int i = -1; i <= 1; ++i ) {
for ( int j = -1; j <= 1; ++j ) {
- set(x+(-1)*i,y+(-1)*j,(-1)*i,(-1)*j, get(x,y,i,j) + omega * (equilibrium(d,v,i,j) - get(x,y,i,j)));
+ set(x+(-1)*i,y+(-1)*j,(-1)*i,(-1)*j, get(x,y,i,j) + relaxationFrequency * (equilibrium(d,v,i,j) - get(x,y,i,j)));
}
}
}
diff --git a/src/shader/code/vertex.glsl b/src/shader/code/vertex.glsl
index 676b68e..7851292 100644
--- a/src/shader/code/vertex.glsl
+++ b/src/shader/code/vertex.glsl
@@ -10,7 +10,8 @@ out VS_OUT {
uniform uint nX;
uniform uint nY;
-const float displayAmplifier = 3.0;
+const float velocityDisplayAmplifier = 3.0;
+const int qualityDisplayRestrictor = 6;
float unit(float x) {
return 1.0/(1.0+exp(-x));
@@ -40,6 +41,14 @@ bool isWallFrontier(int material) {
return material == 2 || material == 3;
}
+float restrictedQuality(float quality) {
+ if ( quality < 0.0 ) {
+ return 0.0;
+ } else {
+ return min(1.0, quality / qualityDisplayRestrictor);
+ }
+}
+
void main() {
const vec2 idx = fluidVertexAtIndex(gl_VertexID);
@@ -58,9 +67,9 @@ void main() {
vs_out.color = vec3(0.0, 0.0, 0.0);
} else {
vs_out.color = mix(
- vec3(-0.5, 0.0, 1.0),
- vec3( 1.0, 0.0, 0.0),
- displayAmplifier * norm(VertexPosition.xy)
+ vec3(0.0, 1.0, 0.0),
+ vec3(1.0, 0.0, 0.0),
+ restrictedQuality(VertexPosition.y)
);
}
}
--
cgit v1.2.3