The main objective of this project is to write a MATLAB code to calculate the horizontal average of the stream-wise velocity 𝐮, the velocity fluctuations 𝐮′, 𝐯′ and 𝐰′ based on Reynold’s decomposition, the horizontal average of (𝐮′⋅𝐮′), (𝐯′⋅𝐯′), (𝐰′⋅𝐰′) and (−𝐮′⋅𝐰). In addition, the objective is to plot the z profile of 𝐮̅ on a standard linear profile and a log-law profile to ensure the results are correct, and to plot the z profiles of all four Reynold’s stress quantities on a linear to linear plot.

​

Project Information and Results:
The application of the 3D data set relates to turbulent flow between two plates; the top plate is moving at constant speed while the bottom plate is fixed. The turbulent plot flow profile on this application would look like a S-shape, but on this project, we only perform the analysis on the lower half plate, which ranges vertically from z = 0 to z = 65, so the profile is unknown.
The provided instantaneous 3D field data shared contains the variables x, y, z, u, v, and w. In this case x, y and z are the coordinate components while u, v and w are the velocity components. All 6 values represented on each row in the original data set correspond to one point in the flow field. In this project, we perform the standard practice when having a 3D flow in cartesian coordinates. All points show fluctuation; in the x-direction, y- direction and z-direction the velocity is not uniform. This is the principal behavior of turbulent flow; nonuniformity and chaos.

The 3D flow field is divided between horizontal and vertical gridlines. Each horizontal and vertical gridline intercept with each other and generate an intersection, called grid point. All the data points are represented on these grid points. In real life simulation, a continuous function cannot be treated to perform analysis using a computer software since continuous functions have infinite points. Therefore, we need to discretize the domain by only having values defined on these grid points.
To obtain the average horizontal velocity of the flow we need to average every point along the k horizontal plane. For each velocity u, v and w there is a 3-dimensional array with 128 index points as the first dimension, 128 points as the second dimension and 65 points as the third dimension; i = 128 (streamwise direction), j = 128 (spanwise direction) and z = 65 (vertical direction or wall-normal direction). The average on each horizontal plane needs to be averaged, which means that we start with a 3-dimensional array for each fixed k and convert it to a 1-dimensional array that only changes with k. Since we average out all the variation with respect to x and y, in terms of array we average out the variation with respect to i and j, hence we obtain the averaged 1-dimensional array that still conserves the variation in terms of z. After achieving this 1-dimensional array, the z profile is plotted using MATLAB. The vertical profile of the velocity is the vertical coordinate z, and z is discretized based of the index k.

Velocity u is a 3D array in MATLAB, but mathematically it is a function of 3 coordinates, u is the x direction component of velocity, v is the y component velocity and w is the z component velocity. Velocity is also a function of time, but for project purposes we do not take time into consideration. In the original data u, v and w are functions of x, y, and z. This is analytically correct, but in standard research and analysis it is not possible to do the continuous
function using computer software. Therefore, we use a discretized form for the computer to solve a set of linear equations. The mathematical definition of 𝐮̅ is:
𝐮̅= 1𝐴𝑥𝑦⋅∬𝑢(𝑥,𝑦,𝑧)𝐴𝑥𝑦 𝑑𝑥 𝑑𝑦
This mathematical model gives us the average value of each horizontal cross-sectional plane, but it cannot be used in MATLAB. In MATLAB we use the following definition of 𝐮̅:
𝐮̅= 1nx⋅ny⋅(ΣΣu(i,j,k)𝑛𝑦𝑗=1nxi=1)
The summation is for i and j, so the k dependent will be kept. After we average 𝐮̅, we obtain it as a 1-dimensional array with index k only. Normalization of 𝐮̅ and z is needed for plotting purposes, since we want to follow the same format as presented on the textbook: u̅u∗ vs z⋅u∗v
After we normalize 𝐮̅ and z to have non-dimensional axis values, the obtained plots for the linear to linear and log-law profiles can be observed in the figure 1.

Figure 1: Average Horizontal Velocity u

