diff options
author | Adrian Kummerlaender | 2019-04-16 22:45:11 +0200 |
---|---|---|
committer | Adrian Kummerlaender | 2019-04-16 22:45:11 +0200 |
commit | fc02e4c5f8c8bbb014449dc27d7b69992ad6043f (patch) | |
tree | af962916b7b3239c698a4825ba24ea1ddc3377ce /src | |
parent | 89718b245b0246e08f9d0545a636e0cf89642a72 (diff) | |
download | compustream-fc02e4c5f8c8bbb014449dc27d7b69992ad6043f.tar compustream-fc02e4c5f8c8bbb014449dc27d7b69992ad6043f.tar.gz compustream-fc02e4c5f8c8bbb014449dc27d7b69992ad6043f.tar.bz2 compustream-fc02e4c5f8c8bbb014449dc27d7b69992ad6043f.tar.lz compustream-fc02e4c5f8c8bbb014449dc27d7b69992ad6043f.tar.xz compustream-fc02e4c5f8c8bbb014449dc27d7b69992ad6043f.tar.zst compustream-fc02e4c5f8c8bbb014449dc27d7b69992ad6043f.zip |
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
Diffstat (limited to 'src')
-rw-r--r-- | src/shader/code/collide.glsl | 50 | ||||
-rw-r--r-- | src/shader/code/vertex.glsl | 17 |
2 files changed, 57 insertions, 10 deletions
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) ); } } |