Initial commit

.gitattributes

@@ -0,0 +1,2 @@
+# Auto detect text files and perform LF normalization
+* text=auto

Hyberbolic Tangent Profile.m

@@ -0,0 +1,63 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%   MATLAB code to simulate 1D advection-diffusion for a hyperbolic tangent profile
+%   Assumed : 
+%   Delta x (lattice distance) = Delta t (lattice time step) = 1 
+%   c = 1, c_s (speed of sound) = 1/sqrt(3)
+%   Periodic boundary conditions are applied at the corner grid points
+%   Author : Sthavishtha Bhopalam Rajakumar
+%   Updated date : 30-09-2017
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clc;clear
+
+%defined parameters in the question
+d = 5e-8; %Diffusion coefficient
+delta = 0.01; %parameter in the hyperbolic tangent equation
+u_max = 0.1; %macroscopic fluid velocity
+n = 2000; %number of grids
+t_max = 4000; %maximum time/iteration to be reached
+c_s = 1./sqrt(3); %speed of sound
+ma = u_max./c_s; %mach number
+tau = d./(c_s*c_s); %relaxation time
+beta = 1./(2*tau+1); %factor considering the effect of relaxation time
+
+%initializing initial populations
+rho = zeros(n,1); %particle densities
+for i=1:n
+    rho(i) = 1.0+0.5*(1 - tanh(((i-1)/(n-1) - 0.2)/delta));
+end
+feq = zeros(3,n); %equilibrium populations
+feq(1,:) = 2*rho(:).*(2 - sqrt(1 + ma*ma))./3; %defining equilibrium populations
+feq(2,:) = rho(:).*((u_max - c_s*c_s)/(2*c_s*c_s) + sqrt(1 + ma*ma))./3;
+feq(3,:) = rho(:).*((-u_max - c_s*c_s)/(2*c_s*c_s) + sqrt(1 + ma*ma))./3;
+f = feq;
+
+for t = 1:t_max     %t = time/iteration
+
+    if t~=1
+        rho(:) = f(1,:) + f(2,:) + f(3,:); %density kinetic equation
+        rho(1) = 2.0; %dirichlet boundary conditions - specified densities
+        rho(n) = 1;
+        feq(1,:) = 2*rho(:).*(2 - sqrt(1 + ma*ma))./3; %defining equilibrium populations
+        feq(2,:) = rho(:).*((u_max - c_s*c_s)/(2*c_s*c_s) + sqrt(1 + ma*ma))./3;
+        feq(3,:) = rho(:).*((-u_max - c_s*c_s)/(2*c_s*c_s) + sqrt(1 + ma*ma))./3;
+        f(:,1) = feq(:,1); f(:,n) = feq(:,n);%equilibirum populations at boundaries due to bcs
+        f(1,2:n-1) = f(1,2:n-1) - 2*beta*(f(1,2:n-1) - feq(1,2:n-1)); %performing collision 
+        f(2,2:n-1) = f(2,2:n-1) - 2*beta*(f(2,2:n-1) - feq(2,2:n-1));
+        f(3,2:n-1) = f(3,2:n-1) - 2*beta*(f(3,2:n-1) - feq(3,2:n-1));
+    end
+
+    %advection/streaming
+    for i=n:-1:2
+        f(2,i) = f(2,i-1);
+    end
+
+    for i=1:n-1
+        f(3,i) = f(3,i+1);
+    end
+end
+%post-processing
+x = 1:n;
+rho_act(x) = 1.0+0.5*(1 - tanh(((x-1)/(n-1) - 0.2)/delta));
+plot(x,rho_act(x),'k--',x,rho(x),'k-');
+

README.md