The figure exhibited above represents the log-law profile (left hand side plot) and the standard linear profile (right hand side plot) of the horizontal average velocity u when plotted against z (both axis values are normalized). As is can be observed from the graph, both profiles display a very similar ascendant behavior; as the z value is increased, the 𝐮̅ value is also
increased, which technically means that the further away the flow is from the bottom plate, the higher the magnitude of the velocity vector until the vertical axis point reaches the mid-point between top and bottom plates. Since this project suppressed the top plate and only shows the bottom half plate, the maximum point should be at the highest z value. This makes sense since the plate would generate some friction on the flow slowing it down the horizontal axis. As it can be seen on the figure above, the increase in horizontal velocity is not constant, since turbulent flow shows non-constant velocity behavior.
To obtain the Reynold stress profiles of the horizontal averages (𝐮′⋅𝐮′), (𝐯′⋅𝐯′), (𝐰′⋅𝐰′) and (−𝐮′⋅𝐰), the first step was obtaining the velocity fluctuations 𝐮′, 𝐯′ and 𝐳′ based on Reynold’s decomposition. These values were calculated using the following formulas and solving for their respective values: 𝐮=u̅+u′ 𝐯=v̅+v′ 𝐰=w̅+w′
In this case 𝐮 being the velocity given in the data, 𝐮̅ being the velocity average and 𝐮′ being the velocity fluctuation. Then, the variables shown below were specified as following: 𝐑𝐱𝐱=u′⋅u′
𝐑𝐲𝐲=v′⋅v′
𝐑𝐳𝐳=w′⋅w′
𝐑𝐱𝐲=−u′⋅w′
These values represent the local value of 𝐑𝐱𝐱,𝐑𝐲𝐲,𝐑𝐳𝐳 and 𝐑𝐱𝐲 for each point. Therefore the averages of (𝐮′⋅𝐮′), (𝐯′⋅𝐯′), (𝐰′⋅𝐰′) and (−𝐮′⋅𝐰) are essentially equal to the average of 𝐑𝐱𝐱, 𝐑𝐲𝐲,𝐑𝐳𝐳 and 𝐑𝐱𝐲. Before the data was ready to be used in MATLAB, it needed to be molded; 𝐮′ needed to be another 3-dimensional array and have the same number of dimensions as u (given), and then the average of 𝐮′,𝐯′ and 𝐰′ values could be taken in the horizontal plane to then compare them to z. To obtain the Reynold’s stress in MATLAB a 3D loop for i, j and k was written. For 𝐮′,𝐯′ and 𝐰′ all I did was subtract 𝐮̅(k) from each (x,y) value in u (given). Inside the loop I defined:

𝐮′(𝐢,𝐣,𝐤)= u(i,j,k)−u̅(k)

𝐯′(𝐢,𝐣,𝐤)= v(i,j,k)−v̅(k)

𝐰′(𝐢,𝐣,𝐤)= w(i,j,k)−w̅(k)
𝐑𝐱𝐱(𝐢,𝐣,𝐤)=u′(i,j,k)⋅u′(i,j,k)

𝐑𝐲𝐲(𝐢,𝐣,𝐤)=v′(i,j,k)⋅v′(i,j,k)
𝐑𝐳𝐳(𝐢,𝐣,𝐤)=w′(i,j,k)⋅w′(i,j,k)
𝐑𝐱𝐳(𝐢,𝐣,𝐤)=−u′(i,j,k)⋅w′(i,j,k)

The normalization of 𝐑𝐱𝐱,𝐑𝐲𝐲,𝐑𝐳𝐳 and 𝐑𝐱𝐲 was required for plotting purposes. In this case, friction velocity 𝒖∗ on the denominator was squared to match the velocity squared units of the numerator. The following formulas were used:

𝑢′⋅𝑢′̅𝑢∗2
𝑣′⋅𝑣′̅𝑢∗2
𝑤′⋅𝑤′̅𝑢∗2
−𝑢′⋅𝑤′̅𝑢∗2
After we normalize 𝐑𝐱𝐱,𝐑𝐲𝐲,𝐑𝐳𝐳 and 𝐑𝐱𝐲 to have non-dimensional axis values, the obtained plots for the linear to linear profile can be observed in Figure 2. 

Figure 2: Reynold's Stress Profiles Rxx, Ryy, Rzz and Rxz

As it can be observed in the figure above, (𝐮′⋅𝐮′) shows the highest Reynold Stress values and ( 𝐯′⋅𝐯′), (𝐰′⋅𝐰′) and (−𝐮′⋅𝐰) show a very similar behavior with lower stress values. This result makes sense since 𝐮′ is the fluctuation of the flow velocity in the stream-wise direction and it should exhibit the highest Reynold’s stress value.

 

The peak value of Rxx occurs at approximately z⋅u∗v = 14.43 and 𝑢′⋅𝑢′̅𝑢∗2=8.175.

In addition, when z⋅u∗v = 0 (right on the bottom flat plate’s surface) there is no flow velocity
or Reynold’s stress. Therefore, I can assume the Reynold’s stress component is related to the turbulent motion of the flow.