Tutorial 22 – Automatic Shape Optimization

This tutorial is intended for the full version of the toolbox.

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:

  1. We need to prepare a script that is able to solve our flow problem.
  2. We need to use one of MATLAB’s optimization algorithms which will run our basic CFD solver multiple times.

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: https://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:

  1. The CFD solver for the base case must be a Matlab function, not a script. The function shall take design variables as input arguments.
  2. The CFD solver must return value of the objective function upon completion of computations.
  3. The CFD solver must be equipped with few additonal lines of code to deform the original geometry according to current design parameters.
  4. Optionally, it is often helpful to record all designs generated during optimization and archive their performance (objective function values and maybe also flow field results).

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:

  • initial value of the design parameter,
  • lower and upper bounds on the design parameter (chosen such that they provide correct geometry deformations – too big values could lead to mesh destruction or deterioration of its quality – self intersecting elements with negative volumes are most common source of errors in case of too strong deformations),
  • maximum number of objective function evaluations (it is recommended to set this parameter to a moderately big value – the default situation is that the optimizer should converge and stop the iterative optimization procedure after a certain number of trials; however, if the optimizer fails to converge, which happens in more advanced technical problems with many constraints, the process will stop anyway after hitting maximum allowed number of function evaluations. Therefore, it is valuable to record all designs that have been tested, because we can then pick the one that proved to be the best so far.

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.