Updated air density-related functionalities in Linearization, Certification, AeroDyn, and WindTurbineDesign.m

Author: akpamososuryo

subfunctions/AeroDyn.m

@@ -1,4 +1,4 @@
-function AeroDyn(Blade,Airfoil,Tower,mode)
+function AeroDyn(Blade,Airfoil,Tower,mode,AirDensity)
 
 if contains(mode,'Linearize')
     AFAeroMod = 1;
@@ -22,7 +22,7 @@ function AeroDyn(Blade,Airfoil,Tower,mode)
 fprintf(fid, 'False           TwrAero            - Calculate tower aerodynamic loads? (flag)\n');
 fprintf(fid, 'False          FrozenWake         - Assume frozen wake during linearization? (flag) [used only when WakeMod=1 and when linearizing]\n');
 fprintf(fid, '======  Environmental Conditions  ===================================================================\n');
-fprintf(fid, '      1.225   AirDens            - Air density (kg/m^3)\n');
+fprintf(fid, '      %4.3f   AirDens            - Air density (kg/m^3)\n', AirDensity);
 fprintf(fid, '  1.464E-05   KinVisc            - Kinematic air viscosity (m^2/s)\n');
 fprintf(fid, '        335   SpdSound           - Speed of sound (m/s)\n');
 fprintf(fid, '======  Blade-Element/Momentum Theory Options  ====================================================== [used only when WakeMod=1]\n');

subfunctions/Certification.m

@@ -35,6 +35,7 @@ function Certification_OpeningFcn(hObject, eventdata, handles, varargin)
 handles.Drivetrain = varargin{5};
 handles.Control = varargin{6};
 handles.CertificationSettings = varargin{7};
+handles.AirDensity = varargin{8};
 
 % Disable simulation time and wind speed if a stepped wind is requested
 if handles.CertificationSettings.Wind.Type == 2
@@ -478,6 +479,7 @@ function Start_Callback(hObject, eventdata, handles)
 Drivetrain = handles.Drivetrain;
 Control = handles.Control;
 CertificationSettings = handles.CertificationSettings;
+AirDensity = handles.AirDensity;
 TMax = CertificationSettings.Run.Time;
 T = CertificationSettings.Wind.T;
 Ly = CertificationSettings.Wind.Ly;
@@ -519,11 +521,11 @@ function Start_Callback(hObject, eventdata, handles)
 % Steady state curves
 disp('Determining steady state rotational speeds and pitch angles...')
 if CertificationSettings.Wind.Type == 2
-    [~, ~, OmegaU, PitchAngle] = SteadyState(Blade, Airfoil, Drivetrain, Control, CertificationSettings.Wind.Step);
+    [~, ~, OmegaU, PitchAngle] = SteadyState(Blade, Airfoil, Drivetrain, Control, CertificationSettings.Wind.Step, AirDensity);
     OmegaU = OmegaU*ones(length(CertificationSettings.Run.WindSpeed));
     PitchAngle = PitchAngle*ones(length(CertificationSettings.Run.WindSpeed));
 else
-    [~, ~, OmegaU, PitchAngle] = SteadyState(Blade, Airfoil, Drivetrain, Control, CertificationSettings.Run.WindSpeed);
+    [~, ~, OmegaU, PitchAngle] = SteadyState(Blade, Airfoil, Drivetrain, Control, CertificationSettings.Run.WindSpeed, AirDensity);
 end
 RPM = OmegaU * 60/(2*pi);
 
@@ -535,7 +537,7 @@ function Start_Callback(hObject, eventdata, handles)
 % Turbine input files
 TMax = CertificationSettings.Run.Time;
 FASTinput(Control.DT, TMax);
-AeroDyn(Blade,Airfoil,Tower,string(CertificationSettings.Mode.Type));
+AeroDyn(Blade,Airfoil,Tower,string(CertificationSettings.Mode.Type), AirDensity);
 
 % Send to base workspace and make structures available for Simulink and run the simulation
 assignin('base', 'FAST_InputFileName', [pwd, filesep 'subfunctions' filesep 'inputfiles' filesep 'FAST.fst']);

subfunctions/Linearization.m

@@ -26,6 +26,7 @@ function Linearization_OpeningFcn(hObject, eventdata, handles, varargin)
 handles.Nacelle = varargin{4};
 handles.Control = varargin{5};
 handles.Drivetrain = varargin{6};
+handles.AirDensity = varargin{7};
 
 % Propose values for the wind speed range
 set(handles.WindSpeed_From, 'String', int2str(ceil(handles.Control.WindSpeed.Cutin)))
