BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Project

Statistical Analysis of Turbulent Flow Simulation

VVenkatesh Sundaresan3 Dec 2019
MATLAB

Turbulent flow has practical considerations in many engineering applications including urban development and green energy. The turbulence associated with air pollution dispersal has been studied to improve urban development planning. Studying turbulence in the wake of wind turbines can help improve the design and overall efficiency of a wind farm. Understanding how turbulence develops is essential to improvements in human health and technology.

A classical experiment investigating turbulent flow uses staggered cubical obstacles along the lower wall to induce a rough flow. This report analyzes a direct numerical simulation of this experiment. The DNS was completed using a Reynolds number of 2800, with the fully developed flow. The cube is positioned at a height h of 0.125 units from the bottom wall. The bounds of the domain in the vertical direction are (-1.125 to 1.0 units) from the lower to the upper wall. The spanwise direction ranges from 0 to 2 units, and the streamwise direction ranges from 0 to 8 units.

The vorticity equation can be found by taking the curl of the Navier-Stokes Equation. Start with modified pressure (p) notation,

`p=P+rhoPsi`

Where P is the pressure, the gravitational potential `Psi` = gz, where g is the gravitational acceleration and z is the vertical coordinate. The Navier-Stokes equation as seen in Turbulent Flows simplifies to the following,

`(DU)/(Dt)` = `-∇prho + v∇2(Du)/(Dt)` = `-∇(p)/rho + nu∇^2U`

Furthermore, the rate of change of the velocity can be defined by the partial derivative of the Lagrangian field and the substantial derivative of the Eulerian field

`DU/Dt ≡ ∂U/∂t + U∙ ∇ U = - ∇ p/ρ + v∇^2 U`

Noting there is non-zero vorticity, the vorticity equation can now be found by taking the curl of the Navier-Stokes Equation

`∇ × (∂U/∂t + U∙ ∇ U = -∇ p/ρ + v∇^2 U )`

`(∇ ×∂U/∂t) + (∇ ×( U∙ ∇ U) )= (-∇ ×∇ p/ρ) + (∇ ×v∇^2 U)`

Using vector identities `∇*∇`· scalar = 0 and the assumption the flow is incompressible, the pressure term is eliminated. The equation is reduced to

`(∇ ×∂U/∂t) + (∇ ×( U∙ ∇ U) )= v∇^2 (∇ ×U)`

Turbulent Flows defines vorticity as `omega = ∇*U`, using the definition of vorticity the following substitutions can be made,

`∂ω/∂t+ (∇ ×( U∙ ∇ U) )= v∇^2 ω`

Lastly using a vector identity for the second term,

`∂ω/∂t+ (∇×∇(u^2/2)-∇×(U×ω))= v∇^2 ω`

Further simplification reveals the vorticity equation,

`Dω/Dt = (ω ∙ ∇)U + v∇^2 ω`

Velocity and Vorticity Components:

In Fig. 1, velocity fluctuations are depicted at three different wall-normal levels, i.e. at the top of cubes edge (y/H = -1), next to the centerline (y/H = 0.2) and close to the top wall (y/H = 0.9). At the lower wall (y/H = 1, Figure 1a), the fluctuation of the streamwise velocity is clearly dominated by the wake behind each cube showing rapid alternation of positive and negative oscillations, whereas the wall-normal fluctuations (Figure 1d) are close to zero due to the boundary condition. Spanwise fluctuations (Figure 1g) are still present since periodic boundary conditions are imposed along z, but they are relatively smaller than the ones in the streamwise direction because no force is acting in this direction. At y/H = 0.2 (Figures 1b, 1e, 1h), close to the centerline of the tube, the fluctuation shows the same pattern in all the directions as the turbulence is almost isotropic. Weak and scattered oscillations characterize this area since no boundary is sensed by the flow. Finally, at y/H = 0.9 (Figure 1c), the presence of a smooth wall induces the streamwise velocity component to show typical high and low momentum regions elongated in the x-direction. In fact, by mass balance, the vertical mixing of momentum induced by the streamwise vorticity, associated with a strong vertical gradient, causes accelerations and deceleration in the ow. On the other hand, still because of the wall, wall-normal fluctuations (Figure 1f) are minimal due to the vicinity of the solid wall. Oscillations in the z-direction (Figure 1i) do not show any preferred path since no mean flow is driving the turbulence in this direction.

