Airflow Simulation in Empty / Occupied Rooms and Environments

     One fine morning, I decided to code the Navier–Stokes equations [read if you are bored 🤣] . This post has the results of  flow simulation inside close environments of various aspect ratios. As is customary with all my CFD work using commercial and 🏡made CFD codes, this too is inspired by the free lectures of Dr. Lorena Barba.

The first case has an aspect ratio of 1:1 while for the second case, the aspect ratio is at 3:1. The airflow is at 0.4555 m/s for both cases. Both cases are isothermal at 293 K. There is no turbulence 🌬model [free code ] 🤑.

The smallest resolved scale (~4x smallest mesh size) for 1:1 case is ~ 0.0045 m and for the 3:1 case is at 0.0112 m. Time scales (~4x time-step size) are at 0.0004 s and 0.0004 s, respectively. Fig. 1 shows results for 1:1 aspect ratio while the results for 3:1 case are shown in Fig. 2. For both cases, flow enters from top-left and exists at bottom-right of the rooms. The boundary conditions are taken from [1]. I compared the results with Fluent simulations I ran at same boundary conditions and stopped my simulations when eye-balling didn't revealed any difference 😆. What you expect? This is a blog not a journal 😝.

Fig. 1, 1:1 aspect ratio

Fig. 2, 3:1 aspect ratio

 

Fig. 3, 3:1 aspect ratio with partition

Fig. 4, Flow inside ducts

     If you want to hire me as your PhD student in the research projects related to turbo-machinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

References

     [1] Horikiri, Kana & Yao, Yufeng. Validation Study of Convective Airflow in an Empty Room, "Recent Researches in Energy, Environment, Devices, Systems, Communications and Computers", ISBN: 978-960-474-284-4

Lid-Driven Cavity MATLAB Code

     MATLAB code for 2D Lid-Driven Cavity. Includes labeled commands, plotting and is less than 100 lines of code. Resume possible, to resume comment close all, clear, clc, u, v, p and then run the code. Results are available here:

%% clear and close

close all

clear

clc

%% define spatial and temporal grids

h = 1/10; % grid spacing

cfl = h; % cfl number

L = 1; % cavity length

D = 1; % cavity depth

Nx = round((L/h)+1); % grid points in x-axis

Ny = round((D/h)+1); % grid points in y-axis

nu = 0.000015111; % kinematic viscosity

Uinf = 0.0015111; % free stream velocity / lid velocity

dt = h * cfl / Uinf; % time step

travel = 2; % times the disturbance travels entire length of computational domain

TT = travel * L / Uinf; % total time

ns = TT / dt; % number of time steps

l_square = 1; % length of square

Re = l_square * Uinf / nu; % Reynolds number

rho = 1.2047; % fluid density

%% initialize flowfield

u = zeros(Nx,Ny); % x-velocity

v = zeros(Nx,Ny); % y-velocity

p = zeros(Nx,Ny); % pressure

i = 2:Nx-1; % spatial interior nodes in x-axis

j = 2:Ny-1; % spatial interior nodes in y-axis

[X, Y] = meshgrid(0:h:L, 0:h:D); % spatial grid

maxNumCompThreads('automatic'); % select CPU cores

%% solve 2D Navier-Stokes equations

for nt = 1:ns

    pn = p;

    p(i, j) = ((pn(i+1, j) + pn(i-1, j)) * h^2 + (pn(i, j+1) + pn(i, j-1)) * h^2) ./ (2 * (h^2 + h^2)) ...

        - h^2 * h^2 / (2 * (h^2 + h^2)) * (rho * (1 / dt * ((u(i+1, j) - u(i-1, j)) / (2 * h) + (v(i, j+1) - v(i, j-1)) / (2 * h)))); % pressure poisson

    p(1, :) = p(2, :); % dp/dx = 0 at x = 0

    p(Nx, :) = p(Nx-1, :); % dp/dx = 0 at x = L

    p(:, 1) = p(:, 2); % dp/dy = 0 at y = 0

    p(:, Ny) = 0; % p = 0 at y = D 

    un = u;

    vn = v;

    u(i, j) = un(i, j) - un(i, j) * dt / (2 * h) .* (un(i+1, j) - un(i-1, j)) ...

        - vn(i, j) * dt / (2 * h) .* (un(i, j+1) - un(i, j-1)) - dt / (2 * rho * h) * (p(i+1, j) - p(i-1, j)) ...

        + nu * (dt / h^2 * (un(i+1, j) - 2 * un(i, j) + un(i-1, j)) + dt / h^2 * (un(i, j+1) - 2 * un(i, j) + un(i, j-1))); % x-momentum

    u(1, :) = 0; % u = 0 at x = 0

    u(Nx, :) = 0; % u = 0 at x = L

    u(:, 1) = 0; % u = 0 at y = 0

    u(:, Ny) = Uinf; % u = Uinf at y = D

    v(i, j) = vn(i, j) - un(i, j) * dt / (2 * h) .* (vn(i+1, j) - vn(i-1, j)) ...

        - vn(i, j) * dt / (2 * h) .* (vn(i, j+1) - vn(i, j-1)) - dt / (2  * rho * h) * (p(i, j+1) - p(i, j-1)) ...

        + nu * (dt / h^2 * (vn(i+1, j) - 2 * vn(i, j) + vn(i-1, j)) + dt / h^2 * (vn(i, j+1) - 2 * vn(i, j) + vn(i, j-1))); % y-momentum

    v(1, :) = 0; % v = 0 at x = 0

    v(Nx, :) = 0; % v = 0 at x = L

    v(:, 1) = 0; % v = 0 at y = 0

    v(:, Ny) = 0; % v = 0 at y = D

