MATLAB|Splart-Allmaras湍流模型及MATLAB编程~ LBM|数值模拟

还是老规矩先宣传一下QQ群群：格子玻尔兹曼救星：293267908。就是为了早点毕业建的群。
代码完整版在群里lbm_matlab-master。文章请搜索：
(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.
今天研究了一下Splart-Allmaras，因为论文要研究一下大涡。首先映入眼帘的是一句话：
文章图片

......
文章图片

好吧，康康代码：

% For the Spalart-Allmaras model in lid-driven cavity. % input variables are 2d matrices. % d: wall distances. % nut is the turbulent viscosity that gets added onto the nominal %viscosity to produce the total effective viscosity. % nu: nonminal (scalar) viscosity% apply bc for lid-driven cavity. nutilde(1,:) = 0; nutilde(end,:) = 0; nutilde(:,1) = 0; nutilde(:,end) = 0; % Model constants. sigma = 2/3; cb1 = 0.1355; cb2 = 0.622; kappa = 0.4187; cw2 = 0.3; cw3 = 2; cv1 = 7.1; % Derived model constants. cw1 = cb1 / kappa^2 + ( 1 + cb2 ) / sigma; % Derived input variables. ux = upwind_x(u,u,dh); uy = upwind_y(u,v,dh); vx = upwind_x(v,u,dh); vy = upwind_y(v,v,dh); % 这个公式是四种情况相加，i=x,j=y，i=j=y，i=x,j=y，i=y,j=x S = 2*ux.^2 + 2*vy.^2 + (uy + vx).^2 - 2/3*(ux + vy).^2; X = nutilde / nu; fv1 = X.^3 ./ ( X.^3 + cv1^3 ); Stilde = S.^0.5 .* ( 1 ./ X + fv1 ); r = tanh( nutilde ./ ( kappa^2*Stilde.*d.^2 ) ) / tanh(1.0); g = r + cw2 * ( r.^6 - 6 ); fw = g.*( ( 1 + cw3^6 ) ./ ( g.^6 + cw3^6 ) ).^(1/6); ntx = upwind_x(nutilde,u,dh); % ？ nty = upwind_y(nutilde,v,dh); ntxx = spatial_difference_x(nutilde,dh); ntyy = spatial_difference_y(nutilde,dh); dnutilde = -( u.*ntx + v.*nty ) ... + (nu + nutilde) / sigma .* ( ntxx + ntyy ) ... + (cb2 + 1) / sigma * ( ntx.^2 + nty.^2 ) ... + cb1 * Stilde .* max(nutilde,nu/100) ... - cw1 * fw .* ( nutilde ./ d ).^2; nutilde = nutilde + dt*dnutilde; % apply bc for lid-driven cavity. nutilde(1,:) = 0; nutilde(end,:) = 0; nutilde(:,1) = 0; nutilde(:,end) = 0; nut = nutilde.*fv1; %动力涡粘性系数 omega = turbulent_relaxation(nu, nut);

1、

文章图片
是四种情况相加，i=x,j=y，i=j=y，i=x,j=y，i=y,j=x，其中的k是指的xy相等的情况，所以S = 2*ux.^2 + 2*vy.^2 + (uy + vx).^2 - 2/3*(ux + vy).^2;
【MATLAB|Splart-Allmaras湍流模型及MATLAB编程~】2、

文章图片
中的(tanh)函数表达式：

文章图片
。
其MATLAB实现为：

figure('NumberTitle', 'off', 'Name', 'Tanh函数'); x=-5:0.1:5; y=2./(1+exp(-2*x))-1; plot(x,y); xlabel('X轴'); ylabel('Y轴'); %坐标轴表示对象标签 grid on; %显示网格线 axis on; %显示坐标轴 axis([-5,5,-1,1]); %x,y的范围限制 title('Tanh函数'); ———————————————— 版权声明：本文为CSDN博主「拦路雨g」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。原文链接：https://blog.csdn.net/lanluyug/article/details/76791271

文章图片

3、迎风格式 upwind difference

function d = upwind_x(s, u, dh) % takes upwind difference % indexing: s(j,i)d = 0*s; % boundaries. d(:,1) = ( s(:,2) - s(:,1) ) / dh; d(:,end) = ( s(:,end) - s(:,end-1) ) / dh; % interior. mask = u(:,2:end-1) >= 0; d(:,2:end-1) = ... (1-mask) .* s(:, (2:end-1)+1 ) ... + (-1+2*mask) .* s(:, 2:end-1 ) ... - mask .* s(:, (2:end-1)-1 );

对上下边界进行一阶差分；对方腔中部采用二阶差分。当然也有其他偏微分方程数值解法。
4、

文章图片

dnutilde = -( u.*ntx + v.*nty ) ... + (nu + nutilde) / sigma .* ( ntxx + ntyy ) ... + (cb2 + 1) / sigma * ( ntx.^2 + nty.^2 ) ... + cb1 * Stilde .* max(nutilde,nu/100) ... - cw1 * fw .* ( nutilde ./ d ).^2;

nutilde在边界处都设为0。
代码中max(nutilde,nu/100)指的是非边界区域最少要加入一个比较小的粘度，否则最终会使得nutilde到处都==0.
5、最终nut = nutilde.*fv1; omege就得到了：