Figure 1: Velocity fluctuations in horizontal sections at different y locations: y/H = 1 (figures a), d),g)), y/H = 0.2 (figures b), e), h)), and y/H = 0.9 (figures c), f), i)). Figures a), b), c) correspond to streamwise fluctuations, figures d), e), f) to wall-normal fluctuations and figures g), h), i) to spanwise fluctuations. Figures a), d), g) report the locations of the cubes as black patches. The colormap is bounded at `pm`0.3 for all the cases.
Figure 2: Vorticity components in horizontal sections at different y locations: y/H = 1 (figures a), d), g)), y/H = 0.2 (figures b), e), h)), and y/H = 0.9 (figures c), f), i)). Figures a), b), c) correspond to streamwise vorticity `omega_x`, figures d), e), f) to wall-normal vorticity `omega_y`, and figures g), h), i) to spanwise vorticity `omega_z`. The colormap is bounded in `pm`5 for all the cases.

In Figure 2, the vorticity components are plotted at different wall-normal locations. At the lower wall (y/H = -1, Figures 2a, 2d, 2g), we can directly observe the influence of the cubes on the turbulence generation. Especially, even though the actual pattern is an instantaneous realization, a strong alternate vortex shedding is recognized after each cube in the plot of `omega_y` (Figure 2d). The averaged vorticity field would probably make this pattern more clear. A large amount of negative `omega_z` (Figure 2g) is expected since the latter is defined as: `omega_z` = `-(delu)/(dely)`, hence the strong shear next to the upper side of the cubes motivates the visualized pattern. Consistently with the plot of velocity fluctuations at y/H = 0.2, in case of isotropic turbulence, no preferred direction is detected in the vorticity field, thus the three distributions (Figures 2b, 2e, 2h) are similar. The small velocity gradient next to the centerline (see Figure 1b, 1e, 1h) motivates the low magnitude of the vorticity here. Close to the upper wall (y/H = 0.9), it is noteworthy that the distribution of `omega_x` is characterized by streaks with opposite sign moving in the spanwise direction (Figure 2c). In fact, here `omega_x` is dominated by the term `(delv)/(delz)`: the secondary motion induced by the streamwise vorticity causes a spanwise alternation of positive and negative fluctuations of v, thus leading to the actual pattern of `omega_x`. At the same y-location, `omega_y` reveals couples of streamwise elongated vortices of opposite sign (Figure 2c). This is motivated by the strong impact of `(delu)/(delz)`in the evaluation of `omega_y`: the long, coherent streaks of u0 (Figure 1c) with alternating sign imply the generation of a strong spanwise gradient leading to the actual pattern of spanwise vorticity. Finally, the large positive `omega_z` close to the top boundary (Figure 2i) is originated by the presence of a strong wall-normal gradient in streamwise velocity (as its alter ego close to the bottom wall).

Figure 3: Left: Visualizations of horizontal planes of vortex stretching at y/H = -0.9 (a), y/H = 0.2(b) and y/H = 1 (c). Right: Visualizations of horizontal planes of `(delu)/(delx)` at y/H = -0:9 (d), y/H = 0:2 (e) and y/H = 1.

Vortex stretching is depicted in Figures 3a, 3b, 1c at three different y-locations (reported in external units), while the x-component of the streamwise velocity gradient is reported for completeness in Figures 3d, 3e, 3f. Each of the wall-normal locations depicts a radical different x-z pattern. At y/H = -1, on top of the cubes, the presence of a rough surface underneath the flow causes a highly irregular path of `(delu)/(delx)`(Figure 1d). This phenomenon is reflected in the vortex stretching, whose pattern is coherent for intervals smaller than the cube spacing (Figure 3a). On the other hand, when the assumed y-location is close to the centerline of the channel (i.e. y/H = 0.2 as in Figures 1b and 1e) the vortex stretching is greatly reduced with respect to the values at the boundaries since the vorticity is less strong here than at the solid edges of the domain (compare Figures 2b, e, h with Figures 2a, d, g). The pattern of `(delu)/(delx)` close to the centerline is locally different from 0 since we are visualizing an instantaneous flow rather than averaged quantities. However, the time-averaged value of the streamwise gradient is expected to be zero since the flow is fully developed. Finally, when y/H = 0.9 (Figures 3c and 3f), the presence of a flat wall induces streamwise coherent paths in the vortex stretching related to the local behavior of the vorticity, which is organized in streaks typically populating the turbulent boundary layer. In any case, it is noted that the vorticity exhibits longer stripes with respect to the vortex stretching in the streamwise direction. This is motivated by the rapid change of sign of `(delU)/(delx)`related to the local action of the turbulent mixing. In fact, the ejection and sweep events induced by the streamwise vortices are associated respectively to deceleration and acceleration of the flow by mass conservation, thus inducing the changes in sign of Figure 1f.

