Tutorial 22 – Automatic Shape Optimization

Introduction

In this tutorial we are going to show, how to define and perform automatic shape optimization in QuickerSim CFD Toolbox. The problem consists of two main steps:

In this tutorial we are going to use the case from Tutorial 13 as our basic flow problem. Our intention will be to deflect trailing edge of the NACA 0012 airfoil to maximize its lift to drag ratio for a given Reynolds number and angle of attack. It is recommended that you first get familiar with the case and simulation code from Tutorial 13 (see here for details: http://www.quickersim.com/cfd-toolbox-for-matlab/tutorials). Geometry parametrization that we are going to introduce allows to deflect the trailing edge of the airfoil in a similar manner to the one shown below.

Defining Optimization Problem

To define automatic optimization we need to choose one of the MATLAB’s optimization algorithms. We will use fminbnd function in this example. It can optimize only one design parameter, but it works in a bounded interval and is available nativily in MATLAB. Hence, each MATLAB user will be able to run the example. Let us comment on two additional requirements. Having an optimization algorithm which works in a bounded design space is always recommended. Otherwise, huge changes to design variables which an unbounded optimizer could generate may lead to mesh destruction and other enormous, unphysical changes to original geometry. Users having access to the Optimization Toolbox can choose more advanced optimization algorithms and perform multivariable optimization. However, the building blocks are always the same, so we believe that this tutorial is generic enought to allow each user to adjust this tutorial to his applications with Matlab Optimization Toolbox examples at hand.

Once optimization algorithm has been chosen, we can now proceed to modifications that we need to make to our CFD solver. The modifications are as follows:

Thus, we are now going to modify original code of Tutorial 13 and introduce all above mentioned additional adjustments. The whole code presented below shall be place in a separate file (cfdSolver.m) and modifications with respect to original code of Tutorial 13 are clearly marked with comments.

% MODIFICATION 1: Add function header with input arguments
function f = cfdSolver(optParameter)

[p,e,t] = importMeshGmsh('naca0012.msh');

% Add boundary layer mesh (refer to Tutorial 5 for boundary layer meshing)
layers.h1 = 0.001;
layers.q = 1.3;
layers.n = 5;
wallIds = [9];
[p,e,t] = extrudeLayers2D(p,e,t,wallIds,layers);
% MODIFICATION 2: Add code to parametrize and deform geometry
% (See Tutorial 6 for details of mesh deformation functions)
% Deform mesh according to optParameter set by optimizer
rbf = [0.25 1.2 0]';
[TM, pinit] = initDeformMeshRBF(p,rbf);
p = deformMeshRBF(pinit, TM, [0.0 optParameter]);
% Convert mesh to second order and display mesh
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);
% displayMesh2D(p,t);
% Define fluid (air kinematic viscosity in m^2/s)
nu = 1.5e-5;
% Init solution, define convergence and under-relaxation parameters
[u, convergence] = initSolution(p,t,[0 0],0);
maxres = 1e-6;
maxiter = 25;
alpha = 0.6;
% Compute wall distance and duplicate velocity vector
d = computeWallDistance(p,e,t,wallIds);
us = u;
% Define inlet velocity profile
vinf = 1.5; % 1.5 m/s yields Re = 100 000 at airfoil chord of 1 m.
aoa = 5; % angle of attack in degrees
aoa = aoa/180*pi; % calculate angle of attack in radians
vel = vinf*[cos(aoa) sin(aoa)]; % define velocity vector
% Iterate nonlinear terms
for iter = 1:maxiter
    % Compute turbulent viscosity (turbulent intensity of 4% and Lt = 0.5m)
    nut = turbulentViscosity(p,t,u,d,nu,indices,'CITM', 0.04, 0.5);
% Under-relax velocity and assemble matrix and right hand side
     u = alpha*u+(1-alpha)*us;
     [NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu+nut,u(indices.indu),u(indices.indv),'supgDoublyAsymptotic');
% Apply boundary conditions
     [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,7,'inlet', vel);
     [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,9,'wall');
% Compute and plot residuals
     [stop, convergence] = computeResiduals(NS,F,u,size(p),convergence,maxres);
     plotResiduals(convergence,2);
% Break if solution converged
     if(stop)
         break;
     end
% Solve equations
     us = u;
     u = NS\F;
 end
% Compute lift and drag force coefficients

% Compute force at the airfoil (Physical Line id = 9)
Force = computeForce(p,e,t,u,nu,9);

% Project the force to direction tangential and perpendicular to freestream
% velocity (drag and lift)
D = Force*[cos(aoa) sin(aoa)]';
L = Force*[-sin(aoa) cos(aoa)]';

% Compute drag and lift coeeficients by dividing the force by dynamic
% pressure (unit fluid density is assumed for incompressible flow
% formulation)
Cd = D/(0.5*norm(vel)^2)
Cl = L/(0.5*norm(vel)^2)

% MODIFICATION 3: Return objective function value f and display current design
f = -Cl/Cd; % return negative value (optimizer minimizes the function, where we want to maximize L/D)
disp(['Current design parameter = ' num2str(optParameter) ', L/D = ' num2str(-f)]);

% OPTIONAL MODIFICATION: Export flow field results from each trial
pressure = generatePressureData(u,p,t);
exportToGmsh2D(['pressure_opt_' num2str(optParameter) '.msh'],pressure,p,t,'Pressure');
exportToGmsh2D(['velocity_opt_' num2str(optParameter) '.msh'],u(indices.indu),p,t,'vx');

% OPTIONAL MODIFICATION: Add current design to a file with all designs
designs = fopen('alldesigns.txt','at');
fprintf(designs,'%f\t%f\n',optParameter,Cl/Cd);
fclose(designs);

% MODIFICATION 4: Add end of the function
end

Once we have converted our solver script to a Matlab function, it is nice to call this function from Matlab command line for undeformed geometry (with zero argument in our case) to check, if everything reflects original simulation. To run the solver type:

cfdSolver(0)

Setup of Optimization Algorithm

The final step is to write a simple Matlab script, which will start and run the whole automatic optimization procedure. We need to set the following settings:

The simple code of optimization script is provided below

% Initial alpha value
alpha = 0.0;

% Lower and upper bound on alpha during optimization
alphamin = -0.2;
alphamax = 0.1;

% Run optimization algorithm
x = fminbnd(@(alpha) cfdSolver(alpha),-0.2,0.1,optimset('MaxFunEval',50));

Results

In this case optimization procedure converges to optimal solution within just 10 simulation runs. The L/D ratio is improved from value of 7.43 for undeformed geometry to L/D = 9.14 for optimized geometry, yielding improvement of almost 23%. Figure below illustrates convergence of the optimization process.

The two following figures show original (undeformed) and optimized geometry.