end

%% post-processing

velocity_magnitude = sqrt(u.^2 + v.^2); % velocity magnitude

% Visualize velocity vectors and pressure contours

hold on

contourf(X, Y, velocity_magnitude', 64, 'LineColor', 'none'); % contour plot

set(gca,'FontSize',40)

% skip = 20;

% quiver(X(1:skip:end, 1:skip:end), Y(1:skip:end, 1:skip:end),... % Velocity vectors

%     u1(1:skip:end, 1:skip:end)', v1(1:skip:end, 1:skip:end)', 1, 'k','LineWidth', 0.1);

hh = streamslice(X, Y, u', v', 5); % Streamlines

set(hh, 'Color', 'k','LineWidth', 0.1);

colorbar; % Add color bar

colormap hsv % Set color map

axis equal % Set true scale

xlim([0 L]); % Set axis limits

ylim([0 D]);

xticks([0 L]) % Set ticks

yticks([0 D]) % Set ticks

clim([0 0.95*max(velocity_magnitude(:))]) % Legend limits

title('Velocity [m/s]');

xlabel('x [m]');

ylabel('y [m]');

Cite as: Fahad Butt (2023). Lid-Driven Cavity (https://fluiddynamicscomputer.blogspot.com/2023/09/lid-driven-cavity-matlab-code.html), Blogger. Retrieved Month Date, Year

Open - Source coded 3D Navier–Stokes equations in C++

     Here are 3D Navier–Stokes equations configured for lid-driven cavity flow. The syntax is C++. The results from this code are shown in Fig. 1. From the pressure-Poisson equation, I removed mixed derivative terms to improve solution stability. 😆

p[i][j][k] = ((p[i+1][j][k] + p[i-1][j][k]) * h * h + (p[i][j+1][k] + p[i][j-1][k]) * h * h + (p[i][j][k+1] + p[i][j][k-1]) * h * h) / (2 * (h * h + h * h + h * h)) - h * h * h * h * h * h / (2 * (h * h + h * h + h * h)) * (rho * (1 / dt * ((u(i+1, j, k) - u(i-1, j, k)) / (2 * h) + (v(i, j+1, k) - v(i, j-1, k)) / (2 * h) + (w(i, j, k+1) - w(i, j, k-1)) / (2 * h))));

p[0][j][k] = p[1][j][k];

p[num_i - 1][j][k] = p[num_i - 2][j][k];

p[i][0][k] = p[i][1][k];

p[i][num_j - 1][k] = p[i][num_j - 2][k];

p[i][j][0] = p[i][j][1];

p[i][j][num_k - 1] = 0.0;

u[i][j][k] = u[i][j][k] - u[i][j][k] * dt / (2 * h) * (u[i+1][j][k] - u[i-1][j][k]) - v[i][j][k] * dt / (2 * h) * (u[i][j+1][k] - u[i][j-1][k]) - w[i][j][k] * dt / (2 * h) * (u[i][j][k+1] - u[i][j][k-1]) - dt / (2 * rho * h) * (p[i+1][j][k] - p[i-1][j][k]) + nu * (dt / (h * h) * (u[i+1][j][k] - 2 * u[i][j][k] + u[i-1][j][k]) + dt / (h * h) * (u[i][j+1] [k] - 2 * u[i][j][k] + u[i][j-1][k]) + dt / (h * h) * (u[i][j][k+1] - 2 * u[i][j][k] + u[i][j][k-1]));

u[0][j][k] = 0.0;

u[num_i - 1][j][k] = 0.0;

u[i][0][k] = 0.0;

u[i][num_j - 1][k] = 0.0;

u[i][j][0] = 0.0;

u[i][j][num_k - 1] = -Uinf;

v[i][j][k] = v[i][j][k] - u[i][j][k] * dt / (2 * h) * (v[i+1][j][k] - v[i-1][j][k]) - v[i][j][k] * dt / (2 * h) * (v[i][j+1][k] - v[i][j-1][k]) - w[i][j][k] * dt / (2 * h) * (v[i][j][k+1] - v[i][j][k-1]) - dt / (2 * rho * h) * (p[i][j+1][k] - p[i][j-1][k]) + nu * (dt / (h * h) * (v[i+1][j][k] - 2 * v[i][j][k] + v[i-1][j][k]) + dt / (h * h) * (v[i][j+1][k] - 2 * v[i][j][k] + v[i][j-1][k]) + dt / (h * h) * (v[i][j][k+1] - 2 * v[i][j][k] + v[i][j][k-1]));

v[0][j][k] = 0.0;

v[num_i - 1][j][k] = 0.0;

v[i][0][k] = 0.0;

v[i][num_j - 1][k] = 0.0;

v[i][j][0] = 0.0;

v[i][j][num_k - 1] = 0.0;

w[i][j][k] = w[i][j][k] - u[i][j][k] * dt / (2 * h) * (w[i+1][j][k] - w[i-1][j][k]) - v[i][j][k] * dt / (2 * h) * (w[i][j+1][k] - w[i][j-1][k]) - w[i][j][k] * dt / (2 * h) * (w[i][j][k+1] - w[i][j][k-1]) - dt / (2 * rho * h) * (p[i][j][k+1] - p[i][j][k-1]) + nu * (dt / (h * h) * (w[i+1][j][k] - 2 * w[i][j][k] + w[i-1][j][k]) + dt / (h * h) * (w[i][j+1][k] - 2 * w[i][j][k] + w[i][j-1][k]) + dt / (h * h) * (w[i][j][k+1] - 2 * w[i][j][k] + w[i][j][k-1]));

w[0][j][k] = 0.0;

w[num_i - 1][j][k] = 0.0;

w[i][0][k] = 0.0;

w[i][num_j - 1][k] = 0.0;

w[i][j][0] = 0.0;

w[i][j][num_k - 1] = 0.0;

Fig. 1, Velocity and pressure iso-surfaces

     Of course constants need to be defined, such as grid spacing in space and time, density, kinematic viscosity. These equations have been validated, as you might have read and here already! Happy coding!

     If you want to hire me as your PhD student in the research projects related to turbo-machinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

Heaving Flat Plate Computational Fluid Dynamics (CFD) Simulation (In-House CFD Code)

     The adventure 🏕 🚵 I started a while ago to make my own CFD 🌬 code / software 💻 for my digital CV and to make a shiny new turbulence mode (one-day perhaps) is going along nicely. This post is about a 2D 10% flat plate undergoing forced heaving motion. Heaving motion is achieved by Eqn. 1.

h_y = H_o*sin(2*π*f_h*t)                                              Eqn. 1

     w.r.t. Eqn. 1 reduced frequency is defined as (f_h*H_o/U∞)), f_h is the frequency of oscillations, while t is the instantaneous time. H_o is the heaving amplitude and U∞ is the free stream velocity. h_y is the position of the flat plate. The animation is shown in Fig. 1.

     The Strouhal number is 0.228 and Reynolds number is at 500. As we can see, the in-house CFD code works very well for this complex CFD simulation. Validation of this work will never be completed 😆. As soon, I will move on to the next project without completing this one. Anyway, discretized Navier-Stokes equations are available here, in both C++ and MATLAB formats if you want to validate this non sense yourself! Good Luck!

The animation from in-house CFD simulation

     If you want to hire me as your PhD student in the research projects related to turbo-machinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

Home-Made Computational Fluid Dynamics (CFD), (Includes CFD code)

     One fine morning, I decided to code the Navier–Stokes equations using the finite-difference method. This post has the results of this adventure 🏞️ (so-far). As is customary with all my CFD work using commercial and home-made CFD codes, this too is inspired by the free lectures of Dr. Lorena Barba.

Internal / External Fluid Dynamics - Lid-Driven Cavity

     I started using Lid-Driven Cavity example. Just because everyone else uses it to validate the code they write. The lid-driven cavity case is giving correct results up to ~Re 1,000 without any turbulence models or wall functions.

     The results shown in Fig. 1 correspond to at Re 1,000. It can be seen that the present code which is just vanilla Navier–Stokes; captures the vortices at both bottom edges well as compared to published data. But why would you want to use a vanilla CFD code which solved only ~Re 1,000?, anyways... Some validation... u-velocity at (0.5, 0.1719) is at -0.3869 m/s as compared to -0.38289 m/s [1]. Meanwhile, v-velocity at (0.2266, 0.5) is 0.3252 as compared to 0.30203. Furthermore, v-velocity at (0.8594, 0.5) is at -0.4128 m/s as compared to -0.44993 m/s. All in all, a good agreement with published data. An animation is uploaded here.  Fig. 1a shows a different variation of lid driven cavity.

Fig. 1, Lid driven cavity post-processing

Fig. 1a, Heated lid driven cavity post-processing

Internal / External Aerodynamics - Backward-Facing Step (BFS)

     The next example is the case of Backward Facing Step (BFS). This case is analogues to flow around a building or a pipe with different diameter or anything with sharp edge along the flow. Let's be honest, flow around a building 🏢 is at much higher Re than the measly ~1,000 this code solves. But for validation, here we are.

 

     The cases are self validating as the re-attachment length of the vortex just after the step is very well documented [4]. Fig. 2 compares the results from the vanilla code with published data. It can be seen that the results are in good agreement.

Fig. 2, Backward Facing Step (BFS) post processing at Re 200

Internal Fluid Dynamics - Aero-Thermal Pipe Flow

     This is the first case I ran as it is self validating. The product of area and velocity should be same at  the outlet as given at the inlet as a boundary condition, which it is! Furthermore, density calculated should be proportional to the temperature calculated, which it is! The results at Re ~3300 are shown in Fig. 3. This case is relevant to bio medical field as Reynolds number in human blood 🩸flow is around same flow regime. Moreover, flow through heated pipes etc. is also relevant.

Fig. 3, Flow through heated pipe / between flat plates

Clot Flow

     Fig. 4 shows an example of flow where a clot appears in a blood 🩸 vessel. Done on this same 🏡made CFD 🌬️ code used throughout this blog. Or in a pipe where oil 🛢️ flows. It can be seen how smooth /less chaotic the flow is without clot (less work for heart ❤️ ). Eat 🦐 🥑 🥙 🥕 🥒 🍲 healthy! Mechanical properties of fluid are taken as follows. 5 ltr/min flowrate. Blood vessel diameter 40 mm. Fluid density at 1,070 Kg/m³ and dynamic viscosity at 4.5 cP.

Fig. 4, Flow inside a blood vessel with and without clot

External Aerodynamics - Flow around a Square Cylinder

     Why not circular cylinder you ask? I am very lazy 😂. It is hard work to draw a circle using code. I am used to CAD. 😇. Again, low Reynolds number are very rare in practical applications but for the sake of completeness, I added this case as well. Fig. 5 shows flow around the cylinder at Re 100. Vorticity is shown in Fig. 5. The results are compared with experimental data. The Strouhal number from the home-made CFD code is at 0.158 as compared to a value of 0.148. The results are within 7%  of published literature [2-3].

Fig. 5, Post-processing

Thank you for reading! If you want to hire me as your most awesome PhD student, please reach out!

Code:

     I present to you the code-able discretized equations for solving fluid flow problems. At first, C++ version is presented followed by the MATLAB version. Equation 1-2 are Poisson equations for pressure. Equation 3 and 4 are x-momentum equations (without source). Equations 5 and 6 are y-momentum equations.

double p_ij = ((pn(i + 1, j) + pn(i - 1, j)) * dy * dy + (pn(i, j + 1) + pn(i, j - 1)) * dx * dx) / (2 * (dx * dx + dy * dy)) - dx * dx * dy * dy / (2 * (dx * dx + dy * dy)) * (rho * (1 / dt * ((u(i + 1, j) - u(i - 1, j)) / (2 * dx) + (v(i, j + 1) - v(i, j - 1)) / (2 * dy)) - ((u(i + 1, j) - u(i - 1, j)) / (2 * dx)) * ((u(i + 1, j) - u(i - 1, j)) / (2 * dx)) - 2 * ((u(i, j + 1) - u(i, j - 1)) / (2 * dy) * (v(i + 1, j) - v(i - 1, j)) / (2 * dx)) - ((v(i, j + 1) - v(i, j - 1)) / (2 * dy)) * ((v(i, j + 1) - v(i, j - 1)) / (2 * dy))));                (1)

p(i, j) = ((pn(i+1, j) + pn(i-1, j)) * dy^2 + (pn(i, j+1) + pn(i, j-1)) * dx^2) ./ (2 * (dx^2 + dy^2)) - dx^2 * dy^2 / (2 * (dx^2 + dy^2)) * (rho * (1/dt * ((u(i+1, j) - u(i-1, j)) / (2 * dx) + (v(i, j+1) - v(i, j-1)) / (2 * dy)) - ((u(i+1, j) - u(i-1, j)) / (2 * dx)).^2 - 2 * ((u(i, j+1) - u(i, j-1)) / (2 * dy) .* (v(i+1, j) - v(i-1, j)) / (2 * dx)) - ((v(i, j+1) - v(i, j-1)) / (2 * dy)).^2));                (2)

double u_ij = un(i, j) - un(i, j) * dt / dx * (un(i, j) - un(i - 1, j)) - vn(i, j) * dt / dy * (un(i, j) - un(i, j - 1)) - dt / (2 * rho * dx) * (p(i + 1, j) - p(i - 1, j)) + nu * (dt / (dx * dx) * (un(i + 1, j) - 2 * un(i, j) + un(i - 1, j)) + dt / (dy * dy) * (un(i, j + 1) - 2 * un(i, j) + un(i, j - 1)));                (3)

u(i, j) = un(i, j) - un(i, j) * dt/dx .* (un(i, j) - un(i-1, j)) - vn(i, j) * dt/dy .* (un(i, j) - un(i, j-1)) - dt / (2 * rho * dx) * (p(i+1, j) - p(i-1, j)) + nu * (dt/dx^2 * (un(i+1, j) - 2 * un(i, j) + un(i-1, j)) + (dt/dy^2 * (un(i, j+1) - 2 * un(i, j) + un(i, j-1))));                (4)

double v_ij = vn(i, j) - un(i, j) * dt / dx * (vn(i, j) - vn(i - 1, j)) - vn(i, j) * dt / dy * (vn(i, j) - vn(i, j - 1)) - dt / (2 * rho * dy) * (p(i, j + 1) - p(i, j - 1)) + nu * (dt / (dx * dx) * (vn(i + 1, j) - 2 * vn(i, j) + vn(i - 1, j)) + dt / (dy * dy) * (vn(i, j + 1) - 2 * vn(i, j) + vn(i, j - 1)));                (5)

v(i, j) = vn(i, j) - un(i, j) * dt/dx .* (vn(i, j) - vn(i-1, j)) - vn(i, j) * dt/dy .* (vn(i, j) - vn(i, j-1)) - dt / (2 * rho * dy) * (p(i, j+1) - p(i, j-1)) + nu * (dt/dx^2 * (vn(i+1, j) - 2 * vn(i, j) + vn(i-1, j)) + (dt/dy^2 * (vn(i, j+1) - 2 * vn(i, j) + vn(i, j-1))));                (6)
    

     Of course constants need to be defined, such as grid spacing in space and time, density, kinematic viscosity. These equations have been validated, as you might have read already! Happy coding!

     If you want to hire me as your PhD student in the research projects related to turbo-machinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

References

[1] U Ghia, K.N Ghia, C.T Shin, "High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method", Journal of Computational Physics, Volume 48, Issue 3, 1982, Pages 387-411, ISSN 0021-9991, https://doi.org/10.1016/0021-9991(82)90058-4
[2] Khademinejad, Taha & Talebizadeh Sardari, Pouyan & Rahimzadeh, Hassan. (2015). Numerical Study of Unsteady Flow around a Square Cylinder in Compare with Circular Cylinder
[3] Ávila, Ítalo & Santos, Gabriel & Ribeiro Neto, Hélio & Neto, Aristeu. (2019). Physical Mathematical and Computational Modeling of the Two-Dimensional Flow Over a Heated Porous Square Cylinder. 10.26678/ABCM.COBEM2019.COB2019-0854
[4] Irisarri, Diego & Hauke, Guillermo. (2019). Stabilized virtual element methods for the unsteady incompressible Navier–Stokes equations. Calcolo. 56. 10.1007/s10092-019-0332-5. 

Datacenter Visualization (Verified and Validated)

     This simulation is done to create an aero-thermal digital twin of a datacenter using CFD. The details of datacenter are taken from [1]. The datacenter CAD model is shown in Fig. 1.

     The simulation employs κ − ε turbulence model with damping functions, SIMPLE-R (modified), as the numerical algorithm and second-order upwind and central approximations as the spatial discretization schemes for the convective fluxes and diffusive terms. The time derivatives are approximated with an implicit first-order Euler scheme. Flow simulation solves the Navier–Stokes equations, which are formulations of mass, momentum, and energy conservation laws for fluid flows. To predict turbulent flows, the Favre-averaged Navier–Stokes equations are used.

Fig. 1, Datacenter CAD

     A Cartesian mesh with octree refinement, cut-cell method and immersed boundary methods is used. Special mesh refinements are deployed in the areas of interest i.e. inlets and outlets and sharp edges of server racks and CRAH units to accurately capture aero-thermal gradients and vortices. The resulting computational mesh has 2,698,156 cells. The computational domain and mesh are shown in Fig. 2.

Fig. 2, Computational mesh and domain

     The results from the numerical analysis were compared with [1]. The results are in excellent agreement with previously published data. The animation in Fig. 3 shows thermal distribution inside datacenter using cut-plots. Fluid velocity distribution is also shown. The cut-plots are superimposed with streamlines and velocity vectors. These post processing features help identify hot-spots and recirculation zones. Design improvements can be made to reduce thee unwanted flow features. Within Fig. 3, Flow trajectories colored by air temperature are also shown. These show path the fluid takes between various inlets and outlets in the datacenter such as the CRAH system supply and return zones and inlets and outlets of servers. Fig. 4 shows various post processing tools available for diagnosing various issues from the aero-thermal perspective. These include iso-surfaces, cut-plots, flow trajectories and surface plots etc.

Fig. 3, The post processing animations

Fig. 4, The post processing images

     A comparison of the results from present simulations with previously published literature is shown in Fig. 5. it can be seen that out results are in close agreement with previously published numerical and experimental data [1]! The locations at which the data is extracted is shown in Fig. 6. Within Fig. 5, solid lines indicate present study, dashed lines indicate published numerical results and filled circles represent published experimental results.

Fig. 5, Comparison of results

Fig. 6, Location of data extraction

     If you want to collaborate on the research projects related to turbomachinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

References

[1] Wibron, Emelie, Anna-Lena Ljung, and T. Staffan Lundström. 2018. "Computational Fluid Dynamics Modeling and Validating Experiments of Airflow in a Data Center" Energies 11, no. 3: 644. https://doi.org/10.3390/en11030644

DARPA Suboff Submarine CFD Simulation (Backed-up by Water Tunnel Data); Update 01: Water Level Simulations (Volume of Fluid)

      This post is about the CFD analysis of the DARPA Suboff at various speeds. The DARPA Suboff model is based on a generic submarine. The submarine geometry is shown in Fig. 1. The submarine in fully submerged in water. Refer to section "Update 01" for submarine sailing at the water level. The submarine geometry is available here. The geometry is made using equations from [1]. Machine Learning available here.

Fig. 1, DARPA Suboff submarine CAD

      The simulations are validated with published literature [2-4]. SolidWorks Flow Simulation Premium software is employed for the simulations. Fig. 2 shows results of drag force at various cruise speeds. It can be seen that the results are in close agreement with the published experimental / numerical data.

Fig. 2, Comparison of results

     The mesh has 656,714 cells in total. With 34,258 cells on the submarine surface. Special mesh refinements are added in the regions of interest i.e. regions with high gradients, the wake and on the control surfaces of the submarine. The mesh for is shown in Fig. 3.

Fig. 3, The mesh and computational domain

     The results from the CFD post processing are presented next. Velocity cut-plots showing velocity distribution and wake of the submarine, surface pressure distribution on the submarine and vorticity around and in the wake of the submarine are also shown in Fig. 4. Within Fig. 4, the black arrows represent the direction of on coming flow.

Fig. 4, The post processing

     Thank you for reading, If you would like to collaborate on projects, please reach out.

Update 01

     This section presents the results from a simulation using volume of fluid method. The water line is just below the sail of the submarine, as shown in Fig. 5. The iso-surfaces are shown in Fig. 6. The black arrows in both Figs. 5-6 represent direction of on coming flow. These simulations allow for visualization of wake of submarines and then methods can be device to reduce the wake.

Fig. 5, The color blue represents water and white represents air

Fig. 6, 3D wake of a sailing submarine

References

[1] Groves, Nancy C. Huang, Thomas T. Chang, Ming S., "Geometric Characteristics of DARPA (Defense Advanced Research Projects Agency) SUBOFF Models (DTRC Model Numbers 5470 and 5471)",  David Taylor Research Center, Bethesda MD, Ship Hydromechanics Dept, ADA210642, 1989 https://apps.dtic.mil/sti/citations/ADA210642

[2] Yu-Hsien Lin and Xian-Chen Li, "The Investigation of a Sliding Mesh Model for Hydrodynamic Analysis of a SUBOFF Model in Turbulent Flow Fields", Journal of Marine Science and Engineering, 8(10), 744, https://doi.org/10.3390/jmse8100744

[3] Liu, H.-L.; Huang, T.T. Summary of DARPA SUBOFF Experimental Program Data; Naval Surface Warfare Center Carderock Division, Hydromechanics Directorate: West Bethesda, MD, USA, 1998.

[4] Özden, Y.A.; Özden, M.C.; Çelik, F. Numerical Investigation of Submarine Tail Form on the Hull Eciency. In Proceedings of the Fifth International Symposium on Marine Propulsors, Espoo, Finland, 12–15 June 2017

Heaving Airfoil Simulation

This post is about a 2D NACA 0012 heaving aerofoil. Heaving motion is achieved by changing the angle of attack on the aerofoil based on the Eqn. 1.

α_e = arctan[2*π*St_a*cos(2*π*f_h*t)]+ α_i               Eqn. 1

w.r.t. Eqn. 1, α_e is the effective angle of attack, St_a is Strouhal number, f_h is the heaving frequency.

The case S1 and H6 from [1] are compared in the animations below.

If you want to collaborate on the research projects related to turbomachinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

References

[1] https://doi.org/10.1121/10.0001419

Hypersonic Flow over a Two Dimensional Heated Cylinder

     This post is about the simulation of hypersonic flow over a heated circular cylinder, in two dimensions.

     Equation 1 is used as a relationship between Mach and the Reynold number.

M= Re*μ*√(R*T) ÷ d*P*√γ     (1)

     w.r.t. equation 1, the parameters represent the following quantities.

     M     Freestream Mach number at 17.6
     Re    Reynolds number at 376,000
     μ     Dynamic viscosity at 1.329045e^-5 Ns.m^-2
     R     Specific gas constant at 286.9 J.(kg.K)^-1
     T     Freestream temperature 200 K
     d     Cylinder diameter at 5.6730225e^-4 m
     P     Freestream pressure at 101325 Pa
     γ     Specific heat ratio at 1.4
     Tw  Wall temperature of cylinder at 500 K
     Pr    Prandtl number at 0.736

     The boundary conditions were taken from [1]. A comparison with [1] is shown in Fig. 1. Inside Fig. 1, the red dotted line with circles represents the data from [1]. The black solid line represents the data from the present simulation. Within Fig. 1, 0° represents the stagnation point. The velocity, pressure, Mach number and temperature contours are shown in Fig. 2.

Fig. 1 A comparison with previous research [1].

Fig. 2, Top Row, L-R: Velocity and pressure contours. Bottom Row, L-R: Mach number and temperature contours.

The computational mesh and the computational domain with boundary conditions visible are shown in Fig. 3-4, respectively. The computational domain had a size of 20D x 20D. The mesh had 836,580 total cells and 944 cells were located at the solid fluid boundary. Several local mesh controls were employed to capture the shockwave properly.

Fig. 3, The computational mesh.

Fig. 4, The computational domain.

     The solution method is Finite Volume method. SIMPLE-R is the solver employed. Implicit central difference scheme for diffusion terms, second-order Upwind scheme for convective terms and first-order implicit for temporal terms are used. The mesh created uses the Cartesian mesh with Immersed Boundary method.

     Reference:

     [1] https://arc.aiaa.org/doi/10.2514/6.2004-2371 

     Thank you for reading. If you would like to collaborate on research projects, please reach out. I am looking for a PhD position, any guidance would be appreciated.

Automotive Computational Fluid Dynamics (CFD) Analysis

     This post is about the numerical analysis of an Ahmed body. It was a new experience because of the area between the car's floor and the road which is different from most of the numerical analysis performed in open channel aeronautics and turbo-machinery.

     The numerical analysis was performed using the commercial software, SolidWorks Flow Simulation. The software employs κ − ε turbulence model with damping functions, SIMPLE-R (modified), as the numerical algorithm and second-order upwind and central approximations as the spatial discretization schemes for the convective fluxes and diffusive terms. The time derivatives are approximated with an implicit first-order Euler scheme. Flow simulation solves the Navier–Stokes equations, which are formulations of mass, momentum, and energy conservation laws for fluid flows. To predict turbulent flows, the Favre-averaged Navier–Stokes equations are used.

     The software generates Cartesian mesh using immersed boundary method. The mesh had a cell size of 0.035 m in the far field regions within the computational domain. A fine mesh was need between the road the the car's floor to make sure the interaction of the car's floor with the road was captured accurately. Therefore the mesh between the car's floor and the road was refined to have a cell size of 0.00875 m. Another mesh control was applied around the body to refine the mesh with a cells size of 0.0175 m to capture the trailing vortices. The resulting mesh had 209,580 total cells, among those cells, 31,783 cells were at the solid fluid boundary. The computational domain size was ~1L x 1.12L x 3L where L being the vehicle's length. The computational domain along with the computational mesh is shown in Fig. 1.

Fig. 1 Mesh, computational domain and the boundary conditions.

     The red arrows within the Fig. 1 represents the inlet boundary condition of ambient (free-stream) velocity and the blue arrows represent the outlet boundary condition of the ambient pressure. The green arrows represents the co-ordinates axes direction.

     The results from the numerical analysis were compared with [1-3]. The results are within 10% of the experimental results. The velocity (superimposed by the velocity streamlines) and pressure profiles around the car body at various free-stream velocities is shown in Fig. 2.

Fig. 3 Velocity and pressure plots. From the top, Row 1, L-R; ambient velocity of 30 and 40 m/s. Row 2, L-R; ambient velocity of 60 and 80 m/s. Row 3, ambient velocity of 105 m/s.

     It was a good experience learning about automotive CFD after spending a long time in aeronautic/turbo-machinery CFD. Thank you for reading. Please share my work. If you would like to collaborate on a project please reach out.

[1] F.J.Bello-Millán, T.Mäkelä, L.Parras, C.delPino, C.Ferrera, "Experimental study on Ahmed's body drag coefficient for different yaw angles", Journal of Wind Engineering and Industrial Aerodynamics, Volume 157, October 2016, Pages 140-144.

[2] Guilmineau E., Deng G.B., Queutey P., Visonneau M. (2018) Assessment of Hybrid LES Formulations for Flow Simulation Around the Ahmed Body. In: Deville M. et al. (eds) Turbulence and Interactions. TI 2015. Notes on Numerical Fluid Mechanics and Multidisciplinary Design, vol 135. Springer, Cham.

[3] A. Thacker, S.Aubrun, A.Leroy, P.Devinant, "Effects of suppressing the 3D separation on the rear slant on the flow structures around an Ahmed body", Journal of Wind Engineering and Industrial Aerodynamics, Volumes 107–108, August–September 2012, Pages 237-243.

Computational Fluid Dynamics Analysis of a Symmetrical Wing, Update 01

     This post is about the computational fluid dynamics analysis of a wing. The wing analyzed employed the NACA 0021 section throughout. The wing had a span of 4 m and a chord length of 1 m. The Reynolds number was kept at 3,000,000. The software employed was SolidWorks Flow Simulation Premium.

     The mesh had a total of 385,064 cells of which 84,826 cells were in contact with the wing surface, as shown in Fig. 1. The results are, indeed, mesh independent. Mesh controls were employed to refine the mesh near the wing surface. The computational domain employed was of cylindrical shape.

 

Fig. 1, The computational mesh around the wing.

 

     The velocity variation at various angles of attack around the wing cross-section is shown in Fig. 3 while the pressure variation on the wing surface is shown in Fig. 4. The results were validated against experiments conducted by [1].

 

Fig. 2, Velocity variation around the wing at 0-25 degree AOA, 5 degree increments.

 

Fig. 3, Pressure variation at the wing surface at 0-25 degree AOA, 5 degree increments.

     The purpose of this blog is maintain my online portfolio. I did this analysis because I realized I haven't written anything of this nature before. All of my previous simulations and/or blog entries were from the propulsion, renewable energy and turbo-machinery areas.
 

     Update 01

     CAD files are available here.

 

    
     Thank you for reading. If you would like to collaborate on research projects, please feel free to contact.

     [1] Fernando A. Rocha, Adson A. de Paula, Marcos d. Sousa, André V. Cavalieri, and Vitor G. Kleine, "Lift enhancement by wavy leading edges at Reynolds numbers between 700,000 and 3,000,000," Proceedings of the 2018 Applied Aerodynamics Conference, AIAA AVIATION Forum, Atlanta, GA, 2018.

     https://doi.org/10.2514/6.2018-3816

1D Transient Diffusion (MATLAB code)

clear; clc;% clear the screen and memory

Nx=500; %number of space nodes

Nt=10000; %number of time nodes

Lx=0.3; %length of space (m)

Lt=10; %length of physical time (s)

dx=Lx/(Nx-1); %grid spacing

dt=Lt/(Nt-1); %time step

c=0.000097; %speed of wave (constant)

a=c*dt/dx.^2;

u=100*ones(Nx,1); %initialization of matrix for the property under investigation

x=zeros(Nx,1); %initialization of space

for i=1:Nx-1 %space loop

    x(i+1)=x(i)+dx;

end

u(1)=-20; %boundary condition

u(Nx)=50; %boundary condition

for t=0:dt:Lt %time loop

    un=u; %u(i,t)

    for i=2:Nx-1 %solution loop, backward in space forward in time

        u(i)=((1-2*a)*un(i))+a*(un(i+1)+un(i-1)); %discretized equation, u(i,t+1)

    end

    plot(x,u) %plotting

    title({t;'Time Elaplsed (s)'},'FontSize',38.5);

    pause(0.001)

end

 

 

1D non-Linear Convection (MATLAB code)

clear; clc;% clear the screen and memory

Nx=41; %number of space nodes

Nt=82; %number of time nodes

Lx=2; %length of space (m)

Lt=1; %length of time (s)

dx=Lx/(Nx-1); %grid spacing

dt=Lt/(Nt-1); %time step

a=dt/dx;

u=ones(Nx,1); %initialization of solution matrix

x=zeros(Nx,1); %initialization of space

hold on

for i=1:Nx-1 %space loop

    x(i+1)=x(i)+dx;

end

for i=0.5/dx:1/dx %initial conditions

    u(i)=2;

end

for t=1:Nt %time loop

    un=u; %u(i,t)

    for i=2:Nx %solution loop, backward in space forward in time

        u(i)=un(i)-u(i)*a*(un(i)-un(i-1)); %discretized equation, u(i,t+1)

    end

    plot(x,u,'k') %plotting

    pause(0.1)

end