Root-mean-square wall-normal profiles:

Figure 4: a) Wall-normal profiles of the r.m.s. of three velocity fluctuations. b) Wall normal profiles of the r.m.s. fluctuation of three vorticity components.

Given a quantity q(x, y, z) having zero average, the wall-normal pro le of root-mean-square (r.m.s.) may be evaluated as:

where Nx and Nz are respectively the numbers of grid points along x and z. Thus, for each wall-normal location, the average value over the x-z plane has been removed from each velocity and vorticity components and the r.m.s. have been calculated by means of the above equation. The results are depicted in Figure 4. The r.m.s. of u0 in Figure 4a is null at both boundaries, then two maxima are rapidly reached. The left peak is about y/H = -1, i.e. around the crests of the cubes, confirming the augment of turbulence intensity due to the wall roughness already inferred from Figures 1a, 1d, 1g. A comparison between the two peak values reveals that a rough wall accounts for a more intense turbulent motion with respect to a smooth boundary. In any case, attention must be paid since an overestimation is done when the r.m.s. is calculated without removing the dispersive component of velocity. As -1:125 `le` y/H `le` -1, the local mean velocity is less than the overall averaged value because of the recirculation zones attached to the cubes. Thus, this modulation of the mean velocity leads to an overestimation of the fluctuating part if it is not properly removed (e.g. by means of a triple decomposition, see [Hussain & Reynolds (1970)]). However, as shown in [Squire et al. (2016)], the progressive increase of roughness leads to a transition from friction drag (in case of the smooth wall) to form drag (related to the pressure field created around the friction elements) which destroys the near-wall cycle. Thus, the r.m.s. associated with the sole turbulent fluctuation is expected to be comparable to the smooth wall case or less than it.

Still observing Figure 4a, the r.m.s. peak in wall-normal fluctuations is less than `u_(rms)` on both sides since the solid wall prevents the development of turbulence along y. However, the obtained profiles of velocity r.m.s. are in good agreement with [Kim, Moin & Moser (1987)] in the proximity of the smooth wall. On the other hand, their r.m.s. profiles collapse very well at y/H = 0 whereas for the present case the clustering is poor and the region of minima is located about y/H = 0.2. The former issue may be due to the fact that an instantaneous eld is treated here rather than an average one, so the so-obtained second-order statistics do not represent converged values. The shift in the minima location is expected since an augment of turbulent motion is caused by the rough wall, thus making the pro les non-symmetric and shifting the r.m.s. minima away from the rough wall.

In Figure 4b, at the upper wall (y/H = 1), < `omega_z`,rms > has a nonzero value at the wall because of the presence of `(delu)/(dely)`, while no fluctuation in spanwise and streamwise vorticity is possible due to the presence of the solid boundary. As the vortex lines are tilted and lifted away from the wall, the r.m.s. values of all three components become nearly equal ([Balint, Wallace & Vukoslavcevic (1991)]) leading to a quasi-isotropic turbulent state. For the bottom wall (y/H = -1.25), two peaks of equal magnitude for `omega_x`,rms and `omega_z`,rms are found. This is explained considering that, locally, the cubes impose vortex breaks and irregularity both in x and in z with the same probability, thus the statistical behavior of `omega_x` and `omega_z` is the same. Conversely, at y/H = -1, we have: `omega_z`,rms > `omega_x`,rms since here < `omega_x` >`approx` 0 while < `omega_z` > is negative because of the positive `(delu)/(dely)` term. Thus, the modulating effect of the shear implies an overestimation of the fluctuations in `omega_z`. In the same location, it is noticeable that < `omega_y`,rms > reaches a local maximum. Recalling the vortex-shedding pattern of Figure 2d, wall-normal vortices released in the wake of each cube motivate the high value of the r.m.s. here.

P.D.F. analysis of the velocity components

The probability density function (p.d.f.) for the three velocity components has been evaluated at three different wall-normal locations. 20 bins are chosen for each component. Observing Figure 5, we find out that all the velocity components show a nearly symmetric distribution about a mean value. Thus, we model the p.d.f. of U using the following Gaussian distribution:

Mean (`mu`) and standard deviation (`sigma`) of the Gaussian distributions are summarized in Table 1.