@@ -71,6 +72,7 @@ function SetFullLoad_Callback(hObject, eventdata, handles)
 Airfoil = handles.Airfoil;
 Drivetrain = handles.Drivetrain;
 Control = handles.Control;
+AirDensity = handles.AirDensity;
 
 % Find rated wind speed by comparing demanded torque at rotational speed C
 disp('Searching for rated wind speed...')
@@ -91,7 +93,7 @@ function SetFullLoad_Callback(hObject, eventdata, handles)
     for iter = 1:100
         U_new = U1 + abs((Control.Torque.Demanded-Qr1)/(Qr1-Qr2))*(U2-U1);
         [~, CQr] = PerformanceCoefficients(Blade, Airfoil, Control.Pitch.Fine, (OmegaC/Drivetrain.Gearbox.Ratio)*Blade.Radius(end)/U_new);
-        Qr_new = 0.5*CQr*1.225*U_new^2*pi*Blade.Radius(end)^3*Drivetrain.Gearbox.Efficiency/Drivetrain.Gearbox.Ratio;
+        Qr_new = 0.5*CQr*AirDensity*U_new^2*pi*Blade.Radius(end)^3*Drivetrain.Gearbox.Efficiency/Drivetrain.Gearbox.Ratio;
         if abs((Qr_new - Control.Torque.Demanded)/Control.Torque.Demanded) < 0.005
             success = true;
             break
@@ -167,6 +169,7 @@ function Linearize_Callback(hObject, eventdata, handles)
 Nacelle = handles.Nacelle;
 Control = handles.Control;
 Drivetrain = handles.Drivetrain;
+AirDensity = handles.AirDensity;
 
 % Run modal analysis for 0 rpm
 disp('Finding non-rotating mode shapes...')
@@ -199,7 +202,7 @@ function Linearize_Callback(hObject, eventdata, handles)
 
 % Steady state curves
 disp('Determining steady state rotational speeds and pitch angles...')
-[~, ~, OmegaU, PitchAngle] = SteadyState(Blade, Airfoil, Drivetrain, Control, WindSpeeds);
+[~, ~, OmegaU, PitchAngle] = SteadyState(Blade, Airfoil, Drivetrain, Control, WindSpeeds, AirDensity);
 RPM = OmegaU * 60/(2*pi);
 
 % Avoid some common errors with linearization
@@ -224,7 +227,7 @@ function Linearize_Callback(hObject, eventdata, handles)
 FAST_InputFileName = [pwd, filesep 'subfunctions' filesep 'inputfiles' filesep 'FAST.fst'];
 
 % Turbine input files
-AeroDyn(Blade,Airfoil,Tower,LinMode);
+AeroDyn(Blade,Airfoil,Tower,LinMode,AirDensity);
 ServoDyn(Drivetrain,Control,LinMode);
 
 % Preload the OutList
@@ -304,11 +307,10 @@ function Linearize_Callback(hObject, eventdata, handles)
     % Set input
     Blade = handles.Blade;
     Drivetrain = handles.Drivetrain;    
-
-    rhoAir = 1.225; % Density of air [kg/m3]
+    AirDensity = handles.AirDensity;
 
     % Calculate aerodynamic torque
-    Q = 0.5*CQr*rhoAir*windSpeed^2*pi*Blade.Radius(end)^3*Drivetrain.Gearbox.Efficiency/Drivetrain.Gearbox.Ratio;
+    Q = 0.5*CQr*AirDensity*windSpeed^2*pi*Blade.Radius(end)^3*Drivetrain.Gearbox.Efficiency/Drivetrain.Gearbox.Ratio;
 
 %% Wind speed steps
 function WindSpeed_From_Callback(hObject, eventdata, handles)

subfunctions/WindTurbineDesign.m

@@ -3025,7 +3025,8 @@ function Linearization_Callback(hObject, eventdata, handles)
     handles.Tower, ...
     handles.Nacelle, ...
     handles.Control, ...
-    handles.Drivetrain);
+    handles.Drivetrain,...
+    handles.AirDensity);
 
 % Enable window
 for i = 1:length(buttons)
@@ -3058,7 +3059,8 @@ function Certification_Callback(hObject, eventdata, handles)
     handles.Nacelle, ...
     handles.Drivetrain, ...
     handles.Control, ...
-    handles.CertificationSettings);
+    handles.CertificationSettings,...
+    handles.AirDensity);
 
 % Enable window
 for i = 1:length(buttons)