Delaunay Triangulation, is one of the most interesting methods in computational engineering, it is mainly used for domain Discretization in which the domain could be an area, volume, polygon etc… Examples of application of Delaunay triangulation are finite element methods and computational geometry.
Best to my knowledge, the term ‘Delaunay’ is originated from the famous French painter, Robert Delaunay and his wife, who used Orphism in their art work. What is Orphism? Why it is related to triangulation? Me and you don’t have to know
In this tutorial, I will explain how to generate Delaunay triangulation around NACA-0012 airfoil in MATLAB, thanks to CGAL (Computational Geometry Algorithms Library), the algorithms is done for us, so we just have to prepare our data and then execute few MATLAB commands to achieve our goal.
Here are the main steps to generate the Delaunay triangulation in MATLAB:
1- Draw a boundary Box.
2- Draw the selected geometry
3- Define the in/out region
4- Insert seed points
5- Connects these points to form the triangulation
Boundary Box
The boundary box is a simple rectangular shape that defines the boundaries of the domain, it must have enough room to accommodate the airfoil or any selected geometry, and so any point outside this domain is ignored. In fluid dynamics schemes, the boundary box is used to define the boundary conditions, to make it easier to understand, let us give a case of an airfoil inside a wind tunnel, the boundary box will represent the wind tunnel and we use it to define the boundary condition i.e wind flow direction and velocity, density, temperature etc…
Now, to form the dynamic Delaunay triangulation i.e. the density of the triangles in the domain area can be scaled, we have to define the boundary edges/lines with equally spaced points, in which the number of points can be controlled, for example, a horizontal line with length of 10 at y=0 can have the following coordinates:
X = [0,5,10] y = [0,0,0]
or
x= [0,2,4,6,8,10] y = [0,0,0,0,0,0]
or
x= [0,1,2,3,4,5,6,7,8,9,10] y=[0,0,0,0,0,0,0,0,0,0,0]
All these coordinates will form the same line although the points densities are different, for a dynamic structure, we must be able to control the number of these points.
So we will start defining the boundary box vertices as follow:
| Code: |
%========================
% GEOMETRY DEFINITION
%========================
left_edge = -2;
right_edge = 2;
bottom_edge = -2;
top_edge = 2;
sub_points = 10;
%Outer Lower Boundary
XLOWER = linspace(left_edge,right_edge,sub_points);
YLOWER = ones(size(XLOWER));
LOWER = [XLOWER; -YLOWER];
%Outer Upper Boundary
% XUPPER = linspace(2,-2,25);
XUPPER = linspace(right_edge,left_edge,33);
YUPPER = ones(size(XUPPER));
UPPER = [XUPPER; YUPPER];
%Left Wall Boundary
WLEFTY = linspace(top_edge,bottom_edge,sub_points);
WLEFTX = ones(size(WLEFTY)).*-2;
WALLLEFT = [WLEFTX; WLEFTY];
%Right Wall Boundary
% WRIGHTY = linspace(-1,1,25);
WRIGHTY = linspace(-0.916667,1,33);
WRIGHTX = ones(size(WRIGHTY)).*2;
WALLRIGHT = [WRIGHTX; WRIGHTY];
%OUTER BOUNDARY
OuterProfile = [LOWER, WALLRIGHT, UPPER, WALLLEFT]';
|
In this part of the code, the
*_edge defines the boundary box size, the [i]sub_point[/b] variables, is the number of N pointes in each edge/line, we change this variable to increase or decrease the resolution of the mesh.
Then we define the X and Y point using the
linspace function, in which we input the start point and the end point of the line, and the number of the intermediate points which are spaced evenly across the line. Take the example of
XLOWER which an array contains the X point of the lower edge of the bourndary box, the X coordinates starts from
left edge to the
right edge and have number of
sub_points spaced evenly, please don’t pay attention to my choice of variables name, they are just RANDOM!. For the Y coordinates of the lower edge, since the line is horizontal, i.e all the Y coordinates are same, we use the function
ones which construct an array of
1 with same size as
XLOWER, then multiply it by our desired Y coordinates value and store it into
YLOWER variable. An example
| Code: |
XLOWER = [-2, -1.56, -1.11, -0.67, -0.22, 0.22, 0.67, 1.11, 1.56, 2]
YLOWER = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
|
If we multiply YLOWER by 2, then all the array element will be 2 and so, the line will shift upwards.
Final step in this part is consolidate X and Y coordinates in into single array (
LOWER), so one array will hold the X,Y coordinates for a single edge, then re-consolidate all the edges in to one array (
OuterProfile) which will represent the boundary box
In the next step will go through adding the airfoil coordinates, and prepare our geometry profile for Delaunay triangulation! So stay tuned