Table 1: Summary of the mean and standard deviation values for the p.d.f. of the velocity components at different y-locations.

For the p.d.f. analysis of the streamwise velocity u (Figure 5a-c), the mean velocity increases from the lower layer (y/H = -0.94) to the upper layer (y/H = 0.2), which makes sense because of the wall boundary condition with mounted cubes exerts a slowing effect on the flow. This is reflected in the values of `mu-`columns in Table 1. More, the cubes lead to turbulence production. Thus, `sigma` values are much higher near the lower wall (y/H = -0.94 and y/H = -0.8) than close to the centerline (y/H = 0.2). When the flow reaches the centerline (y/H = 0.2), the streamwise velocity tends to be higher without any barrier and then fluctuations are reduced, as confirmed by the `sigma-`value of u at y/H = 0.2. However, referring to the p.d.f. analysis of the wall-normal velocity v (figure 5d-f) and the spanwise velocity w (figure 5g-i), we find the mean velocity is almost zero in the two directions regardless of the different y-locations, which means that the streamwise velocity u dominates this ow. This is also con firmed by the mean values of v and w in Table 1. Hence, consistency with the r.m.s. analysis of x2.3 is found for all the components along the wall-normal direction.

Figure 5: Probability density function of the three velocity components at different y locations: y/H = -0.94 (Figures a, d, g), y/H = -0.8 (Figures b, e, h), and y/H = 0:2 (Figures c, f, i). Streamwise velocity p.d.f. is reported in Figures a, b, c. Wall-normal velocity p.d.f. is depicted if Figures d, e, f and p.d.f. of the spanwise component are given in Figures g, h, i. Vertical blue bars represent the number of occurrences for each interval, while the Gaussian is reported with a red line.

MATLAB Codes

%% The following code is broken into sections
% % The first section calculates the required parameters and the rest of
% % the sections are for different plots.
%% Variable definitions:
% % yp1 = coordinates in the spanwise direction
% % yp2 = coordinates in the wall normal direction
% % yp3 = coordinates in the streamwise direction
% % q1 = spanwise velocity (w)
% % q2 = wall normal velocity (v)
% % q3 = streamwise velocity (u)
% % w1 = spanwise vorticity
% % w2 = wall normal vorticity
% % w3 = streamwise vorticity
% % w_stretch = vorticity stretching in streamwise direction
% % index = indices for y = -0.94, -0.8 and 0.2
% % cdf_qij, pdf_qij, v_range_qij = CDF, PDF and range of values for i-th 
% %                                 velocity at j-th y-value.

%% Calculation of different parameters
clear all; close all;
clc;

fn='field.data_0210.nc';
xgrid=load('xgrid.plo');
ygrid=load('ygrid.plo');
zgrid=load('zgrid.plo');

yp1=xgrid(:,2);
yp2=ygrid(:,2);
yp3=zgrid(:,2);

n1m=length(yp1); % spanwise
n2m=length(yp2); % wall normal
n3m=length(yp3); % streamwise

q3=ncread(fn,'q3');
q2=ncread(fn,'q2');
q1=ncread(fn,'q1');
p = ncread(fn,'pr');

y_bottom = 1; y_top = n2m; y_mid = 62; y_stag = 26;

[w1,w2,w3] = calc_vorticity(q1,q2,q3,yp1,yp2,yp3);
w_mag = sqrt(w1.^2 + w2.^2 +w3.^2);
w_stretch = zeros(n1m,n2m,n3m);
for i=1:n1m
    if i==1 || i==n1m
        i1 = n1m-1; i2 = 2;
    else
        i1 = i-1; i2 = i+1;
    end
    for j=1:n2m
        if j==1 || j==n2m
            j1 = n2m-1; j2 = 2;
        else
            j1 = j-1; j2 = j+1;
        end
        for k=1:n3m
            if k==1 || k==n3m
                k1 = n3m-1; k2 = 2;
            else
                k1 = k-1; k2 = k+1;
            end
            dqdy3 = ((yp2(j2)-yp2(j))^2*(q3(i,j,k)-q3(i,j1,k)) + (yp2(j)-yp2(j1))^2*(q3(i,j2,k)-q3(i,j,k)))/...
                    ((yp2(j2)-yp2(j))*(yp2(j)-yp2(j1))*(yp2(j2)-yp2(j1)));
            w_stretch(i,j,k) = w3(i,j,k)*(q3(i,j,k2)-q3(i,j,k1))/(yp3(k2)-yp3(k1))...
                             + w2(i,j,k)*dqdy3...
                             + w1(i,j,k)*(q3(i2,j,k)-q3(i1,j,k))/(yp1(i2)-yp1(i1));
        end
    end
