High Lift Common Research Model (CRM-HL) CFD Simulation (with validation from Stanford CTR and NASA) Update 01: Formation Flight

     This post is about the CFD analysis of the nominal (2A) configuration of the Boeing / NASA CRM-HL with flap and slat angles at 40 and 37 degrees. The configuration is shown in Fig. 1.

Fig. 1, CRM-HL CAD

     The dimensions are mentioned in [1] and is available here and here. The numerical simulations are validated with published literature [3, 4]. SolidWorks Flow Simulation Premium software is employed for the CFD simulations. The flight conditions are Mach 0.2 at 289.444 K and 170093.66 Pa [2].

     Fig. 2 shows the computational mesh along with the computational domain. The surface mesh is also shown with Fig. 2. It can be seen that the mesh is refined in the areas of interests and in the wake of the aircraft to properly capture the relevant flow features. The Cartesian mesh with immersed boundary method, cut-cell approach and octree refinement is used for creating the mesh. The mesh has 6.5 million cells, of which around 0.71 millions cells are at boundary of the aircraft.

     The simulations employ κ-ε turbulence model with damping functions and two-scales wall functions. SIMPLE-R (modified) as the numerical algorithm. Second order upwind and central approximations as the spatial discretization schemes for the convective fluxes and diffusive terms. The software solves the Favre-averaged Navier-Stokes equations to predict turbulent flow. The simulations are performed to predict three-dimensional steady-state flow over the aircraft

Fig. 2, Dashed lines indicate symmetry boundary conditions

     Angle of attack of 7 and 17 degrees are considered. At 7 degree angle of attack, the drag life and moment coefficients are within 6, 8 and 12% of the published data, respectively. For 17 degree angle of attack, the coefficients are within 4, 7, 33% of the published data [3, 4]. On average, the results of the present simulations are within 6% of the published data for force coefficients and within 12% in terms of pitching moment coefficient. The pitching moment coefficient will improve with refinement of the mesh, as shown by [3]; which I will do if ever I convert this into a manuscript 😂. The post processing from CFD simulations is shown in Fig. 3-4. Within Fig. 3, the iso-surfaces represent vorticity in the direction of flow, colored by pressure. The direction of flow is shown by black arrows. Streamlines colored by vorticity are also visible. It can be seen from Fig. 2 that the simulation captures important flow features such as vortex formation very accurately, in such small number of mesh cells. Fig. 4 shows wing of the velocity iso-surfaces colored by vorticity in the direction of flow, focused around the wing.

Fig. 3, Alternating 7 and 17 degree angles of attacks

Fig. 4, Top-Bottom, 17 and 7 degrees angle of attack

     The simulations are solved using 10 of 12 threads of a 4.0 GHz CPU with 32 GB of total system memory of which almost 30 GB remains in use while the simulations are in progress. To solve each angle of attack, 3 hours and 47 minutes are required.

Update 01

     Formation flight simulation is performed with 7.3 million cells (limited by RAM). Symmetry boundary condition is used to simulate only the area of interest. The mesh and computational domain is shown in Fig. 5. Within Fig. 5, dashed lines indicate symmetry boundary condition.

Fig. 5, Airliners in formation flight mesh and computational domain

     The results show that the leading aircraft has 16.52% more drag than the trailing aircraft. The results also show that the the leading aircraft has 5.4% less lift than the trailing aircraft. However, the leading aircraft produces 2.4% more lift than the airliner flying without the trailing aircraft (flying solo). The leading aircraft also produces 3% less drag than the airliner flying solo. These results are at 7 degree angle of attack. Whether these results are fruitful aero-structurally, remains to be seen. The post processing from simulations are shown in Fig, 6.

Fig. 6, Top-Bottom, 17 and 7 degrees angle of attack

     Within Fig. 7, surface plots represent pressure, velocity streamlines indicate vortices, iso-surfaces show vorticity magnitude in the direction of flight.

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

References

[1] Doug Lacy and Adam M. Clark, "Definition of Initial Landing and Takeoff Reference Configurations for the High Lift Common Research Model (CRM-HL)", AIAA Aviation Forum, AIAA 2020-2771, 2020 10.2514/6.2020-2771

[2] 4th AIAA CFD High Lift Prediction Workshop Official Test Cases, https://hiliftpw.larc.nasa.gov/Workshop4/OfficialTestCases-HiLiftPW-4-2021_v15.pdf

