FLOW SYSTEM ANALYSIS Part – 1: Unsteady Kernel function for 3D incompressible potential flow
Starting with the definition of acceleration potential, the equation for kernel function is derived and the steps are detailed as below:
1|Page
2|Page
3|Page
4|Page
5|Page
Part – 2: Tapered Wing – Analysis
Generate Grid points
g Input Variables n i s s e c o r p e r P
Doublet at 3 locations
Compute Kappa
Parabolic Approximation
Determine AIC
M L D
Compute I
M L D
Compute Lift, Moment, and its g n coefficents i s s e c o r P t s o P
The algorithm of the program sequence is briefly explained above. The grid points are generated based on the slope of the lines and similarity principle. The process to determine AIC involve the processes detailed in the referenced paper. All relevant sub-routines are given by different files to aid debugging process. The list of files used in the program is as follows: Project_Final. m –
Main File
getKappa.m, get I1.m, getI2.m, getI1pos.m, getI2pos.m – Subroutines For various cases of reduced frequencies, the lift, moment and coefficients are computed and the results are as follows:
For M = 0.0 Case No
Lift
Moment
CL
CM
-6.4807e+02 + 8.5323e+02i
-0.0547 + 0.0375i
-0.0720 + 0.0948i
1
--32.7946 +22.5023i
2
-34.3305 +33.6495i
-6.3530e+02 + 1.1559e+03i
-0.0572 + 0.0561i
-0.0706 + 0.1284i
3
-33.7213 +45.3042i
-5.6543e+02 + 1.4620e+03i
-0.0562 + 0.0755i
-0.0628 + 0.1624i
4
-30.9916 +56.9340i
-4.4079e+02 + 1.7575e+03i
-0.0517 + 0.0949i
-0.0490 + 0.1953i
6|Page
For M = 0.5 Case
Lift
Moment
CL
1
-38.2409 + 6.3813i
-7.8320e+02 + 4.5654e+02i
-0.0736 + 0.0123i
-0.1005 + 0.0586i
2
-41.4994 +15.5025i
-8.1313e+02 + 7.0952e+02i
-0.0799 + 0.0298i
-0.1043 + 0.0910i
3
-42.7592 +25.3623i
-7.8959e+02 + 9.7169e+02i
-0.0823 + 0.0488i
-0.1013 + 0.1247i
4
-42.0073 +35.4362i
-7.1395e+02 + 1.2293e+03i
-0.0808 + 0.0682i
-0.0916 + 0.1577i
No
CM
Code used for the analysis is attached as a part of the appendix. The complete results set have been added as an attachment to this file (Case-- 1/2 – a / b / c / d.mat). Plots generated are shown below:
Fig.1 Plot of Pressure & Lift at k=0.086, Mach = 0.0
7|Page
Fig.2 Plot of Pressure & Lift at k=0.010, Mach = 0.0
Fig.3 Plot of Pressure & Lift at k=0.114, Mach = 0.0
8|Page
Fig.4 Plot of Pressure & Lift at k=0.128, Mach = 0.0
Fig.5 Plot of Pressure & Lift at k=0.086, Mach = 0.5
9|Page
Fig.6 Plot of Pressure & Lift at k=0.100, Mach = 0.5
Fig.7 Plot of Pressure & Lift at k=0.114, Mach = 0.5
10 | P a g e
Fig.8 Plot of Pressure & Lift at k=0.128, Mach = 0.5 Inferences:
To account for compressibility, the Prandtl-Glauert coefficient has been multiplied to the coefficients of lift, pressure and moment.
As the Mach number is increased from 0.0 to 0.5 the lift & moment changes as its dependent upon the compressibility nature. The Coefficient of Lift & Moment reach their maximum value much faster as the Mach number is increased from 0.0 to 0.5.
11 | P a g e
References: 1.
DLM
Tools
from
University
of
Minnesota
-
http://gitpaaw.umnaem.webfactional.com/tools.git 2.
https://www.grc.nasa.gov/www/k-12/airplane/lifteq.html
3.
http://www.grc.nasa.gov/WWW/wind/valid/fplam/fplam01/fplam01.ht ml
4.
Effects of Mach number on the maximum lift coefficient of wings of NACA
– 2300 Series Air Foil Sections 5. 6.
Aeroservoelasticty by Ashish Tewari, ISBN 978-1-4939-2367-0 The Effect of Wing Damage on Aero elastic Behaviour by Howard J. Conyers Department of Mechanical Engineering and Materials Science Duke University.
7.
AIAC-2007-109,
AEROELASTIC
ANALYSIS
OF
A
TAPERED
AIRCRAFT WING S. Durmaz , O. Ozdemir Ozgumus , M. O. Kaya Istanbul Technical University, Faculty of Aeronautics and Astronautics, 34469, Maslak, Istanbul, Turkey 8.
Aero-elastic characteristics of tapered plate wings ARTICLE in FINITE ELEMENTS IN ANALYSIS AND DESIGN · FEBRUARY 2015 Impact Factor: 2.02 · DOI: 10.1016/j.finel.2014.09.009
12 | P a g e
APPENDIX
D: \ Nat i onal Uni ver s i t y Of Si ngapor e. . . \ PROJ ECT_ FI NAL. m Page 1 8 Apr i l , 2016 12: 09: 28 PM %% Flutter Analysis for Tapered Wing Geometry %% %% References : DLM Tools from University of Minnesota
- http://gitpaaw.umnaem.
webfactional.com/tools.git %% %% INITIALIZATION
r f = 0. 086: 0. 014: 0. 128; for num = 1: 4 r ho = 2. 38e- 3; M = 0. 5; k = r f ( num) ; sw = 40; ds = sw/ 8; sl ope_u = 0. 1;
% Reduced Frequencies for all cases % % Loop for Various Reduced Frequencies % % Density of Air % % Mach Number % % For individual cases % % Span Length % % Increment Spanwise % % Slope dy/dx = 4/40 using median principle
%
sl ope_l = - sl ope_u;
% Slope of symmetric opp. line
%% GRID COORDINATES
i = 0; j = 0; coor di nat e = zeros( 63, 4) ; panel = zer os( 48, 5) ; for t = 0: ds: sw; i = i +1; xl ( i , : ) = s l ope_ u* t ; % Lower Boundary % xu( i , : ) = sl ope_l *t + 19; %#ok<*SAGROW> % Upper Boundary % dx( i , : ) = ( xu( i , 1) - x l ( i , 1) ) / 6; % Increment % xl _ 1( i , : ) = xl ( i , 1) - dx( i , 1) ; xu_ 1( i , : ) = xu( i , 1) - dx( i , 1) ; for q = xl _ 1( i , 1) : dx( i , 1) : xu_ 1( i , 1) ; j = j +1; n( j , : ) = j ; y( j , : ) = t ; x_1( j , : ) = q + dx( i , 1) ; coor di nat e( j , 1) = n( j , 1) ; % Node Number % c oor di nat e( j , 2) = x_ 1( j , 1) ; % Chordwise Coordinate % c oor di nat e( j , 3) = y( j , 1) ; % Spanwise Coordinate % end end
h = 0; %% PANELS for
PB = 1: 1: 48; if( PB==7) h = 7; end if
( PB==13) h = 14;
end if(
PB==19) h = 21;
end if(
PB==25) h = 28;
end if(
PB==31) h = 35;
D: \ Nat i onal Uni ver s i t y Of Si ngapor e. . . \ PROJ ECT_ FI NAL. m Page 2 8 Apr i l , 2016 12: 09: 28 PM end if(
PB==37) h = 42;
end if(
PB==43) h = 49;
end if
( PB~=7| | PB~=13| | PB~=19| | PB~=25| | PB~=31| | PB~=37| | PB~=43) h = h+1;
end
panel ( PB, 1) panel ( PB, 2) panel ( PB, 3) panel ( PB, 4) panel ( PB, 5)
= = = = =
PB; h; h + 1; h + 7; h + 8;
% Panel Number % % Coordinate 1 % % Coordinate 2 % % Coordinate 3 % % Coordinate 4 %
end
%% CELL CENTRES OF PANELS
i _1 = 0; j _1 = 0; cor di nat e_c = zer os( 48, 4) ; for t 1 = ds/ 2: ds: sw- ds/ 2; i _1 = i _1 + 1; xl _ c( i _ 1, : ) = s l ope_ u* t 1; % Upper Boundary % xu_c( i _1, : ) = sl ope_l *t 1 + 19; % Lower Boundary % dx_ c( i _ 1, : ) = ( xu_ c( i _ 1, 1) - xl _ c( i _ 1, 1) ) / 6; % Increment % xl 1_ c( i _ 1, : ) = xl _ c( i _ 1, 1) - dx_ c( i _ 1, 1) + dx_ c( i _ 1, 1) / 2; xu1_ c( i _ 1, : ) = xu_ c( i _ 1, 1) - dx_ c( i _ 1, 1) - dx_ c( i _ 1, 1) / 2; for q_c = xl 1_c(i _1, 1) : dx_c(i _1, 1) : xu1_c(i _1, 1) ; j _1 = j _1+1; n_ c( j _ 1, : ) =j _ 1; y_ c( j _ 1, : ) = t 1; x_1_c( j _1, : ) = q_c + dx_c( i _1, 1) / 2; cor di nat e_c(j _1, 1) = n_c(j _1, 1) ; % Node Number % cor di nat e_c(j _1, 2) = x_1_c(j _1, 1) ; % Chordwise Coordinate % cor di nat e_c(j _1, 3) = y_c(j _1, 1) ; % Spanwise Coordinate % di f f _ c( j _ 1, 1) = 6* dx_ c( i _ 1, 1) ; end end
%% DOWNWASH POSITION
PO = zer os( PB, 3) ; N = si ze( PO, 1) ; PO( : , 1) = cor di nat e_c(: , 2) +( 3/ 8) *( coor di nat e( panel ( : , 3) , 2) - coor di nat e( panel ( : , 2) , 2) + ... c oor di nat e( panel ( : , 5) , 2) - c oor di nat e( panel ( : , 4) , 2) ) ; PO( : , 2) = coor di nat e( panel ( : , 2) , 3) + ( coor di nat e( panel ( : , 4) , 3) - coor di nat e( panel ( : , 2) , 3) ) / 2; PO( : , 3) = 0; %% DOUBLET - START OF WING
PS = zer os( PB, 3) ; PS( : , 1) = coor di nat e( panel ( : , 2) , 2) +( 1/ 4) *( coor di nat e( panel ( : , 3) , 2) - coor di nat e( panel ( : , 2) , 2) ) ; PS( : , 2) = coor di nat e( panel ( : , 2) , 3) ; PS( : , 3) = 0;
D: \ Nat i onal Uni ver s i t y Of Si ngapor e. . . \ PROJ ECT_ FI NAL. m Page 3 8 Apr i l , 2016 12: 09: 28 PM %% DOUBLET - TIP OF WING
PE = zer os( PB, 3) ; PE( : , 1) = coor di nat e( panel ( : , 4) , 2) + ( ( coor di nat e( panel ( : , 5) , 2) c oor di nat e( panel ( : , 4) , 2) ) / 4) ; PE( : , 2) = coor di nat e( panel ( : , 4) , 3) ; PE( : , 3) = 0; %%
...
DOUBLET - HALF SPAN OF WING
PM = zer os( PB, 3) ; PM( : , 1) = ( PS( : , 1) + PE( : , 1) ) / 2; PM( : , 2) = ( PS( : , 2) + PE( : , 2) ) / 2; PM( : , 3) = ( PS( : , 3) + PE( : , 3) ) / 2; %% DEFINING SPAN AND CHORD OF EACH PANEL
Z_1 = zer os( 48, 1) ; SP = ( 0. 5*( coor di nat e( panel ( : , 4) , 3) - coor di nat e( panel ( : , 2) , 3) ) ) ' ; CP = ( ( abs( ( coor di nat e( panel ( : , 2) , 2) - coor di nat e( panel ( : , 3) , 2) ) ) + abs ( ( c oor di nat e( panel ( : , 4) , 2) - c oor di nat e( panel ( : , 5) , 2) ) ) ) / 2) ' ;
...
%% DLM NORMALWASH MATRIX - UNNSTEADY CASE %% ROOT DOUBLET LOCATION
x01 = bsxf un( @mi nus, PO( : , 1) , PS( : , 1) ' ) ; y01 = bsxf un( @mi nus, PO( : , 2) , PS( : , 2) ' ) ; z01 = r epmat ( Z_1' , 48, 1) ; %% SEMI-SPAN DOUBLET LOCATION
x02 = bsxf un( @mi nus , PO( : , 1) , PM( : , 1) ' ) ; y02 = bsxf un( @mi nus , PO( : , 2) , PM( : , 2) ' ) ; z02 = r epmat ( Z_1' , 48, 1) ; %% TIP OF DOUBLET LOCATION
x03 = bsxf un( @mi nus, PO( : , 1) , PE( : , 1) ' ) ; y03 = bsxf un( @mi nus, PO( : , 2) , PE( : , 2) ' ) ; z03 = r epmat ( Z_1' , 48, 1) ; cosGamma = ones( 48, 48) ; si nGamma = zer os( 48, 48) ; %% KAPPA
Ki _w = get Kappa( x01, y01, z01, cosGamma, si nGamma, k, M) ; Ki _0 = get Kappa( x01, y01, z01, cosGamma, si nGamma, 0, M) ; Ki = Ki _w - Ki _0; Km_w = get Kappa( x02, y02, z02, cos Gamma, si nGamma, k, M) ; Km_0 = get Kappa( x02, y02, z02, cos Gamma, si nGamma, 0, M) ; Km = Km_w - Km_0; K0_w = get Kappa( x03, y03, z03, cos Gamma, si nGamma, k, M) ; K0_0 = get Kappa( x03, y03, z03, cos Gamma, si nGamma, 0, M) ; K0 = K0_w - K0_0; %% PARABOLIC APPROXIMATION
e1 = ( r epmat ( SP, N, 1) ) ; A = ( Ki - 2*Km+K0) . / ( 2*e1. ^2) ; B = ( K0- Ki ) . / ( 2*e1) ; C = Km; n0 = y02. *cos Gamma; zet a0 = z02. *cosGamma; r 2 = sqr t ( ( n0. ^2) + ( zet a0. ^2) ) ; %% NORMALWASH MATRIX FACTOR "I"
I = ( A. *( 2*e1) ) + ( ( 0. 5*B+n0. *A) . *l og( ( r 2. ^2 - 2*n0. *e1 + e1. ^2) . / ( r 2. ^2 + 2*n0. *e1 + e1. ^2) ) ) ... +( ( ( n0. ^2 - zet a0. ^2) . *A+n0. *B + C) . / abs( zet a0) . *atan( 2*e1. *abs( zet a0) . / ( r 2. ^2 -
D: \ Nat i onal Uni ver s i t y Of Si ngapor e. . . \ PROJ ECT_ FI NAL. m Page 4 8 Apr i l , 2016 12: 09: 28 PM e1. ^2) ) ) ; i nd = f i nd( zet a0==0) ; I 2 = ( ( A. *( 2*e1) ) + ( ( 0. 5*B+n0. *A) . *l og( ( ( n0- e1) . / ( n0+e1) ) . ^2) ) +( ( n0. ^2) . *A+n0. *B + C) . *( ( 2*e1) . / ( n0. ^2 - e1. ^2) ) ) ; I ( i nd) = I 2( i nd) ; D = r epmat ( CP, N, 1) . *I / ( pi *8) ;
...
%% VLM NORMALWASH MATRIX - QUASI - STEADY STATE
bet a = sqr t ( 1- ( M^2) ) ; %% AIC MATRIX
AI C = i nv( D) ; wbar =ones( 48, 1) ; %% PRESSURE DISTRIBUTION
P = AI C\ wbar ; P_R = r eal ( P) ; P_I = i mag( P) ; P_ABS = abs( P) ; P_PR = vec2mat ( P_R, 8) ; P_PI = vec2mat ( P_I , 8) ; f i gur e( num) subpl ot ( 2, 2, 1) s cat t er ( y, x_ 1) ; xl abel ( ' Span' ) ; yl abel ( ' Chor d' ) ; t i t l e( ' Gr i d poi nt s ' ) ; subpl ot ( 2, 2, 2) meshz( P_PR) ; xl abel ( ' Span Di vi si ons' ) ; yl abel ( ' Chor d Di vi s i ons ' ) ; z l abel ( ' P( R) ' ) ; t i t l e( ' Pr es sur e - Real Par t ' ) ; subpl ot ( 2, 2, 3) meshz( P_PI ) ; xl abel ( ' Span Di vi si ons' ) ; yl abel ( ' Chor d Di vi s i ons ' ) ; zl abel ( ' P( I ) ' ) ; t i t l e( ' Pr es sur e - I MG Par t ' ) ; Li f t = 0; MO = 0; for i = 1: 48 Li ( i , : ) = P( i ) * 2* SP( i ) * CP( i ) ; Li f t = Li f t + Li ( i , : ) ; MO = MO + Li f t *0. 4*CP( i ) ;
% delP / Q %
% Non dimensional Lift - Lift / Q %
end
L_D = vec2mat ( abs( Li ) , 8) ; subpl ot ( 2, 2, 4) meshz( L_ D) ; xl abel ( ' Span Di vi si ons' ) ; yl abel ( ' Chor d Di vi s i ons ' ) ; z l abel ( ' abs ( L) ' ) ; t i t l e( ' Li f t - Absol ut e' ) ; S = 0. 5*( ( coor di nat e( 7, 2) - coor di nat e( 1, 2) ) + ( coor di nat e( 63, 2) - coor di nat e( 57, 2) ) ) *sw; % Surface Area %
cavg = 0. 5*( ( coor di nat e( 7, 2) - coor di nat e( 1, 2) ) + ( coor di nat e( 63, 2) - coor di nat e
D: \ Nat i onal Uni ver s i t y Of Si ngapor e. . . \ PROJ ECT_ FI NAL. m Page 5 8 Apr i l , 2016 12: 09: 28 PM ( 57, 2) ) ) ; % Avg. Chord Length % CL = Li f t / ( bet a*S) ; CM = MO/ ( bet a*S*cavg) ; di sp( M) ; di s p( k) ; di spl ay( num) ; di s pl ay( Li f t ) ; di spl ay( MO) ; di spl ay( CL) ; di spl ay( CM) ; end
% Non dimensional Coefficent of Lift % % Non dimensional Coefficient of Moment %