end

q1_mean = zeros(n2m,1);
q2_mean = zeros(n2m,1);
q3_mean = zeros(n2m,1);
w1_mean = zeros(n2m,1);
w2_mean = zeros(n2m,1);
w3_mean = zeros(n2m,1);

for i=1:n2m
    q1_mean(i) = mean(q1(:,i,:),'all');
    q2_mean(i) = mean(q2(:,i,:),'all');
    q3_mean(i) = mean(q3(:,i,:),'all');
    w1_mean(i) = mean(w1(:,i,:),'all');
    w2_mean(i) = mean(w2(:,i,:),'all');
    w3_mean(i) = mean(w3(:,i,:),'all');
end

index = [34 42 65];
n_tot = n1m*n3m;
[cdf_q31,pdf_q31,v_range_q31] = calc_distribution(q3,index(1),n_tot);
[cdf_q21,pdf_q21,v_range_q21] = calc_distribution(q2,index(1),n_tot);
[cdf_q11,pdf_q11,v_range_q11] = calc_distribution(q1,index(1),n_tot);
[cdf_q32,pdf_q32,v_range_q32] = calc_distribution(q3,index(2),n_tot);
[cdf_q22,pdf_q22,v_range_q22] = calc_distribution(q2,index(2),n_tot);
[cdf_q12,pdf_q12,v_range_q12] = calc_distribution(q1,index(2),n_tot);
[cdf_q33,pdf_q33,v_range_q33] = calc_distribution(q3,index(3),n_tot);
[cdf_q23,pdf_q23,v_range_q23] = calc_distribution(q2,index(3),n_tot);
[cdf_q13,pdf_q13,v_range_q13] = calc_distribution(q1,index(3),n_tot);

[q3_rms] = calc_rms(q3,n1m,n2m,n3m);
[q2_rms] = calc_rms(q2,n1m,n2m,n3m);
[q1_rms] = calc_rms(q1,n1m,n2m,n3m);
[w1_rms] = calc_rms(w1,n1m,n2m,n3m);
[w2_rms] = calc_rms(w2,n1m,n2m,n3m);
[w3_rms] = calc_rms(w3,n1m,n2m,n3m);
[w_stretch_rms] = calc_rms(w_stretch,n1m,n2m,n3m);

%% Plotting vorticity and streamwise velcity fluctuations

