art2模型 matlab,Splart-Allmaras湍流模型及MATLAB编程~

还是老规矩先宣传一下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;

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就得到了：

tau = 3 * ( nu_lb + nut ) + 0.5;

omega = 1 ./ tau;