[3] K. Goc, S., T. Bose and P. Moin, "Large-eddy simulation of the NASA high-lift common research model", Center for Turbulence Research, Stanford University, Annual Research Briefs 2021

[4] 4TH High Lift Workshop Results, ADS, https://new.aerodynamic-solutions.com/news/18

1D Laplace's equation using Finite Difference Method

This post is about FDM for Laplace Equation with various boundary conditions.

MATLAB Code (1D, Dirichlet Boundary Conditions)

%% initialize the workspace, clear the command window

clear; clc

%% finite difference 1D laplace dirichlet boundary conditions %% d2u/dx2 = 0 %% u(o) = 10, u(L) = 4 %% Ax=b%%

N = 4 ; %number of grid points

L = 1; %length of domain

dx = L/(N-1); %element size

%% initialize variables %%

l = linspace(0,L,N); %independent

u=zeros(1,N); %dependent

%% boundary conditions %%

u(1)=10;

u(end)=4;

%% b vector %%

b=zeros(N-2,1);

b(1) = b(1) - u(1);

b(end) = b(end) - u(end);

%% A matrix

A = -2*eye(N-2,N-2);

for i=1:N-2

    if i<N-2

        A(i,i+1)=1;

    end

    if i>1

        A(i,i-1)=1;

    end

end

%% solve for unknowns %%

x = A\b;

%% fill the u vector with unknowns %%

u(2:end-1) = x;

%% plot the results %%

hold on; grid on; box on, grid minor

set(gca,'FontSize',40)

set(gca, 'FontName', 'Times New Roman')

ylabel('u','FontSize',44)

xlabel('l','FontSize',44)

plot(l,u,'-o','color',[0 0 0],'LineWidth',2,'MarkerSize',20)

MATLAB Code (1D, Mixed Boundary Conditions)

%% initialize the workspace, clear the command window

clear; clc

%% finite difference 1D laplace mixed boundary conditions %% d2u/dx2 = 0 %% u(o) = 10, du/dx(L) = 4 %% Ax=b%%

N = 5; %number of grid points

L = 1; %length of domain

dx = L/(N-1); %element size

a = 4;

%% initialize variables %%

l = linspace(0,L,N); %independent

u=zeros(1,N); %dependent

%% dirichlet boundary condition %%

u(1) = 10;

%% b vector %%

b=zeros(N-1,1);

b(1) = b(1) - u(1);

b(end) = b(end) + dx*a; %neumann boundary condition added to b vector

%% A matrix

A = -2*eye(N-1,N-1);

for i=1:N-1

    if i<N-1

        A(i,i+1)=1;

    end

    if i>1

        A(i,i-1)=1;

    end

end

A(N-1,N-1) = -1; %neumann boundary condition added to A matrix

%% solve for unknowns %%

x = A\b;

% fill the u vector with unknowns %%

u(2:end) = x;

%% plot the results %%

hold on; grid on; box on, grid minor

set(gca,'FontSize',40)

set(gca, 'FontName', 'Times New Roman')

ylabel('u','FontSize',44)

xlabel('l','FontSize',44)

