-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrae2822_sim.m
99 lines (78 loc) · 2.97 KB
/
rae2822_sim.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% RAE 2822 Fall 2024 %%%%
%%%% Code for Writing .geo File %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear all
%% Code to obtain y value for RAE 2822 airfoil trial case
% Freestream conditions/Inputs definition
z = 8.567e+3; % (m) Sea-level conditions
FTd = FTdclass(z);
M = 0.8;
FV = FVclass('M',M);
[props] = thermoprops(FTd,'alt',FV,'M');
L = 1; % (m) Airfoil chord length
var_y = 0; % Calculate y for a specified yplus value
y1 = 0;
yplus = 5;
[y1,yplus,cf,tau_w] = wallyp(props,M,L,var_y,y1,yplus);
%% Airfoil
% Domain extent multipliers [x-extent;y-extent]:
% C-type mesh
% y-extent: Radius of the inlet arc
% x-extent: Outlet standoff distance
domainext = 5*[L;L];
% Filename for airfoil co-ordinates file
filename = 'rae2822airfoil.txt';
% Function to get airfoil co-ordinates and find out airfoil attributes
[airfoildata]=airfoil_cond(filename,'Sharp',domainext);
%% Boundary Layer calculations
% Progression ratio
blr = 1/0.95;
% BL thickness
blthick = 0.01*L; % (m) 1% of chord length
% Element sizes
blmin = y1; % (m)
nbl = round(log(1-blthick*(1-blr)/y1)/log(blr));
blmax = blmin*blr^(nbl); % (m)
%% Meshing parameters
% Airfoil Curve lengths
le_length = airfoildata.le_length;
u_length = airfoildata.u_length;
cxt_length = sum(u_length(find(u_length(:,1)==airfoildata.cx_ind(1)):end,2));
l_length = airfoildata.l_length;
cxb_length = sum(l_length(find(l_length(:,1)==airfoildata.cx_ind(2)):end,2));
% te_length = airfoildata.te_length;
te_ind = airfoildata.te_ind;
% Domain curve lengths
i_arc_length = pi*domainext(2);
tev_length = domainext(2);
% Number of elements on curves
rle = 1;
nle = 120;
rv = 1/0.95; % Progression
nv = 140;
rairfoil = 0.1; % Bump
rcxh = 2; % Bump
nairfoil = 200;
rwake = 1/0.95; % Progression
rfarteo = 0.2; % Bump
nwake = 200;
nte = 5; % Finest elements in the mesh
meshdata = airfoilmeshdataclass(nle,rle,nv,rv,nairfoil,rairfoil,rcxh,nwake,...
rwake,rfarteo,nte,blr,blthick,blmin,blmax);
%% Y+ calculation
eymin_max = 5*L*(1-rv)/(1-rv^nv);
eymin_min = min(5*L-airfoildata.A(find(airfoildata.A(1:te_ind,2)==max(airfoildata.A(1:te_ind,2)),1),2),...
5*L-airfoildata.A(find(airfoildata.A(te_ind:end,2)==min(airfoildata.A(te_ind:end,2)),1),2))...
*(1-rv)/(1-rv^nv);
var_y = 1; % Calculate yplus for existing y value
y1 = eymin_max;
yplus_max = 0;
[y1,yplus_max,~,~] = wallyp(props,M,L,var_y,y1,yplus_max);
y2 = eymin_min;
yplus_min = 0;
[y2,yplus_min,~,~] = wallyp(props,M,L,var_y,y2,yplus_min);
%% Create .geo file for Gmsh
% filename = 'rae2822.geo';
% d2gfile_airfoil(airfoildata,meshdata,filename);