Tutorial 6 – Shape Optimization

Posted on July 17, 2017 - Tutorials

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

Introduction

This tutorial introduces geometry manipulation functions available in QuickerSim CFD Toolobx to allow for automatic shape deformation during optimization process. We have defined two distinct approaches to geometry parametrization. They differ in the following manner:

  1. Methodology 1 applies, if you want to deform boundaries of your domain to design optimal shape of a body/channel etc. In that case, you are about to place control points on the geometry to be able to define shape deformation which will follow displacement of the control points and will smoothly adapt original geometry to these changes.
  2. Methodology 2 is intended for use in cases, in which you deal with fluid domain, which contains some distinct bodies. Imagine that this applies when you for example want to optimize your design with respect to optimal position or orientation of these bodies inside of the domain. This means that you will be able to translate and rotate these bodies and computational mesh will be fitted, but shape of each of these bodies will remain the same. They will only be differently placed and/or oriented in space. Methodoly 2 is currently supported in 2-D only.

Methodology 1

Assume that we have a 3-D mesh loaded in MATLAB’s workspace. We wiil define two control points on geometry and illustrate meaning of the parameters.

The first methodology uses parametrization with radial basis functions. This means that each of the radial basis functions has a certain value in each of the nodes in the mesh. This values are used as weights, which weight how strongly a given node shall be displaced. The weight at the control point equals 1 and at other points the weight field is spanned by the radial basis function (RBF). A simple RBF in 1 dimension is defined by Gaussian distribution $f(\xi) = exp(-\frac{\xi^2}{2\sigma^2})$ In the above $\xi$ stands for distance of a given mesh node from the control point and parameter $\sigma$ defines extension range of the basis functions. From engineering point of view: the greater $\sigma$, the bigger radius of influence of the given RBF will be and the smoother the deformed geometry will be.

Below, we plot a 1-D radial basis function for $\sigma = 1$ and $\sigma = 0.25$.

Now, we present an example code which parametrizes and deforms geometry. Control points of subsequent radial basis functions are placed at $(0, 0.038, 0.03)$ and $(-0.02, 0.01 0.1)$. Parameters $\sigma$ are 0.05 and 0.03, respectively.

clc;
clear;

load('dieMesh85k.mat');
figure(1);
displayMesh(p,t);


% Build radial basis functions (first parameter in every row is sigma
% parameter, the next 2 in case of 2-D geometry and 3 in case of 3-D
% geometry parameters are coordinates of the control points).
rbf = [0.05 0.0 0.038 0.03; 0.03 -0.02 0.01 0.1]';


% Initialize radial basis functions and return output arguments. TM stores
% a transformation matrix and pinit coordinates of the original, undeformed
% geometry.
[TM, pinit] = initDeformMeshRBF(p,rbf);


% Define displacement of each of the control points - usually, this comes
% in vector of design variables from an optimization algorithm.
disp = [0 0.02 0; 0 -0.02 0];


% Call function, which actually does the deformation.
p = deformMeshRBF(pinit,TM,disp);


% Display the result
figure(2);
displayMesh(p,t);


% Below we show influence of the $\sigma$ parameter on geometry
% deformation. The only change in the above algorithm is change of the
% $\sigma$ parameter of the second RBF from $\sigma = 0.03$ to $\sigma =
% 0.01$.

CAD Export after Shape Optimization

An optimized geometry can be exported to a CAD environment in the form of a *.stl file. This can be done with the following function. Please, be cautious with extensive use of this function, because it is slow and will be optimized in future releases.

exportStl(p,e,'optimalDesign.stl');

% If you want to export to *.stl only part of the model (e.g. boundary with
% id = 2, you can extract this boundary from e array as follows:
ebound2 = e(:,e(7,:)==2);
exportStl(p,ebound2,'sideFace.stl');

Methodology 2

In this methodology we first need to define parts of the boundary which should be moving and which not. We import a mesh file, where inlet, outlet and channel walls are denoted with ids 10, 11, 12 and the airfoil which we want to place in different place and orientation has id = 13.

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


% The result in CAD software is visible below.


% We now need to init deformation function by specifying which boundaries
% will be stationary and which are about to be transformed.
[wdef, pinit] = initDeformMesh2D(p,e,t,13,[10 11 12]);


% In the next step we define initial position of the rotation axis of the
% object. Here we will assume that the axis is placed at x = 0.25, y = 0.
% Additionally we define, how strongly we want to displace the object. Let
% us displace it by dx = -0.2, dy = 0.1 and rotate clockwise by 0.6 rad. The
% information about initial axis position go into the first row of the
% following matrix and information about current position and rotation into
% the second row:
aP = [0.25 0 0; -0.2 0.1 -0.6];
p = deformMesh2D(pinit,wdef,aP);
% Display deformed mesh
displayMesh2D(p,t);


% Please note that very sharp displacements or rotations can lead to mesh
% destruction, which can form elements with negative volumes. An example of
% such an erroneous behavior is depicted in figure below.

 

 

Tags