plot(l,u','-o','color',[0 0 0],'LineWidth',2,'MarkerSize',20)

Computational Fluid Dynamics Simulation of a Swimming Fish (Includes UDF)

      This post is about the simulation of a swimming fish. The fish body is made of NACA 0020 and 0015 aero-foils (air-foils). The fluke is made of NACA 0025 aero-foil (air-foil), as shown in Fig. 1. the CAD files with computational domain modelled around the fish is available here.

Fig. 1, The generic fish CAD model

      The motion of the fish's body is achieved using a combination of two user-defined functions (UDF). The UDFs use DEFINE_GRID_MOTION script mentioned below, for the head/front portion. This is taken from the ANSYS Fluent software manual, available in its original form here. The original UDF is modified for present use as required. To move the mesh, dynamic mesh option within ANSYS Fluent is enabled; with smoothing and re-meshing options. The period of oscillation is kept at 2.0 s. The Reynolds number of flow is kept at 100,000; which is typical for a swimming fish.

/**********************************************************

 node motion based on simple beam deflection equation

 compiled UDF

 **********************************************************/

#include "udf.h"

DEFINE_GRID_MOTION(undulating_head,domain,dt,time,dtime)

{

  Thread *tf = DT_THREAD(dt);

  face_t f;

  Node *v;

  real NV_VEC(omega), NV_VEC(axis), NV_VEC(dx);

  real NV_VEC(origin), NV_VEC(rvec);

  real sign;

  int n;

  

  /* set deforming flag on adjacent cell zone */

  SET_DEFORMING_THREAD_FLAG(THREAD_T0(tf));

  sign = 0.15707963267948966192313216916398 * cos (3.1415926535897932384626433832795 * time);

  

  Message ("time = %f, omega = %f\n", time, sign);

  

  NV_S(omega, =, 0.0);

  NV_D(axis, =, 0.0, 1.0, 0.0);

  NV_D(origin, =, 0.7, 0.0, 0.0);

  

  begin_f_loop(f,tf)

    {

      f_node_loop(f,tf,n)

        {

          v = F_NODE(f,tf,n);

          /* update node if x position is greater than 0.02

             and that the current node has not been previously

             visited when looping through previous faces */

          if (NODE_X(v) > 0.05 && NODE_X(v) < 0.7 && NODE_POS_NEED_UPDATE (v))

            {

              /* indicate that node position has been update

                 so that it's not updated more than once */

              NODE_POS_UPDATED(v);

              omega[1] = sign * pow (NODE_X(v), 0.5);

              NV_VV(rvec, =, NODE_COORD(v), -, origin);

              NV_CROSS(dx, omega, rvec);

              NV_S(dx, *=, dtime);

              NV_V(NODE_COORD(v), +=, dx);

            }

        }

    }

  end_f_loop(f,tf);

}

      The computational mesh, as shown in Fig. 2, uses cut-cell method with inflation layers. The mesh has 2,633,133 cells. The near wall y+ is kept at 5. The Spalart-Allmaras turbulence model is used to model the turbulence. The second order upwind scheme is used to discretize the momentum and modified turbulent viscosity equations. The time-step for this study is kept at 100th/period of oscillation.

Fig. 2, The mesh and zoom in view of the trailing edge.

      The animation showing fish motion is shown in Fig. 3. Within Fig. 3, the left side showcases the velocity iso-surfaces coloured by pressure and the vorticity iso-surfaces coloured by velocity magnitude is shown on the right.

Fig. 3, The animation.

      Another animation showing the fish motion is shown in Fig.4. Within Fig. 4, the left side shows surface pressure while the right side shows pressure iso-surfaces coloured by vorticity.

Fig. 4, The animation.

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

Flapping Aerofoil For Propulsion

     This post is about a 2D NACA 0012 aerofoil undergoing forced flapping motion for propulsion purposes. Heaving motion is achieved by applying a vertical velocity on the aerofoil based on the Eqn. 1. Similarly the pitching motion is achieved by applying a rotational velocity, governed by Eqn. 2.

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

ω = -2*π*f_h*ϑ*sin[(2*π*f_h*t) + 1.5708]                               Eqn. 2

     w.r.t. Eqn. 1-2 reduced frequency is defined as (2*π*f_h*H_o/U∞)), f_h is the frequency of oscillations, while ω, t and ϑ_o represent rotational velocity, instantaneous time and maximum pitching angle. H_o is the heaving amplitude and U∞ is the free stream velocity.

     The flapping motion is achieved by a combination of the heaving and pitching. In this particular simulation, the aerofoil is in the propulsion mode, meaning the feathering parameter χ is less in magnitude than 1.0. Feathering parameter is defined by Eqn. 3.

χ = ϑ/arctan(h0*2*π*f_h/U∞)                                  Eqn. 3

     The boundary conditions employed for the simulation are at Re 1,000, K = 1.41, H_o = aerofoil chord length, χ = 0.5489 and f_h = 0.003391 Hz. The animation of the pressure, vorticity and velocity contours is shown in Fig. 1.

Fig. 1, Flow animation, fluid flow direction is from left to right.

     The results of present simulation are compared with [1]. In terms of maximum lift, a maximum deviation of 5% is observed as compared to [1], as shown in Fig. 2. The maximum lift coefficient for available data is ~4.224 while the maximum lift coefficient from the present simulation is ~4.057. The average thrust produced is within 2% of [1]. Average thrust coefficient per cycle from [1] is 0.9957 while the result from present simulation reveals the thrust coefficient to be 1.0098.

Fig. 2, A comparison of coefficient of lift.

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

References

[1] https://doi.org/10.1017/jfm.2017.508