@@ -0,0 +1,8 @@
+**LBM 1D**
+
+This repository hosts the exercise solutions of the course 151-0213-00L [**Fluid Dynamics with the Lattice Boltzmann Method**](http://www.vvz.ethz.ch/lerneinheitPre.do?semkez=2017W&lerneinheitId=116549&lang=en) during the Autumn semester 2017 at ETH Zurich.
+
+**Contents**
+- Steep gaussian profile (1D Advection-diffusion)
+- Hyperbolic tangent profile (1D Advection-diffusion)
+- Shock tube (N-S 1D)

Shock Tube simulation.m

@@ -0,0 +1,77 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%   MATLAB code to simulate 1D NSE in a shock tube
+%   Assumed : 
+%   Delta x (lattice distance) = Delta t (lattice time step) = 1 
+%   c = 1, c_s (speed of sound) = 1/sqrt(3)
+%   Periodic boundary conditions are applied at the corner grid points
+%   Author : Sthavishtha Bhopalam Rajakumar
+%   Updated date : 30-09-2017
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clc;clear
+
+%defined parameters in the question
+nu = 0.06; %Diffusion coefficient
+% u_max = 0; %macroscopic fluid velocity
+n = 800; %number of grids
+t_max = 500 ; %maximum time/iteration to be reached
+c_s = 1./sqrt(3); %speed of sound
+% ma= u_max./c_s; %mach number
+tau = nu./(c_s*c_s); %relaxation time
+beta = 1./(2*tau+1); %factor considering the effect of relaxation time
+c = [0 1 -1]; %direction velocities
+
+%initializing initial populations
+rho = zeros(n,1); %particle densities
+u = zeros(n,1); %particle velocities
+for i=1:n %initial densities
+    if i <= n/2
+        rho(i) = 1.5;
+    else
+        rho(i) = 1.0;
+    end
+end
+feq = zeros(3,n); %equilibrium populations
+feq(1,:) = 2*rho(:)/3; %defining equilibrium populations
+feq(2,:) = rho(:)./6;
+feq(3,:) = rho(:)./6;
+f = feq; %initial population is in equilbrium
+
+for t = 1:t_max   %t = time/iteration
+
+    if t~=1
+        rho(:) = f(1,:) + f(2,:) + f(3,:); %density kinetic equation
+        u(:) = (f(1,:)*c(1) + f(2,:)*c(2) + f(3,:)*c(3))./(rho(:))'; %momentum kinetic equation
+        ma = u./(c_s); %local mach number
+        feq(1,:) = 2*rho(:).*(1 - (u.^2/(2*c_s*c_s)))./3; %defining equilibrium populations - rewritten
+        feq(2,:) = rho(:).*((1 + (u(:))./(c_s*c_s) + u.^2/(c_s*c_s)))./6;
+        feq(3,:) = rho(:).*((1 - (u(:))./(c_s*c_s) + u.^2/(c_s*c_s)))./6;
+%         feq(1,:) = 2*rho(:).*(2 - sqrt(1 + ma.^2))./3; %defining equilibrium populations
+%         feq(2,:) = rho(:).*((-1./2 + u.*c(2)/(c_s*c_s) + sqrt(1 + ma.^2)))./3;
+%         feq(3,:) = rho(:).*((-1./2 - u.*c(2)/(c_s*c_s) + sqrt(1 + ma.^2)))./3;
+        f(1,:) = f(1,:) - 2*beta*(f(1,:) - feq(1,:)); %performing collision 
+        f(2,:) = f(2,:) - 2*beta*(f(2,:) - feq(2,:));
+        f(3,:) = f(3,:) - 2*beta*(f(3,:) - feq(3,:));
+    end
+
+    %advection across direction of c
+    for i=n:-1:2
+        f(2,i) = f(2,i-1);
+    end
+
+    for i=1:n-1
+        f(3,i) = f(3,i+1);
+    end
+
+    %bounce-back
+    f(2,1) = f(3,1);
+    f(3,n) = f(2,n);
+
+end
+% for i=1:n
+%     re = u.*i/nu;
+% end
+x = 1:n;
+plot(x,rho(x),'k-');
+
+

Steep Gaussian - Advection Diffusion.m

@@ -0,0 +1,66 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%   MATLAB code to simulate 1D advection-diffusion for a steep gaussian profile
+%   Assumed : 
+%   Delta x (lattice distance) = Delta t (lattice time step) = 1 
+%   c = 1, c_s (speed of sound) = 1/sqrt(3)
+%   Periodic boundary conditions are applied at the corner grid points
+%   Author : Sthavishtha Bhopalam Rajakumar
+%   Updated date : 30-09-2017
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clc;clear
+
+%defined parameters in the question
+d = 5e-8; %Diffusion coefficient
+u_max = 0.1; %macroscopic fluid velocity
+n = 1000; %number of grids
+t_max = 4000; %maximum time/iteration to be reached
+c_s = 1./sqrt(3); %speed of sound
+ma = u_max./c_s; %mach number
+tau = d./(c_s*c_s); %relaxation time
+beta = 1./(2*tau+1); %factor considering the effect of relaxation time
+
+%initializing initial populations
+rho = zeros(n,1); %particle densities
+for i=1:n
+    rho(i) = 1.0+0.5*exp(-5000*power(((i-1)/(n-1) - 0.25),2)); %initial density variation
+end
+feq = zeros(3,n); %equilibrium populations
+feq(1,:) = 2*rho(:).*(2 - sqrt(1 + ma*ma))./3; %defining equilibrium populations
+feq(2,:) = rho(:).*((u_max - c_s*c_s)/(2*c_s*c_s) + sqrt(1 + ma*ma))./3;
+feq(3,:) = rho(:).*((-u_max - c_s*c_s)/(2*c_s*c_s) + sqrt(1 + ma*ma))./3;
+f = feq; %initial population is in equilibrium 
+
+for t = 1:t_max     %t = time/iteration
+
+    if t~=1
+        rho(:) = f(1,:) + f(2,:) + f(3,:); %density kinetic equation
+        feq(1,:) = 2*rho(:).*(2 - sqrt(1 + ma*ma))./3; %defining equilibrium populations
+        feq(2,:) = rho(:).*((u_max - c_s*c_s)/(2*c_s*c_s) + sqrt(1 + ma*ma))./3;
+        feq(3,:) = rho(:).*((-u_max - c_s*c_s)/(2*c_s*c_s) + sqrt(1 + ma*ma))./3;
+        f(1,:) = f(1,:) - 2*beta*(f(1,:) - feq(1,:)); %performing collision 
+        f(2,:) = f(2,:) - 2*beta*(f(2,:) - feq(2,:));
+        f(3,:) = f(3,:) - 2*beta*(f(3,:) - feq(3,:));
+    end
+
+    temp_left = f(3,1); %periodic boundary conditions
+    temp_right = f(2,n);
+
+    %advection across direction of c
+    for i=n:-1:2
+        f(2,i) = f(2,i-1);
+    end
+
+    for i=1:n-1
+        f(3,i) = f(3,i+1);
+    end
+
+    f(2,1) = temp_right;
+    f(3,n) = temp_left;
+
+end
+%postprocessing
+x = 1:n;
+rho_act(x) = 1+0.5*exp(-5000*power(((x-1)/(n-1) - 0.25),2));
+plot(x,rho_act(x),'k--',x,rho(x),'k-');
+