SPALART–ALLMARAS TURBULENCE MODEL--Sa大涡的方腔流代码--感谢Robert Lee

还是老规矩先宣传一下QQ群群： 格子玻尔兹曼救星：293267908。

%文章来源：

AN IMPLEMENTATION OF THE SPALART–ALLMARAS TURBULENCE MODEL IN A MULTI-DOMAIN
  LATTICE BOLTZMANN METHOD FOR SOLVING TURBULENT AIRFOIL FLOWS.  (2015) Nicolas Pellerin, Sebastien Leclaire, and Marcelo Reggio.
中文解释清参考博客：https://blog.csdn.net/weixin_37783345/article/details/104609941

% A Lattice Boltzmann (single relaxation time) D2Q9 solver,
% with the Spalart Allmaras turbulence model, on a lid-driven cavity.
% Cell centers (nodes) are placed on the boundaries.
% Author: Robert Lee
% Email: rlee32@gatech.edu

% UNDER CONSTRUCTION
%E:\matlabwork\lbm_matlab-master\navier_stokes
%文章来源：(2015) Nicolas Pellerin, Sebastien Leclaire, and Marcelo Reggio. 
%  AN IMPLEMENTATION OF THE SPALART–ALLMARAS TURBULENCE MODEL IN A MULTI-DOMAIN 
% LATTICE BOLTZMANN METHOD FOR SOLVING TURBULENT AIRFOIL FLOWS.
%中文：https://www.doc88.com/p-2901681623020.html

% A Lattice Boltzmann (single relaxation time) D2Q9 solver,
% with the Spalart Allmaras turbulence model, on a lid-driven cavity. 
% Cell centers (nodes) are placed on the boundaries. 
% Author: Robert Lee
% Email: rlee32@gatech.edu

clear;close all;clc;
 
% Algorithm steps:
% Initialize meso (f)
% Apply meso BCs
% Determine macro variables and apply macro BCs
% Loop:
%   Collide
%   Apply meso BCs
%   Stream
%   Apply meso BCs?
%   Determine macro variables and apply macro BCs

% D2Q9 collisions on 2-d matrix.
w = zeros(9,1);
w(1) = 4/9;
w(2:5) = 1/9;
w(6:9) = 1/36;
c = zeros(9,2);
c(1,:) = [0, 0];
c(2,:) = [1, 0];
c(3,:) = [0, 1];
c(4,:) = [-1, 0];
c(5,:) = [0, -1];
c(6,:) = [1, 1];
c(7,:) = [-1, 1];
c(8,:) = [-1, -1];
c(9,:) = [1, -1];

% Physical parameters.
L_p = 4; %1.1; % Cavity dimension. 
U_p = 1; %1.1; % Cavity lid velocity.
nu_p = 1.2e-3; % 1.586e-5; % Physical kinematic viscosity.
rho0 = 1;
% Discrete/numerical parameters.
nodes = 100;
dt = .002;
timesteps = 10000;
nutilde0 = 1e-5; % initial nutilde value (should be non-zero for seeding).

% Derived nondimensional parameters.
Re = L_p * U_p / nu_p;
disp(['Reynolds number: ' num2str(Re)]);
% Derived physical parameters.
t_p = L_p / U_p;
disp(['Physical time scale: ' num2str(t_p) ' s']);
% Derived discrete parameters.