figure(1)
subplot(221)
contourf(yp3,yp1,squeeze(w3(:,y_bottom,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Vorticity at y = -1.125','FontSize',22)
colormap(jet(15)); colorbar;
subplot(222)
contourf(yp3,yp1,squeeze(w3(:,y_stag,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Vorticity at y = -1','FontSize',22)
colormap(jet(15)); colorbar;
subplot(223)
contourf(yp3,yp1,squeeze(w3(:,y_mid,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Vorticity at y = 0','FontSize',22)
colormap(jet(15)); colorbar;
subplot(224)
contourf(yp3,yp1,squeeze(w3(:,y_top,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Vorticity at y = 1','FontSize',22)
colormap(jet(15)); colorbar;

figure(2)
subplot(221)
contourf(yp3,yp1,squeeze(w2(:,y_bottom,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Wall Normal Vorticity at y = -1.125','FontSize',22)
colormap(jet(15)); colorbar;
subplot(222)
contourf(yp3,yp1,squeeze(w2(:,y_stag,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Wall Normal Vorticity at y = -1','FontSize',22)
colormap(jet(15)); colorbar;
subplot(223)
contourf(yp3,yp1,squeeze(w2(:,y_mid,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Wall Normal Vorticity at y = 0','FontSize',22)
colormap(jet(15)); colorbar;
subplot(224)
contourf(yp3,yp1,squeeze(w2(:,y_top,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Wall Normal Vorticity at y = 1','FontSize',22)
colormap(jet(15)); colorbar;

figure(3)
subplot(221)
contourf(yp3,yp1,squeeze(w1(:,y_bottom,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Spanwise Vorticity at y = -1.125','FontSize',22)
colormap(jet(15)); colorbar;
subplot(222)
contourf(yp3,yp1,squeeze(w1(:,y_stag,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Spanwise Vorticity at y = -1','FontSize',22)
colormap(jet(15)); colorbar;
subplot(223)
contourf(yp3,yp1,squeeze(w1(:,y_mid,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Spanwise Vorticity at y = 0','FontSize',22)
colormap(jet(15)); colorbar;
subplot(224)
contourf(yp3,yp1,squeeze(w1(:,y_top,:)),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Spanwise Vorticity at y = 1','FontSize',22)
colormap(jet(15)); colorbar;

figure(4)
subplot(221)
contourf(yp3,yp1,squeeze(q3(:,y_bottom,:))-q3_mean(y_bottom),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Velocity Fluctuation at y = -1.125','FontSize',22)
colormap(jet(15)); colorbar;
subplot(222)
contourf(yp3,yp1,squeeze(q3(:,y_stag,:))-q3_mean(y_stag),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Velocity Fluctuation at y = -1','FontSize',22)
colormap(jet(15)); colorbar;
subplot(223)
contourf(yp3,yp1,squeeze(q3(:,y_mid,:))-q3_mean(y_mid),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Velocity Fluctuation at y = 0','FontSize',22)
colormap(jet(15)); colorbar;
subplot(224)
contourf(yp3,yp1,squeeze(q3(:,y_top,:))-q3_mean(y_top),50,'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Velocity Fluctuation at y = 1','FontSize',22)
colormap(jet(15)); colorbar;

%% Plotting vortex stretching for streamwise vorticity

figure(5)
subplot(221)
contourf(yp3,yp1,squeeze(w_stretch(:,y_bottom,:)),'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Vorticity Stretching at y = -1.125','FontSize',22)
colormap(parula(10)); colorbar;
subplot(222)
contourf(yp3,yp1,squeeze(w_stretch(:,y_stag,:)),'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Vorticity Stretching at y = -1','FontSize',22)
colormap(parula(10)); colorbar;
subplot(223)
contourf(yp3,yp1,squeeze(w_stretch(:,y_mid,:)),'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Vorticity Stretching at y = 0','FontSize',22)
colormap(parula(10)); colorbar;
subplot(224)
contourf(yp3,yp1,squeeze(w_stretch(:,y_top,:)),'EdgeColor','none');
xlabel('X','FontSize',16); ylabel('Z','FontSize',16);
title('Streamwise Vorticity Stretching at y = 1','FontSize',22)
colormap(parula(10)); colorbar;

%% Plotting RMS of three velocity and vorticity components

figure(6)
plot(yp2,q3_rms,'m-^','LineWidth',2,'MarkerSize',8);
hold on
plot(yp2,q2_rms,'b-o','LineWidth',2,'MarkerSize',7);
hold on
plot(yp2,q1_rms,'r-x','LineWidth',2,'MarkerSize',14);
hold on
plot([-1 -1],[0 0.9],'k--','LineWidth',1.5)
hold off
xlim([-1.125 1]); ylim([0 0.9]); grid on;
xlabel('Y','FontSize',20); ylabel('q_{rms}','FontSize',24);
lgd = legend('Streamwise','Wall Normal','Spanwise');
lgd.FontSize = 18; ax = gca; ax.FontSize = 20;
title('Root Mean Square of Velocity vs. Wall Distance','FontSize',18);

figure(7)
plot(yp2,w3_rms,'m-^','LineWidth',2,'MarkerSize',8);
hold on
plot(yp2,w2_rms,'b-o','LineWidth',2,'MarkerSize',7);
hold on
plot(yp2,w1_rms,'r-x','LineWidth',2,'MarkerSize',12);
hold on
plot([-1 -1],[0 18],'k--','LineWidth',1.5)
hold off
xlim([-1.125 1]); ylim([0 18]); grid on;
xlabel('Y','FontSize',20); ylabel('w_{rms}','FontSize',24);
lgd = legend('Streamwise','Wall Normal','Spanwise');
lgd.FontSize = 18; ax = gca; ax.FontSize = 20;
title('Root Mean Square of Vorticity vs. Wall Distance','FontSize',18);

figure(8)
plot(yp2,w_stretch_rms,'r-^','LineWidth',2,'MarkerSize',8);
hold on
plot([-1 -1],[0 18],'k--','LineWidth',1.5)
xlim([-1.125 1]); ylim([0 18]); grid on;
xlabel('Y','FontSize',20); ylabel('w\_stretch_{rms}','FontSize',24);
title('Root Mean Square of Streamwise Vortex Stretching vs. Wall Distance','FontSize',18);

%% Plotting Velocity PDFs

figure(9)
plot(v_range_q31,pdf_q31,'m-^','LineWidth',2,'MarkerSize',6);
hold on
plot(v_range_q32,pdf_q32,'b-o','LineWidth',2,'MarkerSize',5);
hold on
plot(v_range_q33,pdf_q33,'r-x','LineWidth',2,'MarkerSize',10);
hold off
xlabel('q_3','FontSize',20); ylabel('f(q_3)','FontSize',22);
xlim([-0.1 1]); grid on; ax = gca; ax.FontSize = 20;
lgd = legend('Y = -0.94','Y = -0.8','Y = 0.2'); lgd.FontSize = 18;
title('Streamwise Velocity PDF for different Wall Distances','FontSize',18);

figure(10)
plot(v_range_q21,pdf_q21,'m-^','LineWidth',2,'MarkerSize',6);
hold on
plot(v_range_q22,pdf_q22,'b-o','LineWidth',2,'MarkerSize',5);
hold on
plot(v_range_q23,pdf_q23,'r-x','LineWidth',2,'MarkerSize',10);
hold off
xlabel('q_2','FontSize',20); ylabel('f(q_2)','FontSize',22);
xlim([-0.3 0.4]); grid on; ax = gca; ax.FontSize = 20;
lgd = legend('Y = -0.94','Y = -0.8','Y = 0.2'); lgd.FontSize = 18;
title('Wall Normal Velocity PDF for different Wall Distances','FontSize',18);

figure(11)
plot(v_range_q11,pdf_q11,'m-^','LineWidth',2,'MarkerSize',6);
hold on
plot(v_range_q12,pdf_q12,'b-o','LineWidth',2,'MarkerSize',5);
hold on
plot(v_range_q13,pdf_q13,'r-x','LineWidth',2,'MarkerSize',10);
hold off
xlabel('q_1','FontSize',20); ylabel('f(q_1)','FontSize',22);
xlim([-0.4 0.4]); grid on; ax = gca; ax.FontSize = 20;
lgd = legend('Y = -0.94','Y = -0.8','Y = 0.2'); lgd.FontSize = 18;
title('Spanwise Velocity PDF for different Wall Distances','FontSize',18);

figure(12)
plot(v_range_q31,pdf_q31,'m-^','LineWidth',2,'MarkerSize',6);
hold on
plot(v_range_q21,pdf_q21,'b-o','LineWidth',2,'MarkerSize',5);
hold on
plot(v_range_q11,pdf_q11,'r-x','LineWidth',2,'MarkerSize',10);
hold off
xlabel('V','FontSize',20); ylabel('f(V)','FontSize',22);
xlim([-0.4 1]); grid on; ax = gca; ax.FontSize = 20;
lgd = legend('Streamwise','Wall Normal','Spanwise'); lgd.FontSize = 18;
title('Velocity PDF for different Components at Y = -0.94','FontSize',18);

figure(13)
plot(v_range_q32,pdf_q32,'m-^','LineWidth',2,'MarkerSize',6);
hold on
plot(v_range_q22,pdf_q22,'b-o','LineWidth',2,'MarkerSize',5);
hold on
plot(v_range_q12,pdf_q12,'r-x','LineWidth',2,'MarkerSize',10);
hold off
xlabel('V','FontSize',20); ylabel('f(V)','FontSize',22);
xlim([-0.4 1]); grid on; ax = gca; ax.FontSize = 20;
lgd = legend('Streamwise','Wall Normal','Spanwise'); lgd.FontSize = 18;
title('Velocity PDF for different Components at Y = -0.8','FontSize',18);

figure(14)
plot(v_range_q33,pdf_q33,'m-^','LineWidth',2,'MarkerSize',6);
hold on
plot(v_range_q23,pdf_q23,'b-o','LineWidth',2,'MarkerSize',5);
hold on
plot(v_range_q13,pdf_q13,'r-x','LineWidth',2,'MarkerSize',10);
hold off
xlabel('V','FontSize',20); ylabel('f(V)','FontSize',22);
xlim([-0.4 1]); grid on; ax = gca; ax.FontSize = 20;
lgd = legend('Streamwise','Wall Normal','Spanwise'); lgd.FontSize = 18;
title('Velocity PDF for different Components at Y = 0.2','FontSize',18);
function [w1,w2,w3] = calc_vorticity(q1,q2,q3,yp1,yp2,yp3)
%--------------------------------------------------------------------------
% Function to calculate vorticity based on the central difference scheme
%
% Variable definitions:
%   n1m = number of nodes in the spanwise direction
%   n2m = number of nodes in the wall normal direction
%   n3m = number of nodes in the streamwise direction
%--------------------------------------------------------------------------

    n1m = length(yp1); n2m = length(yp2); n3m = length(yp3);
    w1 = zeros(n1m,n2m,n3m);
    w2 = zeros(n1m,n2m,n3m);
    w3 = zeros(n1m,n2m,n3m);
    for i=1:n1m
        if i==1 || i==n1m
            i1 = n1m-1; i2 = 2;
        else
            i1 = i-1; i2 = i+1;
        end
        for j=1:n2m
            if j==1 || j==n2m
                j1 = n2m-1; j2 = 2;
            else
                j1 = j-1; j2 = j+1;
            end
            for k=1:n3m
                if k==1 || k==n3m
                    k1 = n3m-1; k2 = 2;
                else
                    k1 = k-1; k2 = k+1;
                end
                dqdy1 = ((yp2(j2)-yp2(j))^2*(q1(i,j,k)-q1(i,j1,k)) + (yp2(j)-yp2(j1))^2*(q1(i,j2,k)-q1(i,j,k)))/...
                    ((yp2(j2)-yp2(j))*(yp2(j)-yp2(j1))*(yp2(j2)-yp2(j1)));
                dqdy3 = ((yp2(j2)-yp2(j))^2*(q3(i,j,k)-q3(i,j1,k)) + (yp2(j)-yp2(j1))^2*(q3(i,j2,k)-q3(i,j,k)))/...
                    ((yp2(j2)-yp2(j))*(yp2(j)-yp2(j1))*(yp2(j2)-yp2(j1)));
                w1(i,j,k) = (q2(i,j,k2)-q2(i,j,k1))/(yp3(k2)-yp3(k1)) - dqdy3;
                w2(i,j,k) = (q3(i2,j,k)-q3(i1,j,k))/(yp1(i2)-yp1(i1)) - (q1(i,j,k2)-q1(i,j,k1))/(yp3(k2)-yp3(k1));
                w3(i,j,k) = dqdy1 - (q2(i2,j,k)-q2(i1,j,k))/(yp1(i2)-yp1(i1));
            end
        end
    end
end

ROOT MEAN SQUARE VALUE CALCULATION FUNCTION:

function [q_rms] = calc_rms(q,n1m,n2m,n3m)
%--------------------------------------------------------------------------
% Function for calculating rms value of a quantity 'q' in the horizontal
% planes (rms over streamwise and spanwise direction)
%--------------------------------------------------------------------------
    n_tot = n1m*n3m;
    q_rms = zeros(1,n2m);
    sum = 0;
    for i=1:n2m
        for j=1:n1m
            for k=1:n3m
                sum = sum + q(j,i,k)^2;
            end
        end
        q_rms(i) = sqrt(sum/n_tot);
        sum = 0;
    end
end

CUMULATIVE DISTRIBUTION FUNCTION CALCULATION:

function [cdf,pdf,v_range] = calc_distribution(q,index,n_tot)
%--------------------------------------------------------------------------
% Function to calculate the cumulative distribution function (CDF) and the
% probability density function (PDF).
%
% Variable definitions:
%   n_tot = total number of nodes in streamwise and spanwise direction
%--------------------------------------------------------------------------
    v_range = min(q(:,index,:),[],'all'):0.01:max(q(:,index,:),[],'all');
    cdf = zeros(1,length(v_range));
    for i=1:length(v_range)
        cdf(i) = length(find(q(:,index,:)<=v_range(i)))/n_tot;
    end
    pdf = gradient(cdf,v_range);
end

Conclusions

In this work, instantaneous realizations of velocity and vorticity components in a turbulent channel flow have been examined at different wall-normal locations. Both horizontal-plane visualizations and averaged statistics have been involved to give an insight into the turbulent motion inside the domain. The results show an increased level of turbulence close to the staggered-cubes rough wall with respect to the smooth boundary confirming the role of the wall roughness in turbulence generation. On the other side, a classic smooth-wall turbulent boundary layer has been detected, whose characteristics are in accordance with previous literature results summarized in [Kim, Moin & Moser (1987)]. With respect to that study, the present case shows a lack of symmetry in the wall-normal direction motivated by the different roughness levels on the two bounded sides along. However, even with an instantaneous snapshot of velocity, it has been possible to reconstruct the main features of a channel ow including the effect of the wall roughness on the turbulent process.

1 comment
  1. Wenkai2fluidWhere to download the .nc .plo files?
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.