[1820] | 1 | !> @file wind_turbine_model_mod.f90 |
---|
[1912] | 2 | !------------------------------------------------------------------------------! |
---|
[1819] | 3 | ! This file is part of PALM. |
---|
| 4 | ! |
---|
[2000] | 5 | ! PALM is free software: you can redistribute it and/or modify it under the |
---|
| 6 | ! terms of the GNU General Public License as published by the Free Software |
---|
| 7 | ! Foundation, either version 3 of the License, or (at your option) any later |
---|
| 8 | ! version. |
---|
[1819] | 9 | ! |
---|
| 10 | ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY |
---|
| 11 | ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR |
---|
| 12 | ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. |
---|
| 13 | ! |
---|
| 14 | ! You should have received a copy of the GNU General Public License along with |
---|
| 15 | ! PALM. If not, see <http://www.gnu.org/licenses/>. |
---|
| 16 | ! |
---|
[2257] | 17 | ! Copyright 2009-2017 Carl von Ossietzky Universitaet Oldenburg |
---|
[2101] | 18 | ! Copyright 1997-2017 Leibniz Universitaet Hannover |
---|
[1912] | 19 | !------------------------------------------------------------------------------! |
---|
[1819] | 20 | ! |
---|
| 21 | ! Current revisions: |
---|
| 22 | ! ----------------- |
---|
[1913] | 23 | ! |
---|
[2233] | 24 | ! |
---|
[1913] | 25 | ! Former revisions: |
---|
| 26 | ! ----------------- |
---|
| 27 | ! $Id: wind_turbine_model_mod.f90 2563 2017-10-19 15:36:10Z Giersch $ |
---|
[2563] | 28 | ! Restart runs with wind turbine model are possible now. For this purpose, two |
---|
| 29 | ! new subroutines wtm_write_restart_data and wtm_read_restart_data had to be |
---|
| 30 | ! defined |
---|
| 31 | ! |
---|
| 32 | ! 2553 2017-10-18 08:03:45Z Giersch |
---|
[2553] | 33 | ! Bugfix of vertical loop in wtm_tendencies to account for different turbine |
---|
| 34 | ! heights, bugfix of the interpolation of the u-component concerning the |
---|
| 35 | ! vertical index and further small adjustments of the programming style |
---|
| 36 | ! |
---|
| 37 | ! 2410 2017-09-06 08:16:30Z Giersch |
---|
[2410] | 38 | ! Revise error message PA0462 |
---|
| 39 | ! |
---|
| 40 | ! 2349 2017-08-10 15:44:04Z Giersch |
---|
[2349] | 41 | ! Add parameter pitch_rate to namelist and revise/add error messages |
---|
| 42 | ! |
---|
| 43 | ! 2343 2017-08-08 11:28:43Z Giersch |
---|
[2341] | 44 | ! Unit correction in Doxygen comments |
---|
[2323] | 45 | ! |
---|
[2343] | 46 | ! 2323 2017-07-26 12:57:38Z Giersch |
---|
| 47 | ! Change unit number of file WTM_DATA from 201 to 90 |
---|
| 48 | ! |
---|
[2323] | 49 | ! 2322 2017-07-26 08:30:28Z Giersch |
---|
[2322] | 50 | ! Bugfix of error message and assign error numbers |
---|
| 51 | ! |
---|
| 52 | ! 2257 2017-06-07 14:07:05Z witha |
---|
[2257] | 53 | ! Bugfix: turb_cl_tab and turb_cd_tab were set to zero before being allocated |
---|
[1844] | 54 | ! |
---|
[2257] | 55 | ! 2233 2017-05-30 18:08:54Z suehring |
---|
| 56 | ! |
---|
[2233] | 57 | ! 2232 2017-05-30 17:47:52Z suehring |
---|
| 58 | ! Adjustments to new topography concept |
---|
| 59 | ! |
---|
[2153] | 60 | ! 2152 2017-02-17 13:27:24Z lvollmer |
---|
| 61 | ! Bugfix in subroutine wtm_read_blade_tables |
---|
| 62 | ! Addition of a tip loss model |
---|
| 63 | ! |
---|
[2016] | 64 | ! 2015 2016-09-28 08:45:18Z lvollmer |
---|
| 65 | ! Bugfix of pitch control |
---|
| 66 | ! |
---|
[2001] | 67 | ! 2000 2016-08-20 18:09:15Z knoop |
---|
| 68 | ! Forced header and separation lines into 80 columns |
---|
| 69 | ! |
---|
[1930] | 70 | ! 1929 2016-06-09 16:25:25Z suehring |
---|
| 71 | ! Bugfix: added preprocessor directives for parallel and serial mode |
---|
| 72 | ! |
---|
[1917] | 73 | ! 1914 2016-05-26 14:44:07Z witha |
---|
[1821] | 74 | ! Initial revision |
---|
| 75 | ! |
---|
[1819] | 76 | ! |
---|
| 77 | ! Description: |
---|
| 78 | ! ------------ |
---|
| 79 | !> This module calculates the effect of wind turbines on the flow fields. The |
---|
| 80 | !> initial version contains only the advanced actuator disk with rotation method |
---|
| 81 | !> (ADM-R). |
---|
| 82 | !> The wind turbines include the tower effect, can be yawed and tilted. |
---|
| 83 | !> The wind turbine model includes controllers for rotational speed, pitch and |
---|
| 84 | !> yaw. |
---|
| 85 | !> Currently some specifications of the NREL 5 MW reference turbine |
---|
| 86 | !> are hardcoded whereas most input data comes from separate files (currently |
---|
| 87 | !> external, planned to be included as namelist which will be read in |
---|
| 88 | !> automatically). |
---|
| 89 | !> |
---|
| 90 | !> @todo Revise code according to PALM Coding Standard |
---|
| 91 | !> @todo Implement ADM and ALM turbine models |
---|
| 92 | !> @todo Generate header information |
---|
[1917] | 93 | !> @todo Implement further parameter checks and error messages |
---|
[1819] | 94 | !> @todo Revise and add code documentation |
---|
| 95 | !> @todo Output turbine parameters as timeseries |
---|
| 96 | !> @todo Include additional output variables |
---|
[1864] | 97 | !> @todo Revise smearing the forces for turbines in yaw |
---|
| 98 | !> @todo Revise nacelle and tower parameterization |
---|
| 99 | !> @todo Allow different turbine types in one simulation |
---|
[1819] | 100 | ! |
---|
| 101 | !------------------------------------------------------------------------------! |
---|
| 102 | MODULE wind_turbine_model_mod |
---|
| 103 | |
---|
| 104 | USE arrays_3d, & |
---|
[2553] | 105 | ONLY: tend, u, v, w, zu, zw |
---|
[1819] | 106 | |
---|
[2553] | 107 | USE constants, & |
---|
| 108 | ONLY: pi |
---|
[1819] | 109 | |
---|
| 110 | USE control_parameters, & |
---|
[2563] | 111 | ONLY: dt_3d, dz, message_string, simulated_time, wind_turbine, & |
---|
| 112 | initializing_actions |
---|
[1819] | 113 | |
---|
| 114 | USE cpulog, & |
---|
| 115 | ONLY: cpu_log, log_point_s |
---|
| 116 | |
---|
| 117 | USE grid_variables, & |
---|
| 118 | ONLY: ddx, dx, ddy, dy |
---|
| 119 | |
---|
| 120 | USE indices, & |
---|
| 121 | ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, & |
---|
[2232] | 122 | nzb, nzt, wall_flags_0 |
---|
[1819] | 123 | |
---|
| 124 | USE kinds |
---|
| 125 | |
---|
| 126 | USE pegrid |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | IMPLICIT NONE |
---|
| 130 | |
---|
[1864] | 131 | PRIVATE |
---|
[1819] | 132 | |
---|
| 133 | ! |
---|
| 134 | !-- Variables specified in the namelist wind_turbine_par |
---|
| 135 | |
---|
[1864] | 136 | INTEGER(iwp) :: nairfoils = 8 !< number of airfoils of the used turbine model (for ADM-R and ALM) |
---|
[1839] | 137 | INTEGER(iwp) :: nturbines = 1 !< number of turbines |
---|
[1819] | 138 | |
---|
[1839] | 139 | LOGICAL :: pitch_control = .FALSE. !< switch for use of pitch controller |
---|
| 140 | LOGICAL :: speed_control = .FALSE. !< switch for use of speed controller |
---|
| 141 | LOGICAL :: yaw_control = .FALSE. !< switch for use of yaw controller |
---|
[2152] | 142 | LOGICAL :: tl_cor = .FALSE. !< switch for use of tip loss correct. |
---|
[1819] | 143 | |
---|
[1839] | 144 | REAL(wp) :: segment_length = 1.0_wp !< length of the segments, the rotor area is divided into |
---|
| 145 | !< (in tangential direction, as factor of MIN(dx,dy,dz)) |
---|
| 146 | REAL(wp) :: segment_width = 0.5_wp !< width of the segments, the rotor area is divided into |
---|
| 147 | !< (in radial direction, as factor of MIN(dx,dy,dz)) |
---|
| 148 | REAL(wp) :: time_turbine_on = 0.0_wp !< time at which turbines are started |
---|
| 149 | REAL(wp) :: tilt = 0.0_wp !< vertical tilt of the rotor [degree] ( positive = backwards ) |
---|
[1819] | 150 | |
---|
[1912] | 151 | REAL(wp), DIMENSION(1:100) :: dtow = 0.0_wp !< tower diameter [m] |
---|
[2553] | 152 | REAL(wp), DIMENSION(1:100) :: omega_rot = 0.9_wp !< inital or constant rotor speed [rad/s] |
---|
[1912] | 153 | REAL(wp), DIMENSION(1:100) :: phi_yaw = 0.0_wp !< yaw angle [degree] ( clockwise, 0 = facing west ) |
---|
| 154 | REAL(wp), DIMENSION(1:100) :: pitch_add = 0.0_wp !< constant pitch angle |
---|
| 155 | REAL(wp), DIMENSION(1:100) :: rcx = 9999999.9_wp !< position of hub in x-direction |
---|
| 156 | REAL(wp), DIMENSION(1:100) :: rcy = 9999999.9_wp !< position of hub in y-direction |
---|
| 157 | REAL(wp), DIMENSION(1:100) :: rcz = 9999999.9_wp !< position of hub in z-direction |
---|
| 158 | REAL(wp), DIMENSION(1:100) :: rnac = 0.0_wp !< nacelle diameter [m] |
---|
| 159 | REAL(wp), DIMENSION(1:100) :: rr = 63.0_wp !< rotor radius [m] |
---|
| 160 | REAL(wp), DIMENSION(1:100) :: turb_cd_nacelle = 0.85_wp !< drag coefficient for nacelle |
---|
| 161 | REAL(wp), DIMENSION(1:100) :: turb_cd_tower = 1.2_wp !< drag coefficient for tower |
---|
[1839] | 162 | |
---|
[1819] | 163 | ! |
---|
| 164 | !-- Variables specified in the namelist for speed controller |
---|
| 165 | !-- Default values are from the NREL 5MW research turbine (Jonkman, 2008) |
---|
| 166 | |
---|
[1912] | 167 | REAL(wp) :: rated_power = 5296610.0_wp !< rated turbine power [W] |
---|
[1839] | 168 | REAL(wp) :: gear_ratio = 97.0_wp !< Gear ratio from rotor to generator |
---|
[2341] | 169 | REAL(wp) :: inertia_rot = 34784179.0_wp !< Inertia of the rotor [kg*m2] |
---|
| 170 | REAL(wp) :: inertia_gen = 534.116_wp !< Inertia of the generator [kg*m2] |
---|
[1839] | 171 | REAL(wp) :: gen_eff = 0.944_wp !< Electric efficiency of the generator |
---|
| 172 | REAL(wp) :: gear_eff = 1.0_wp !< Loss between rotor and generator |
---|
| 173 | REAL(wp) :: air_dens = 1.225_wp !< Air density to convert to W [kg/m3] |
---|
| 174 | REAL(wp) :: rated_genspeed = 121.6805_wp !< Rated generator speed [rad/s] |
---|
[1912] | 175 | REAL(wp) :: max_torque_gen = 47402.91_wp !< Maximum of the generator torque [Nm] |
---|
[1839] | 176 | REAL(wp) :: slope2 = 2.332287_wp !< Slope constant for region 2 |
---|
| 177 | REAL(wp) :: min_reg2 = 91.21091_wp !< Lower generator speed boundary of region 2 [rad/s] |
---|
| 178 | REAL(wp) :: min_reg15 = 70.16224_wp !< Lower generator speed boundary of region 1.5 [rad/s] |
---|
[1912] | 179 | REAL(wp) :: max_trq_rate = 15000.0_wp !< Max generator torque increase [Nm/s] |
---|
| 180 | REAL(wp) :: pitch_rate = 8.0_wp !< Max pitch rate [degree/s] |
---|
[1839] | 181 | |
---|
[1912] | 182 | |
---|
[1819] | 183 | ! |
---|
| 184 | !-- Variables specified in the namelist for yaw control |
---|
| 185 | |
---|
[1839] | 186 | REAL(wp) :: yaw_speed = 0.005236_wp !< speed of the yaw actuator [rad/s] |
---|
| 187 | REAL(wp) :: max_miss = 0.08726_wp !< maximum tolerated yaw missalignment [rad] |
---|
| 188 | REAL(wp) :: min_miss = 0.008726_wp !< minimum yaw missalignment for which the actuator stops [rad] |
---|
| 189 | |
---|
[1819] | 190 | ! |
---|
| 191 | !-- Set flag for output files TURBINE_PARAMETERS |
---|
| 192 | TYPE file_status |
---|
| 193 | LOGICAL :: opened, opened_before |
---|
| 194 | END TYPE file_status |
---|
| 195 | |
---|
[1912] | 196 | TYPE(file_status), DIMENSION(500) :: openfile_turb_mod = & |
---|
| 197 | file_status(.FALSE.,.FALSE.) |
---|
[1819] | 198 | |
---|
| 199 | ! |
---|
| 200 | !-- Variables for initialization of the turbine model |
---|
| 201 | |
---|
[1839] | 202 | INTEGER(iwp) :: inot !< turbine loop index (turbine id) |
---|
| 203 | INTEGER(iwp) :: nsegs_max !< maximum number of segments (all turbines, required for allocation of arrays) |
---|
| 204 | INTEGER(iwp) :: nrings_max !< maximum number of rings (all turbines, required for allocation of arrays) |
---|
| 205 | INTEGER(iwp) :: ring !< ring loop index (ring number) |
---|
| 206 | INTEGER(iwp) :: rr_int !< |
---|
| 207 | INTEGER(iwp) :: upper_end !< |
---|
[1819] | 208 | |
---|
[1839] | 209 | INTEGER(iwp), DIMENSION(1) :: lct !< |
---|
[1819] | 210 | |
---|
[1912] | 211 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_hub !< index belonging to x-position of the turbine |
---|
| 212 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_smear !< index defining the area for the smearing of the forces (x-direction) |
---|
| 213 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_hub !< index belonging to y-position of the turbine |
---|
| 214 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_smear !< index defining the area for the smearing of the forces (y-direction) |
---|
| 215 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_hub !< index belonging to hub height |
---|
| 216 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_smear !< index defining the area for the smearing of the forces (z-direction) |
---|
| 217 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nrings !< number of rings per turbine |
---|
| 218 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nsegs_total !< total number of segments per turbine |
---|
[1819] | 219 | |
---|
[1912] | 220 | INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nsegs !< number of segments per ring and turbine |
---|
[1819] | 221 | |
---|
[1912] | 222 | ! |
---|
| 223 | !- parameters for the smearing from the rotor to the cartesian grid |
---|
[1864] | 224 | REAL(wp) :: pol_a !< parameter for the polynomial smearing fct |
---|
| 225 | REAL(wp) :: pol_b !< parameter for the polynomial smearing fct |
---|
[1912] | 226 | REAL(wp) :: delta_t_factor !< |
---|
[1864] | 227 | REAL(wp) :: eps_factor !< |
---|
[1839] | 228 | REAL(wp) :: eps_min !< |
---|
| 229 | REAL(wp) :: eps_min2 !< |
---|
| 230 | REAL(wp) :: sqrt_arg !< |
---|
[1819] | 231 | |
---|
[1839] | 232 | ! |
---|
| 233 | !-- Variables for the calculation of lift and drag coefficients |
---|
[1912] | 234 | REAL(wp), DIMENSION(:), ALLOCATABLE :: ard !< |
---|
| 235 | REAL(wp), DIMENSION(:), ALLOCATABLE :: crd !< |
---|
| 236 | REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_r !< radial segment length |
---|
| 237 | REAL(wp), DIMENSION(:), ALLOCATABLE :: lrd !< |
---|
[1864] | 238 | |
---|
[1912] | 239 | REAL(wp) :: accu_cl_cd_tab = 0.1_wp !< Accuracy of the interpolation of |
---|
| 240 | !< the lift and drag coeff [deg] |
---|
[1819] | 241 | |
---|
[1912] | 242 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cd_tab !< table of the blade drag coefficient |
---|
| 243 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cl_tab !< table of the blade lift coefficient |
---|
[1819] | 244 | |
---|
[1912] | 245 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: nac_cd_surf !< 3d field of the tower drag coefficient |
---|
| 246 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tow_cd_surf !< 3d field of the nacelle drag coefficient |
---|
[1819] | 247 | |
---|
| 248 | ! |
---|
| 249 | !-- Variables for the calculation of the forces |
---|
[1912] | 250 | |
---|
[1839] | 251 | REAL(wp) :: cur_r !< |
---|
| 252 | REAL(wp) :: phi_rotor !< |
---|
| 253 | REAL(wp) :: pre_factor !< |
---|
| 254 | REAL(wp) :: torque_seg !< |
---|
| 255 | REAL(wp) :: u_int_l !< |
---|
| 256 | REAL(wp) :: u_int_u !< |
---|
| 257 | REAL(wp) :: u_rot !< |
---|
| 258 | REAL(wp) :: v_int_l !< |
---|
| 259 | REAL(wp) :: v_int_u !< |
---|
| 260 | REAL(wp) :: w_int_l !< |
---|
| 261 | REAL(wp) :: w_int_u !< |
---|
[1912] | 262 | ! |
---|
| 263 | !- Tendencies from the nacelle and tower thrust |
---|
| 264 | REAL(wp) :: tend_nac_x = 0.0_wp !< |
---|
| 265 | REAL(wp) :: tend_tow_x = 0.0_wp !< |
---|
| 266 | REAL(wp) :: tend_nac_y = 0.0_wp !< |
---|
| 267 | REAL(wp) :: tend_tow_y = 0.0_wp !< |
---|
[1819] | 268 | |
---|
[1912] | 269 | REAL(wp), DIMENSION(:), ALLOCATABLE :: alpha_attack !< |
---|
| 270 | REAL(wp), DIMENSION(:), ALLOCATABLE :: chord !< |
---|
| 271 | REAL(wp), DIMENSION(:), ALLOCATABLE :: phi_rel !< |
---|
| 272 | REAL(wp), DIMENSION(:), ALLOCATABLE :: torque_total !< |
---|
| 273 | REAL(wp), DIMENSION(:), ALLOCATABLE :: thrust_rotor !< |
---|
| 274 | REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cl !< |
---|
| 275 | REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cd !< |
---|
| 276 | REAL(wp), DIMENSION(:), ALLOCATABLE :: vrel !< |
---|
| 277 | REAL(wp), DIMENSION(:), ALLOCATABLE :: vtheta !< |
---|
| 278 | |
---|
| 279 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rbx, rby, rbz !< coordinates of the blade elements |
---|
| 280 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rotx, roty, rotz !< normal vectors to the rotor coordinates |
---|
| 281 | |
---|
| 282 | ! |
---|
| 283 | !- Fields for the interpolation of velocities on the rotor grid |
---|
| 284 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_int !< |
---|
| 285 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_int_1_l !< |
---|
| 286 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_int !< |
---|
| 287 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_int_1_l !< |
---|
| 288 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_int !< |
---|
| 289 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_int_1_l !< |
---|
| 290 | |
---|
| 291 | ! |
---|
| 292 | !- rotor tendencies on the segments |
---|
| 293 | REAL(wp), DIMENSION(:), ALLOCATABLE :: thrust_seg !< |
---|
| 294 | REAL(wp), DIMENSION(:), ALLOCATABLE :: torque_seg_y !< |
---|
| 295 | REAL(wp), DIMENSION(:), ALLOCATABLE :: torque_seg_z !< |
---|
| 296 | |
---|
| 297 | ! |
---|
| 298 | !- rotor tendencies on the rings |
---|
| 299 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: thrust_ring !< |
---|
| 300 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: torque_ring_y !< |
---|
| 301 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: torque_ring_z !< |
---|
| 302 | |
---|
| 303 | ! |
---|
| 304 | !- rotor tendencies on rotor grids for all turbines |
---|
| 305 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: thrust !< |
---|
| 306 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: torque_y !< |
---|
| 307 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: torque_z !< |
---|
| 308 | |
---|
| 309 | ! |
---|
| 310 | !- rotor tendencies on coordinate grid |
---|
| 311 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rot_tend_x !< |
---|
| 312 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rot_tend_y !< |
---|
| 313 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rot_tend_z !< |
---|
| 314 | ! |
---|
| 315 | !- variables for the rotation of the rotor coordinates |
---|
| 316 | REAL(wp), DIMENSION(1:100,1:3,1:3) :: rot_coord_trans !< matrix for rotation of rotor coordinates |
---|
| 317 | |
---|
[1839] | 318 | REAL(wp), DIMENSION(1:3) :: rot_eigen_rad !< |
---|
| 319 | REAL(wp), DIMENSION(1:3) :: rot_eigen_azi !< |
---|
| 320 | REAL(wp), DIMENSION(1:3) :: rot_eigen_nor !< |
---|
| 321 | REAL(wp), DIMENSION(1:3) :: re !< |
---|
| 322 | REAL(wp), DIMENSION(1:3) :: rea !< |
---|
| 323 | REAL(wp), DIMENSION(1:3) :: ren !< |
---|
| 324 | REAL(wp), DIMENSION(1:3) :: rote !< |
---|
| 325 | REAL(wp), DIMENSION(1:3) :: rota !< |
---|
| 326 | REAL(wp), DIMENSION(1:3) :: rotn !< |
---|
[1819] | 327 | |
---|
[1839] | 328 | ! |
---|
| 329 | !-- Fixed variables for the speed controller |
---|
[1819] | 330 | |
---|
[1912] | 331 | LOGICAL :: start_up = .TRUE. !< |
---|
[1864] | 332 | |
---|
[1912] | 333 | REAL(wp) :: Fcorner !< corner freq for the controller low pass filter |
---|
| 334 | REAL(wp) :: min_reg25 !< min region 2.5 |
---|
| 335 | REAL(wp) :: om_rate !< rotor speed change |
---|
| 336 | REAL(wp) :: slope15 !< slope in region 1.5 |
---|
| 337 | REAL(wp) :: slope25 !< slope in region 2.5 |
---|
| 338 | REAL(wp) :: trq_rate !< torque change |
---|
| 339 | REAL(wp) :: vs_sysp !< |
---|
| 340 | REAL(wp) :: lp_coeff !< coeff for the controller low pass filter |
---|
[1819] | 341 | |
---|
[1864] | 342 | REAL(wp), DIMENSION(100) :: omega_rot_l = 0.0_wp !< local rot speed [rad/s] |
---|
[2563] | 343 | |
---|
[1839] | 344 | ! |
---|
| 345 | !-- Fixed variables for the yaw controller |
---|
[1819] | 346 | |
---|
[1912] | 347 | REAL(wp), DIMENSION(:) , ALLOCATABLE :: yawdir !< direction to yaw |
---|
| 348 | REAL(wp), DIMENSION(:) , ALLOCATABLE :: phi_yaw_l !< local (cpu) yaw angle |
---|
| 349 | REAL(wp), DIMENSION(:) , ALLOCATABLE :: wd30_l !< local (cpu) long running avg of the wd |
---|
| 350 | REAL(wp), DIMENSION(:) , ALLOCATABLE :: wd2_l !< local (cpu) short running avg of the wd |
---|
| 351 | REAL(wp), DIMENSION(:) , ALLOCATABLE :: wdir !< wind direction at hub |
---|
| 352 | REAL(wp), DIMENSION(:) , ALLOCATABLE :: u_inflow !< wind speed at hub |
---|
| 353 | REAL(wp), DIMENSION(:) , ALLOCATABLE :: wdir_l !< |
---|
| 354 | REAL(wp), DIMENSION(:) , ALLOCATABLE :: u_inflow_l !< |
---|
| 355 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wd30 !< |
---|
| 356 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wd2 !< |
---|
| 357 | LOGICAL, DIMENSION(1:100) :: doyaw = .FALSE. !< |
---|
| 358 | INTEGER(iwp) :: WDLON !< |
---|
| 359 | INTEGER(iwp) :: WDSHO !< |
---|
[1819] | 360 | |
---|
[2563] | 361 | ! |
---|
| 362 | !-- Variables that have to be saved in the binary file for restarts |
---|
| 363 | REAL(wp), DIMENSION(1:100) :: pitch_add_old = 0.0_wp !< old constant pitch angle |
---|
| 364 | REAL(wp), DIMENSION(1:100) :: omega_gen = 0.0_wp !< curr. generator speed |
---|
| 365 | REAL(wp), DIMENSION(1:100) :: omega_gen_f = 0.0_wp !< filtered generator speed |
---|
| 366 | REAL(wp), DIMENSION(1:100) :: omega_gen_old = 0.0_wp !< last generator speed |
---|
| 367 | REAL(wp), DIMENSION(1:100) :: omega_gen_f_old = 0.0_wp !< last filtered generator speed |
---|
| 368 | REAL(wp), DIMENSION(1:100) :: torque_gen = 0.0_wp !< generator torque |
---|
| 369 | REAL(wp), DIMENSION(1:100) :: torque_gen_old = 0.0_wp !< last generator torque |
---|
[1819] | 370 | |
---|
[2563] | 371 | |
---|
[1839] | 372 | SAVE |
---|
[1819] | 373 | |
---|
[1839] | 374 | |
---|
| 375 | INTERFACE wtm_parin |
---|
| 376 | MODULE PROCEDURE wtm_parin |
---|
| 377 | END INTERFACE wtm_parin |
---|
[2563] | 378 | |
---|
| 379 | INTERFACE wtm_write_restart_data |
---|
| 380 | MODULE PROCEDURE wtm_write_restart_data |
---|
| 381 | END INTERFACE wtm_write_restart_data |
---|
| 382 | |
---|
| 383 | INTERFACE wtm_read_restart_data |
---|
| 384 | MODULE PROCEDURE wtm_read_restart_data |
---|
| 385 | END INTERFACE wtm_read_restart_data |
---|
[1912] | 386 | |
---|
| 387 | INTERFACE wtm_check_parameters |
---|
| 388 | MODULE PROCEDURE wtm_check_parameters |
---|
| 389 | END INTERFACE wtm_check_parameters |
---|
[1839] | 390 | |
---|
| 391 | INTERFACE wtm_init_arrays |
---|
| 392 | MODULE PROCEDURE wtm_init_arrays |
---|
| 393 | END INTERFACE wtm_init_arrays |
---|
| 394 | |
---|
| 395 | INTERFACE wtm_init |
---|
| 396 | MODULE PROCEDURE wtm_init |
---|
| 397 | END INTERFACE wtm_init |
---|
| 398 | |
---|
[1912] | 399 | INTERFACE wtm_read_blade_tables |
---|
| 400 | MODULE PROCEDURE wtm_read_blade_tables |
---|
| 401 | END INTERFACE wtm_read_blade_tables |
---|
[1839] | 402 | |
---|
[1819] | 403 | INTERFACE wtm_forces |
---|
| 404 | MODULE PROCEDURE wtm_forces |
---|
[2553] | 405 | END INTERFACE wtm_forces |
---|
| 406 | |
---|
| 407 | INTERFACE wtm_yawcontrol |
---|
[1839] | 408 | MODULE PROCEDURE wtm_yawcontrol |
---|
[2553] | 409 | END INTERFACE wtm_yawcontrol |
---|
[1864] | 410 | |
---|
| 411 | INTERFACE wtm_rotate_rotor |
---|
| 412 | MODULE PROCEDURE wtm_rotate_rotor |
---|
| 413 | END INTERFACE wtm_rotate_rotor |
---|
| 414 | |
---|
| 415 | INTERFACE wtm_speed_control |
---|
| 416 | MODULE PROCEDURE wtm_init_speed_control |
---|
| 417 | MODULE PROCEDURE wtm_speed_control |
---|
| 418 | END INTERFACE wtm_speed_control |
---|
[1819] | 419 | |
---|
| 420 | INTERFACE wtm_tendencies |
---|
| 421 | MODULE PROCEDURE wtm_tendencies |
---|
| 422 | MODULE PROCEDURE wtm_tendencies_ij |
---|
| 423 | END INTERFACE wtm_tendencies |
---|
[1864] | 424 | |
---|
| 425 | |
---|
[1912] | 426 | PUBLIC wtm_check_parameters, wtm_forces, wtm_init, wtm_init_arrays, & |
---|
[2563] | 427 | wtm_parin, wtm_write_restart_data, wtm_tendencies |
---|
[1819] | 428 | |
---|
[2563] | 429 | |
---|
[1819] | 430 | CONTAINS |
---|
| 431 | |
---|
| 432 | |
---|
| 433 | !------------------------------------------------------------------------------! |
---|
| 434 | ! Description: |
---|
| 435 | ! ------------ |
---|
[1839] | 436 | !> Parin for &wind_turbine_par for wind turbine model |
---|
[1819] | 437 | !------------------------------------------------------------------------------! |
---|
[1839] | 438 | SUBROUTINE wtm_parin |
---|
[1819] | 439 | |
---|
| 440 | |
---|
| 441 | IMPLICIT NONE |
---|
[1912] | 442 | |
---|
| 443 | INTEGER(iwp) :: ierrn !< |
---|
[1819] | 444 | |
---|
[1839] | 445 | CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file |
---|
| 446 | |
---|
| 447 | NAMELIST /wind_turbine_par/ air_dens, dtow, gear_eff, gear_ratio, & |
---|
[1864] | 448 | gen_eff, inertia_gen, inertia_rot, max_miss, & |
---|
[1839] | 449 | max_torque_gen, max_trq_rate, min_miss, & |
---|
[1864] | 450 | min_reg15, min_reg2, nairfoils, nturbines, & |
---|
[1912] | 451 | omega_rot, phi_yaw, pitch_add, pitch_control,& |
---|
[1839] | 452 | rated_genspeed, rated_power, rcx, rcy, rcz, & |
---|
| 453 | rnac, rr, segment_length, segment_width, & |
---|
[1864] | 454 | slope2, speed_control, tilt, time_turbine_on,& |
---|
[2349] | 455 | turb_cd_nacelle, turb_cd_tower, pitch_rate, & |
---|
[2152] | 456 | yaw_control, yaw_speed, tl_cor |
---|
[1839] | 457 | |
---|
[1819] | 458 | ! |
---|
[1839] | 459 | !-- Try to find wind turbine model package |
---|
| 460 | REWIND ( 11 ) |
---|
| 461 | line = ' ' |
---|
| 462 | DO WHILE ( INDEX( line, '&wind_turbine_par' ) == 0 ) |
---|
| 463 | READ ( 11, '(A)', END=10 ) line |
---|
| 464 | ENDDO |
---|
| 465 | BACKSPACE ( 11 ) |
---|
| 466 | |
---|
| 467 | ! |
---|
| 468 | !-- Read user-defined namelist |
---|
[1912] | 469 | READ ( 11, wind_turbine_par, IOSTAT=ierrn ) |
---|
[2349] | 470 | |
---|
[1912] | 471 | IF ( ierrn < 0 ) THEN |
---|
[2349] | 472 | message_string = 'no wind_turbine_par-NAMELIST found: ' // & |
---|
| 473 | 'End of file has reached' |
---|
[2322] | 474 | CALL message( 'wtm_parin', 'PA0460', 1, 2, 0, 6, 0 ) |
---|
[2349] | 475 | ELSEIF ( ierrn > 0 ) THEN |
---|
| 476 | message_string = 'errors in wind_turbine_par-NAMELIST: ' // & |
---|
| 477 | 'some variables for steering may not be properly set' |
---|
| 478 | CALL message( 'wtm_parin', 'PA0466', 1, 2, 0, 6, 0 ) |
---|
[1912] | 479 | ENDIF |
---|
| 480 | |
---|
[1839] | 481 | ! |
---|
| 482 | !-- Set flag that indicates that the wind turbine model is switched on |
---|
| 483 | wind_turbine = .TRUE. |
---|
| 484 | |
---|
[2563] | 485 | IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN |
---|
| 486 | CALL wtm_read_restart_data |
---|
| 487 | ENDIF |
---|
| 488 | |
---|
[1912] | 489 | 10 CONTINUE ! TBD Change from continue, mit ierrn machen |
---|
[1839] | 490 | |
---|
| 491 | |
---|
| 492 | END SUBROUTINE wtm_parin |
---|
| 493 | |
---|
[2563] | 494 | |
---|
| 495 | !------------------------------------------------------------------------------! |
---|
| 496 | ! Description: |
---|
| 497 | ! ------------ |
---|
| 498 | !> This routine writes the respective restart data. |
---|
| 499 | !------------------------------------------------------------------------------! |
---|
| 500 | SUBROUTINE wtm_write_restart_data |
---|
| 501 | |
---|
| 502 | IMPLICIT NONE |
---|
| 503 | |
---|
| 504 | WRITE ( 14 ) 'omega_gen ' |
---|
| 505 | WRITE ( 14 ) omega_gen |
---|
| 506 | WRITE ( 14 ) 'omega_gen_f ' |
---|
| 507 | WRITE ( 14 ) omega_gen_f |
---|
| 508 | WRITE ( 14 ) 'omega_gen_f_old ' |
---|
| 509 | WRITE ( 14 ) omega_gen_f_old |
---|
| 510 | WRITE ( 14 ) 'omega_gen_old ' |
---|
| 511 | WRITE ( 14 ) omega_gen_old |
---|
| 512 | WRITE ( 14 ) 'omega_rot ' |
---|
| 513 | WRITE ( 14 ) omega_rot |
---|
| 514 | WRITE ( 14 ) 'phi_yaw ' |
---|
| 515 | WRITE ( 14 ) phi_yaw(:) |
---|
| 516 | WRITE ( 14 ) 'pitch_add ' |
---|
| 517 | WRITE ( 14 ) pitch_add |
---|
| 518 | WRITE ( 14 ) 'pitch_add_old ' |
---|
| 519 | WRITE ( 14 ) pitch_add_old |
---|
| 520 | WRITE ( 14 ) 'torque_gen ' |
---|
| 521 | WRITE ( 14 ) torque_gen |
---|
| 522 | WRITE ( 14 ) 'torque_gen_old ' |
---|
| 523 | WRITE ( 14 ) torque_gen_old |
---|
| 524 | |
---|
| 525 | WRITE ( 14 ) '*** end wtm *** ' |
---|
| 526 | |
---|
| 527 | END SUBROUTINE wtm_write_restart_data |
---|
| 528 | |
---|
| 529 | |
---|
| 530 | !------------------------------------------------------------------------------! |
---|
| 531 | ! Description: |
---|
| 532 | ! ------------ |
---|
| 533 | !> This routine reads the respective restart data. |
---|
| 534 | !------------------------------------------------------------------------------! |
---|
| 535 | SUBROUTINE wtm_read_restart_data |
---|
| 536 | |
---|
| 537 | |
---|
| 538 | IMPLICIT NONE |
---|
| 539 | |
---|
| 540 | CHARACTER (LEN=30) :: variable_chr !< dummy variable to read string |
---|
| 541 | |
---|
| 542 | |
---|
| 543 | READ ( 13 ) variable_chr |
---|
| 544 | DO WHILE ( TRIM( variable_chr ) /= '*** end wtm ***' ) |
---|
| 545 | |
---|
| 546 | SELECT CASE ( TRIM( variable_chr ) ) |
---|
| 547 | |
---|
| 548 | CASE ( 'omega_gen' ) |
---|
| 549 | READ ( 13 ) omega_gen |
---|
| 550 | CASE ( 'omega_gen_f' ) |
---|
| 551 | READ ( 13 ) omega_gen_f |
---|
| 552 | CASE ( 'omega_gen_old' ) |
---|
| 553 | READ ( 13 ) omega_gen_old |
---|
| 554 | CASE ( 'omega_gen_f_old' ) |
---|
| 555 | READ ( 13 ) omega_gen_f_old |
---|
| 556 | CASE ( 'omega_rot' ) |
---|
| 557 | READ ( 13 ) omega_rot |
---|
| 558 | CASE ( 'phi_yaw' ) |
---|
| 559 | READ ( 13 ) phi_yaw |
---|
| 560 | CASE ( 'pitch_add' ) |
---|
| 561 | READ ( 13 ) pitch_add |
---|
| 562 | CASE ( 'pitch_add_old' ) |
---|
| 563 | READ ( 13 ) pitch_add_old |
---|
| 564 | CASE ( 'torque_gen' ) |
---|
| 565 | READ ( 13 ) torque_gen |
---|
| 566 | CASE ( 'torque_gen_old' ) |
---|
| 567 | READ ( 13 ) torque_gen_old |
---|
| 568 | |
---|
| 569 | END SELECT |
---|
| 570 | |
---|
| 571 | READ ( 13 ) variable_chr |
---|
| 572 | |
---|
| 573 | ENDDO |
---|
| 574 | |
---|
| 575 | END SUBROUTINE wtm_read_restart_data |
---|
| 576 | |
---|
| 577 | |
---|
| 578 | !------------------------------------------------------------------------------! |
---|
| 579 | ! Description: |
---|
| 580 | ! ------------ |
---|
| 581 | !> Check namelist parameter |
---|
| 582 | !------------------------------------------------------------------------------! |
---|
[1912] | 583 | SUBROUTINE wtm_check_parameters |
---|
| 584 | |
---|
| 585 | |
---|
| 586 | IMPLICIT NONE |
---|
| 587 | |
---|
| 588 | IF ( ( .NOT.speed_control ) .AND. pitch_control ) THEN |
---|
| 589 | message_string = 'pitch_control = .TRUE. requires '// & |
---|
| 590 | 'speed_control = .TRUE.' |
---|
[2322] | 591 | CALL message( 'wtm_check_parameters', 'PA0461', 1, 2, 0, 6, 0 ) |
---|
[1912] | 592 | ENDIF |
---|
| 593 | |
---|
[2322] | 594 | IF ( ANY( omega_rot(1:nturbines) < 0.0 ) ) THEN |
---|
| 595 | message_string = 'omega_rot < 0.0, Please set omega_rot to ' // & |
---|
[2410] | 596 | 'a value larger than or equal to zero' |
---|
[2322] | 597 | CALL message( 'wtm_check_parameters', 'PA0462', 1, 2, 0, 6, 0 ) |
---|
[1912] | 598 | ENDIF |
---|
| 599 | |
---|
| 600 | |
---|
| 601 | IF ( ANY( rcx(1:nturbines) == 9999999.9_wp ) .OR. & |
---|
| 602 | ANY( rcy(1:nturbines) == 9999999.9_wp ) .OR. & |
---|
| 603 | ANY( rcz(1:nturbines) == 9999999.9_wp ) ) THEN |
---|
| 604 | |
---|
| 605 | message_string = 'rcx, rcy, rcz ' // & |
---|
| 606 | 'have to be given for each turbine.' |
---|
[2322] | 607 | CALL message( 'wtm_check_parameters', 'PA0463', 1, 2, 0, 6, 0 ) |
---|
[1912] | 608 | |
---|
| 609 | ENDIF |
---|
| 610 | |
---|
| 611 | |
---|
| 612 | END SUBROUTINE wtm_check_parameters |
---|
| 613 | |
---|
| 614 | |
---|
[1839] | 615 | !------------------------------------------------------------------------------! |
---|
| 616 | ! Description: |
---|
| 617 | ! ------------ |
---|
| 618 | !> Allocate wind turbine model arrays |
---|
| 619 | !------------------------------------------------------------------------------! |
---|
| 620 | SUBROUTINE wtm_init_arrays |
---|
| 621 | |
---|
| 622 | |
---|
| 623 | IMPLICIT NONE |
---|
| 624 | |
---|
[1864] | 625 | REAL(wp) :: delta_r_factor !< |
---|
| 626 | REAL(wp) :: delta_r_init !< |
---|
| 627 | |
---|
[1839] | 628 | ! |
---|
| 629 | !-- To be able to allocate arrays with dimension of rotor rings and segments, |
---|
[1819] | 630 | !-- the maximum possible numbers of rings and segments have to be calculated: |
---|
| 631 | |
---|
| 632 | ALLOCATE( nrings(1:nturbines) ) |
---|
| 633 | ALLOCATE( delta_r(1:nturbines) ) |
---|
| 634 | |
---|
| 635 | nrings(:) = 0 |
---|
| 636 | delta_r(:) = 0.0_wp |
---|
| 637 | |
---|
| 638 | ! |
---|
| 639 | !-- Thickness (radial) of each ring and length (tangential) of each segment: |
---|
| 640 | delta_r_factor = segment_width |
---|
| 641 | delta_t_factor = segment_length |
---|
| 642 | delta_r_init = delta_r_factor * MIN( dx, dy, dz) |
---|
| 643 | |
---|
| 644 | DO inot = 1, nturbines |
---|
| 645 | ! |
---|
| 646 | !-- Determine number of rings: |
---|
| 647 | nrings(inot) = NINT( rr(inot) / delta_r_init ) |
---|
| 648 | |
---|
| 649 | delta_r(inot) = rr(inot) / nrings(inot) |
---|
| 650 | |
---|
| 651 | ENDDO |
---|
| 652 | |
---|
| 653 | nrings_max = MAXVAL(nrings) |
---|
| 654 | |
---|
| 655 | ALLOCATE( nsegs(1:nrings_max,1:nturbines) ) |
---|
| 656 | ALLOCATE( nsegs_total(1:nturbines) ) |
---|
| 657 | |
---|
| 658 | nsegs(:,:) = 0 |
---|
| 659 | nsegs_total(:) = 0 |
---|
| 660 | |
---|
| 661 | |
---|
| 662 | DO inot = 1, nturbines |
---|
| 663 | DO ring = 1, nrings(inot) |
---|
| 664 | ! |
---|
| 665 | !-- Determine number of segments for each ring: |
---|
[1839] | 666 | nsegs(ring,inot) = MAX( 8, CEILING( delta_r_factor * pi * & |
---|
| 667 | ( 2.0_wp * ring - 1.0_wp ) / & |
---|
| 668 | delta_t_factor ) ) |
---|
[1819] | 669 | ENDDO |
---|
| 670 | ! |
---|
| 671 | !-- Total sum of all rotor segments: |
---|
[1839] | 672 | nsegs_total(inot) = SUM( nsegs(:,inot) ) |
---|
[1819] | 673 | |
---|
| 674 | ENDDO |
---|
| 675 | |
---|
| 676 | ! |
---|
| 677 | !-- Maximum number of segments per ring: |
---|
| 678 | nsegs_max = MAXVAL(nsegs) |
---|
| 679 | |
---|
[1864] | 680 | !! |
---|
| 681 | !!-- TODO: Folgendes im Header ausgeben! |
---|
| 682 | ! IF ( myid == 0 ) THEN |
---|
| 683 | ! PRINT*, 'nrings(1) = ', nrings(1) |
---|
| 684 | ! PRINT*, '--------------------------------------------------' |
---|
| 685 | ! PRINT*, 'nsegs(:,1) = ', nsegs(:,1) |
---|
| 686 | ! PRINT*, '--------------------------------------------------' |
---|
| 687 | ! PRINT*, 'nrings_max = ', nrings_max |
---|
| 688 | ! PRINT*, 'nsegs_max = ', nsegs_max |
---|
| 689 | ! PRINT*, 'nsegs_total(1) = ', nsegs_total(1) |
---|
| 690 | ! ENDIF |
---|
[1819] | 691 | |
---|
| 692 | |
---|
| 693 | ! |
---|
| 694 | !-- Allocate 1D arrays (dimension = number of turbines) |
---|
| 695 | ALLOCATE( i_hub(1:nturbines) ) |
---|
| 696 | ALLOCATE( i_smear(1:nturbines) ) |
---|
| 697 | ALLOCATE( j_hub(1:nturbines) ) |
---|
| 698 | ALLOCATE( j_smear(1:nturbines) ) |
---|
| 699 | ALLOCATE( k_hub(1:nturbines) ) |
---|
| 700 | ALLOCATE( k_smear(1:nturbines) ) |
---|
| 701 | ALLOCATE( torque_total(1:nturbines) ) |
---|
[1912] | 702 | ALLOCATE( thrust_rotor(1:nturbines) ) |
---|
[1819] | 703 | |
---|
| 704 | ! |
---|
| 705 | !-- Allocation of the 1D arrays for yaw control |
---|
| 706 | ALLOCATE( yawdir(1:nturbines) ) |
---|
| 707 | ALLOCATE( u_inflow(1:nturbines) ) |
---|
| 708 | ALLOCATE( wdir(1:nturbines) ) |
---|
| 709 | ALLOCATE( u_inflow_l(1:nturbines) ) |
---|
| 710 | ALLOCATE( wdir_l(1:nturbines) ) |
---|
| 711 | ALLOCATE( phi_yaw_l(1:nturbines) ) |
---|
| 712 | |
---|
| 713 | ! |
---|
| 714 | !-- Allocate 1D arrays (dimension = number of rotor segments) |
---|
| 715 | ALLOCATE( alpha_attack(1:nsegs_max) ) |
---|
| 716 | ALLOCATE( chord(1:nsegs_max) ) |
---|
| 717 | ALLOCATE( phi_rel(1:nsegs_max) ) |
---|
| 718 | ALLOCATE( thrust_seg(1:nsegs_max) ) |
---|
| 719 | ALLOCATE( torque_seg_y(1:nsegs_max) ) |
---|
| 720 | ALLOCATE( torque_seg_z(1:nsegs_max) ) |
---|
[1839] | 721 | ALLOCATE( turb_cd(1:nsegs_max) ) |
---|
| 722 | ALLOCATE( turb_cl(1:nsegs_max) ) |
---|
[1819] | 723 | ALLOCATE( vrel(1:nsegs_max) ) |
---|
| 724 | ALLOCATE( vtheta(1:nsegs_max) ) |
---|
| 725 | |
---|
| 726 | ! |
---|
| 727 | !-- Allocate 2D arrays (dimension = number of rotor rings and segments) |
---|
| 728 | ALLOCATE( rbx(1:nrings_max,1:nsegs_max) ) |
---|
| 729 | ALLOCATE( rby(1:nrings_max,1:nsegs_max) ) |
---|
| 730 | ALLOCATE( rbz(1:nrings_max,1:nsegs_max) ) |
---|
| 731 | ALLOCATE( thrust_ring(1:nrings_max,1:nsegs_max) ) |
---|
| 732 | ALLOCATE( torque_ring_y(1:nrings_max,1:nsegs_max) ) |
---|
| 733 | ALLOCATE( torque_ring_z(1:nrings_max,1:nsegs_max) ) |
---|
| 734 | |
---|
| 735 | ! |
---|
| 736 | !-- Allocate additional 2D arrays |
---|
| 737 | ALLOCATE( rotx(1:nturbines,1:3) ) |
---|
| 738 | ALLOCATE( roty(1:nturbines,1:3) ) |
---|
| 739 | ALLOCATE( rotz(1:nturbines,1:3) ) |
---|
| 740 | |
---|
| 741 | ! |
---|
| 742 | !-- Allocate 3D arrays (dimension = number of grid points) |
---|
[1912] | 743 | ALLOCATE( nac_cd_surf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 744 | ALLOCATE( rot_tend_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 745 | ALLOCATE( rot_tend_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 746 | ALLOCATE( rot_tend_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
[1819] | 747 | ALLOCATE( thrust(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 748 | ALLOCATE( torque_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 749 | ALLOCATE( torque_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
[1912] | 750 | ALLOCATE( tow_cd_surf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
[1819] | 751 | |
---|
| 752 | ! |
---|
| 753 | !-- Allocate additional 3D arrays |
---|
| 754 | ALLOCATE( u_int(1:nturbines,1:nrings_max,1:nsegs_max) ) |
---|
| 755 | ALLOCATE( u_int_1_l(1:nturbines,1:nrings_max,1:nsegs_max) ) |
---|
| 756 | ALLOCATE( v_int(1:nturbines,1:nrings_max,1:nsegs_max) ) |
---|
| 757 | ALLOCATE( v_int_1_l(1:nturbines,1:nrings_max,1:nsegs_max) ) |
---|
| 758 | ALLOCATE( w_int(1:nturbines,1:nrings_max,1:nsegs_max) ) |
---|
| 759 | ALLOCATE( w_int_1_l(1:nturbines,1:nrings_max,1:nsegs_max) ) |
---|
| 760 | |
---|
| 761 | ! |
---|
| 762 | !-- All of the arrays are initialized with a value of zero: |
---|
| 763 | i_hub(:) = 0 |
---|
| 764 | i_smear(:) = 0 |
---|
| 765 | j_hub(:) = 0 |
---|
| 766 | j_smear(:) = 0 |
---|
| 767 | k_hub(:) = 0 |
---|
| 768 | k_smear(:) = 0 |
---|
[1912] | 769 | |
---|
[1819] | 770 | torque_total(:) = 0.0_wp |
---|
[1912] | 771 | thrust_rotor(:) = 0.0_wp |
---|
[1819] | 772 | |
---|
[2563] | 773 | IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN |
---|
| 774 | omega_gen(:) = 0.0_wp |
---|
| 775 | omega_gen_old(:) = 0.0_wp |
---|
| 776 | omega_gen_f(:) = 0.0_wp |
---|
| 777 | omega_gen_f_old(:) = 0.0_wp |
---|
| 778 | pitch_add_old(:) = 0.0_wp |
---|
| 779 | torque_gen(:) = 0.0_wp |
---|
| 780 | torque_gen_old(:) = 0.0_wp |
---|
| 781 | ENDIF |
---|
| 782 | |
---|
[1819] | 783 | yawdir(:) = 0.0_wp |
---|
| 784 | wdir(:) = 0.0_wp |
---|
| 785 | u_inflow(:) = 0.0_wp |
---|
| 786 | |
---|
| 787 | ! |
---|
| 788 | !-- Allocate 1D arrays (dimension = number of rotor segments) |
---|
| 789 | alpha_attack(:) = 0.0_wp |
---|
| 790 | chord(:) = 0.0_wp |
---|
| 791 | phi_rel(:) = 0.0_wp |
---|
| 792 | thrust_seg(:) = 0.0_wp |
---|
| 793 | torque_seg_y(:) = 0.0_wp |
---|
| 794 | torque_seg_z(:) = 0.0_wp |
---|
[1864] | 795 | turb_cd(:) = 0.0_wp |
---|
| 796 | turb_cl(:) = 0.0_wp |
---|
[1819] | 797 | vrel(:) = 0.0_wp |
---|
| 798 | vtheta(:) = 0.0_wp |
---|
| 799 | |
---|
| 800 | rbx(:,:) = 0.0_wp |
---|
| 801 | rby(:,:) = 0.0_wp |
---|
| 802 | rbz(:,:) = 0.0_wp |
---|
| 803 | thrust_ring(:,:) = 0.0_wp |
---|
| 804 | torque_ring_y(:,:) = 0.0_wp |
---|
| 805 | torque_ring_z(:,:) = 0.0_wp |
---|
| 806 | |
---|
| 807 | rotx(:,:) = 0.0_wp |
---|
| 808 | roty(:,:) = 0.0_wp |
---|
| 809 | rotz(:,:) = 0.0_wp |
---|
| 810 | |
---|
[1912] | 811 | nac_cd_surf(:,:,:) = 0.0_wp |
---|
| 812 | rot_tend_x(:,:,:) = 0.0_wp |
---|
| 813 | rot_tend_y(:,:,:) = 0.0_wp |
---|
| 814 | rot_tend_z(:,:,:) = 0.0_wp |
---|
[1819] | 815 | thrust(:,:,:) = 0.0_wp |
---|
| 816 | torque_y(:,:,:) = 0.0_wp |
---|
| 817 | torque_z(:,:,:) = 0.0_wp |
---|
[1912] | 818 | tow_cd_surf(:,:,:) = 0.0_wp |
---|
[1819] | 819 | |
---|
| 820 | u_int(:,:,:) = 0.0_wp |
---|
| 821 | u_int_1_l(:,:,:) = 0.0_wp |
---|
| 822 | v_int(:,:,:) = 0.0_wp |
---|
| 823 | v_int_1_l(:,:,:) = 0.0_wp |
---|
| 824 | w_int(:,:,:) = 0.0_wp |
---|
| 825 | w_int_1_l(:,:,:) = 0.0_wp |
---|
| 826 | |
---|
| 827 | |
---|
[1839] | 828 | END SUBROUTINE wtm_init_arrays |
---|
[1819] | 829 | |
---|
| 830 | |
---|
| 831 | !------------------------------------------------------------------------------! |
---|
| 832 | ! Description: |
---|
| 833 | ! ------------ |
---|
[1839] | 834 | !> Initialization of the wind turbine model |
---|
[1819] | 835 | !------------------------------------------------------------------------------! |
---|
[1839] | 836 | SUBROUTINE wtm_init |
---|
[1819] | 837 | |
---|
[1839] | 838 | |
---|
[1819] | 839 | IMPLICIT NONE |
---|
| 840 | |
---|
| 841 | INTEGER(iwp) :: i !< running index |
---|
| 842 | INTEGER(iwp) :: j !< running index |
---|
| 843 | INTEGER(iwp) :: k !< running index |
---|
[1864] | 844 | |
---|
[1819] | 845 | ! |
---|
[1864] | 846 | !-- Help variables for the smearing function |
---|
| 847 | REAL(wp) :: eps_kernel !< |
---|
| 848 | |
---|
[1839] | 849 | ! |
---|
[1864] | 850 | !-- Help variables for calculation of the tower drag |
---|
| 851 | INTEGER(iwp) :: tower_n !< |
---|
| 852 | INTEGER(iwp) :: tower_s !< |
---|
[1839] | 853 | ! |
---|
[1864] | 854 | !-- Help variables for the calulaction of the nacelle drag |
---|
| 855 | INTEGER(iwp) :: i_ip !< |
---|
| 856 | INTEGER(iwp) :: i_ipg !< |
---|
| 857 | |
---|
| 858 | REAL(wp) :: yvalue |
---|
| 859 | REAL(wp) :: dy_int !< |
---|
| 860 | REAL(wp) :: dz_int !< |
---|
| 861 | |
---|
| 862 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: circle_points !< |
---|
| 863 | |
---|
| 864 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nacb !< |
---|
| 865 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nacl !< |
---|
| 866 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nacr !< |
---|
| 867 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nact !< |
---|
| 868 | |
---|
| 869 | ALLOCATE( index_nacb(1:nturbines) ) |
---|
| 870 | ALLOCATE( index_nacl(1:nturbines) ) |
---|
| 871 | ALLOCATE( index_nacr(1:nturbines) ) |
---|
| 872 | ALLOCATE( index_nact(1:nturbines) ) |
---|
[1819] | 873 | |
---|
| 874 | |
---|
[1864] | 875 | IF ( speed_control) THEN |
---|
| 876 | |
---|
| 877 | CALL wtm_speed_control |
---|
| 878 | |
---|
[2563] | 879 | IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN |
---|
| 880 | |
---|
| 881 | DO inot = 1, nturbines |
---|
| 882 | |
---|
| 883 | IF ( nxl > i_hub(inot) ) THEN |
---|
| 884 | torque_gen(inot) = 0.0_wp |
---|
| 885 | omega_gen_f(inot) = 0.0_wp |
---|
| 886 | omega_rot_l(inot) = 0.0_wp |
---|
| 887 | ENDIF |
---|
| 888 | |
---|
| 889 | IF ( nxr < i_hub(inot) ) THEN |
---|
| 890 | torque_gen(inot) = 0.0_wp |
---|
| 891 | omega_gen_f(inot) = 0.0_wp |
---|
| 892 | omega_rot_l(inot) = 0.0_wp |
---|
| 893 | ENDIF |
---|
| 894 | |
---|
| 895 | IF ( nys > j_hub(inot) ) THEN |
---|
| 896 | torque_gen(inot) = 0.0_wp |
---|
| 897 | omega_gen_f(inot) = 0.0_wp |
---|
| 898 | omega_rot_l(inot) = 0.0_wp |
---|
| 899 | ENDIF |
---|
| 900 | |
---|
| 901 | IF ( nyn < j_hub(inot) ) THEN |
---|
| 902 | torque_gen(inot) = 0.0_wp |
---|
| 903 | omega_gen_f(inot) = 0.0_wp |
---|
| 904 | omega_rot_l(inot) = 0.0_wp |
---|
| 905 | ENDIF |
---|
| 906 | |
---|
| 907 | IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) THEN |
---|
| 908 | IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) ) THEN |
---|
| 909 | |
---|
| 910 | omega_rot_l(inot) = omega_gen(inot) / gear_ratio |
---|
| 911 | |
---|
| 912 | ENDIF |
---|
| 913 | ENDIF |
---|
| 914 | |
---|
| 915 | END DO |
---|
| 916 | |
---|
| 917 | ENDIF |
---|
| 918 | |
---|
[1864] | 919 | ENDIF |
---|
| 920 | |
---|
[1839] | 921 | ! |
---|
[1819] | 922 | !------------------------------------------------------------------------------! |
---|
[1839] | 923 | !-- Calculation of parameters for the regularization kernel |
---|
| 924 | !-- (smearing of the forces) |
---|
[1819] | 925 | !------------------------------------------------------------------------------! |
---|
| 926 | ! |
---|
[1839] | 927 | !-- In the following, some of the required parameters for the smearing will |
---|
| 928 | !-- be calculated: |
---|
[1819] | 929 | |
---|
[1839] | 930 | !-- The kernel is set equal to twice the grid spacing which has turned out to |
---|
| 931 | !-- be a reasonable value (see e.g. Troldborg et al. (2013), Wind Energy, |
---|
[1819] | 932 | !-- DOI: 10.1002/we.1608): |
---|
| 933 | eps_kernel = 2.0_wp * dx |
---|
| 934 | ! |
---|
[1839] | 935 | !-- The zero point (eps_min) of the polynomial function must be the following |
---|
| 936 | !-- if the integral of the polynomial function (for values < eps_min) shall |
---|
| 937 | !-- be equal to the integral of the Gaussian function used before: |
---|
| 938 | eps_min = ( 105.0_wp / 32.0_wp )**( 1.0_wp / 3.0_wp ) * & |
---|
| 939 | pi**( 1.0_wp / 6.0_wp ) * eps_kernel |
---|
[1819] | 940 | ! |
---|
| 941 | !-- Square of eps_min: |
---|
| 942 | eps_min2 = eps_min**2 |
---|
| 943 | ! |
---|
| 944 | !-- Parameters in the polynomial function: |
---|
[1864] | 945 | pol_a = 1.0_wp / eps_min**4 |
---|
| 946 | pol_b = 2.0_wp / eps_min**2 |
---|
[1819] | 947 | ! |
---|
[1839] | 948 | !-- Normalization factor which is the inverse of the integral of the smearing |
---|
| 949 | !-- function: |
---|
| 950 | eps_factor = 105.0_wp / ( 32.0_wp * pi * eps_min**3 ) |
---|
| 951 | |
---|
[1864] | 952 | !-- Change tilt angle to rad: |
---|
| 953 | tilt = tilt * pi / 180.0_wp |
---|
| 954 | |
---|
[1819] | 955 | ! |
---|
[1864] | 956 | !-- Change yaw angle to rad: |
---|
[2563] | 957 | IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN |
---|
| 958 | phi_yaw(:) = phi_yaw(:) * pi / 180.0_wp |
---|
| 959 | ENDIF |
---|
[1819] | 960 | |
---|
[1864] | 961 | |
---|
[1819] | 962 | DO inot = 1, nturbines |
---|
| 963 | ! |
---|
[1864] | 964 | !-- Rotate the rotor coordinates in case yaw and tilt are defined |
---|
| 965 | CALL wtm_rotate_rotor( inot ) |
---|
[1819] | 966 | |
---|
| 967 | ! |
---|
| 968 | !-- Determine the indices of the hub height |
---|
| 969 | i_hub(inot) = INT( rcx(inot) / dx ) |
---|
| 970 | j_hub(inot) = INT( ( rcy(inot) + 0.5_wp * dy ) / dy ) |
---|
| 971 | k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz ) / dz ) |
---|
| 972 | |
---|
| 973 | ! |
---|
| 974 | !-- Determining the area to which the smearing of the forces is applied. |
---|
[1839] | 975 | !-- As smearing now is effectively applied only for distances smaller than |
---|
| 976 | !-- eps_min, the smearing area can be further limited and regarded as a |
---|
| 977 | !-- function of eps_min: |
---|
[1819] | 978 | i_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dx ) |
---|
| 979 | j_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dy ) |
---|
| 980 | k_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dz ) |
---|
[1864] | 981 | |
---|
[1819] | 982 | ENDDO |
---|
| 983 | |
---|
| 984 | ! |
---|
| 985 | !------------------------------------------------------------------------------! |
---|
[1839] | 986 | !-- Determine the area within each grid cell that overlaps with the area |
---|
| 987 | !-- of the nacelle and the tower (needed for calculation of the forces) |
---|
[1819] | 988 | !------------------------------------------------------------------------------! |
---|
| 989 | ! |
---|
| 990 | !-- Note: so far this is only a 2D version, in that the mean flow is |
---|
| 991 | !-- perpendicular to the rotor area. |
---|
| 992 | |
---|
| 993 | ! |
---|
| 994 | !-- Allocation of the array containing information on the intersection points |
---|
| 995 | !-- between rotor disk and the numerical grid: |
---|
| 996 | upper_end = ( ny + 1 ) * 10000 |
---|
| 997 | |
---|
| 998 | ALLOCATE( circle_points(1:2,1:upper_end) ) |
---|
[1839] | 999 | |
---|
| 1000 | circle_points(:,:) = 0.0_wp |
---|
[1819] | 1001 | |
---|
[1839] | 1002 | |
---|
| 1003 | DO inot = 1, nturbines ! loop over number of turbines |
---|
[1819] | 1004 | ! |
---|
[1839] | 1005 | !-- Determine the grid index (u-grid) that corresponds to the location of |
---|
| 1006 | !-- the rotor center (reduces the amount of calculations in the case that |
---|
| 1007 | !-- the mean flow is perpendicular to the rotor area): |
---|
[1819] | 1008 | i = i_hub(inot) |
---|
| 1009 | |
---|
| 1010 | ! |
---|
[1839] | 1011 | !-- Determine the left and the right edge of the nacelle (corresponding |
---|
| 1012 | !-- grid point indices): |
---|
[1819] | 1013 | index_nacl(inot) = INT( ( rcy(inot) - rnac(inot) + 0.5_wp * dy ) / dy ) |
---|
| 1014 | index_nacr(inot) = INT( ( rcy(inot) + rnac(inot) + 0.5_wp * dy ) / dy ) |
---|
| 1015 | ! |
---|
[1839] | 1016 | !-- Determine the bottom and the top edge of the nacelle (corresponding |
---|
| 1017 | !-- grid point indices).The grid point index has to be increased by 1, as |
---|
| 1018 | !-- the first level for the u-component (index 0) is situated below the |
---|
| 1019 | !-- surface. All points between z=0 and z=dz/s would already be contained |
---|
| 1020 | !-- in grid box 1. |
---|
[1819] | 1021 | index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz ) + 1 |
---|
| 1022 | index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz ) + 1 |
---|
| 1023 | |
---|
| 1024 | ! |
---|
| 1025 | !-- Determine the indices of the grid boxes containing the left and |
---|
[1864] | 1026 | !-- the right boundaries of the tower: |
---|
[1819] | 1027 | tower_n = ( rcy(inot) + 0.5_wp * dtow(inot) - 0.5_wp * dy ) / dy |
---|
| 1028 | tower_s = ( rcy(inot) - 0.5_wp * dtow(inot) - 0.5_wp * dy ) / dy |
---|
| 1029 | |
---|
| 1030 | ! |
---|
| 1031 | !-- Determine the fraction of the grid box area overlapping with the tower |
---|
[1864] | 1032 | !-- area and multiply it with the drag of the tower: |
---|
[1839] | 1033 | IF ( ( nxlg <= i ) .AND. ( nxrg >= i ) ) THEN |
---|
[1819] | 1034 | |
---|
| 1035 | DO j = nys, nyn |
---|
| 1036 | ! |
---|
[1839] | 1037 | !-- Loop from south to north boundary of tower |
---|
| 1038 | IF ( ( j >= tower_s ) .AND. ( j <= tower_n ) ) THEN |
---|
| 1039 | |
---|
[1819] | 1040 | DO k = nzb, nzt |
---|
| 1041 | |
---|
| 1042 | IF ( k == k_hub(inot) ) THEN |
---|
| 1043 | IF ( tower_n - tower_s >= 1 ) THEN |
---|
[1839] | 1044 | ! |
---|
[1819] | 1045 | !-- leftmost and rightmost grid box: |
---|
| 1046 | IF ( j == tower_s ) THEN |
---|
[1912] | 1047 | tow_cd_surf(k,j,i) = ( rcz(inot) - & |
---|
| 1048 | ( k_hub(inot) * dz - 0.5_wp * dz ) ) * & ! extension in z-direction |
---|
| 1049 | ( ( tower_s + 1.0_wp + 0.5_wp ) * dy - & |
---|
| 1050 | ( rcy(inot) - 0.5_wp * dtow(inot) ) ) * & ! extension in y-direction |
---|
| 1051 | turb_cd_tower(inot) |
---|
[1819] | 1052 | ELSEIF ( j == tower_n ) THEN |
---|
[1912] | 1053 | tow_cd_surf(k,j,i) = ( rcz(inot) - & |
---|
| 1054 | ( k_hub(inot) * dz - 0.5_wp * dz ) ) * & ! extension in z-direction |
---|
| 1055 | ( ( rcy(inot) + 0.5_wp * dtow(inot) ) - & |
---|
| 1056 | ( tower_n + 0.5_wp ) * dy ) * & ! extension in y-direction |
---|
| 1057 | turb_cd_tower(inot) |
---|
[1819] | 1058 | ! |
---|
| 1059 | !-- grid boxes inbetween |
---|
[1912] | 1060 | !-- (where tow_cd_surf = grid box area): |
---|
[1819] | 1061 | ELSE |
---|
[1912] | 1062 | tow_cd_surf(k,j,i) = ( rcz(inot) - & |
---|
| 1063 | ( k_hub(inot) * dz - 0.5_wp * dz ) ) * & |
---|
| 1064 | dy * turb_cd_tower(inot) |
---|
[1819] | 1065 | ENDIF |
---|
| 1066 | ! |
---|
| 1067 | !-- tower lies completely within one grid box: |
---|
| 1068 | ELSE |
---|
[1912] | 1069 | tow_cd_surf(k,j,i) = ( rcz(inot) - & |
---|
[1839] | 1070 | ( k_hub(inot) * dz - 0.5_wp * dz ) ) * & |
---|
[1912] | 1071 | dtow(inot) * turb_cd_tower(inot) |
---|
[1819] | 1072 | ENDIF |
---|
[1839] | 1073 | ! |
---|
| 1074 | !-- In case that k is smaller than k_hub the following actions |
---|
| 1075 | !-- are carried out: |
---|
[1819] | 1076 | ELSEIF ( k < k_hub(inot) ) THEN |
---|
| 1077 | |
---|
[1839] | 1078 | IF ( ( tower_n - tower_s ) >= 1 ) THEN |
---|
| 1079 | ! |
---|
[1819] | 1080 | !-- leftmost and rightmost grid box: |
---|
| 1081 | IF ( j == tower_s ) THEN |
---|
[1912] | 1082 | tow_cd_surf(k,j,i) = dz * ( & |
---|
[1839] | 1083 | ( tower_s + 1 + 0.5_wp ) * dy - & |
---|
| 1084 | ( rcy(inot) - 0.5_wp * dtow(inot) ) & |
---|
[1912] | 1085 | ) * turb_cd_tower(inot) |
---|
[1819] | 1086 | ELSEIF ( j == tower_n ) THEN |
---|
[1912] | 1087 | tow_cd_surf(k,j,i) = dz * ( & |
---|
[1839] | 1088 | ( rcy(inot) + 0.5_wp * dtow(inot) ) - & |
---|
| 1089 | ( tower_n + 0.5_wp ) * dy & |
---|
[1912] | 1090 | ) * turb_cd_tower(inot) |
---|
[1839] | 1091 | ! |
---|
| 1092 | !-- grid boxes inbetween |
---|
[1912] | 1093 | !-- (where tow_cd_surf = grid box area): |
---|
[1819] | 1094 | ELSE |
---|
[1912] | 1095 | tow_cd_surf(k,j,i) = dz * dy * turb_cd_tower(inot) |
---|
[1819] | 1096 | ENDIF |
---|
[1839] | 1097 | ! |
---|
[1819] | 1098 | !-- tower lies completely within one grid box: |
---|
| 1099 | ELSE |
---|
[1912] | 1100 | tow_cd_surf(k,j,i) = dz * dtow(inot) * & |
---|
| 1101 | turb_cd_tower(inot) |
---|
[1839] | 1102 | ENDIF ! end if larger than grid box |
---|
| 1103 | |
---|
| 1104 | ENDIF ! end if k == k_hub |
---|
| 1105 | |
---|
| 1106 | ENDDO ! end loop over k |
---|
| 1107 | |
---|
| 1108 | ENDIF ! end if inside north and south boundary of tower |
---|
| 1109 | |
---|
| 1110 | ENDDO ! end loop over j |
---|
| 1111 | |
---|
| 1112 | ENDIF ! end if hub inside domain + ghostpoints |
---|
[1819] | 1113 | |
---|
[1839] | 1114 | |
---|
[1912] | 1115 | CALL exchange_horiz( tow_cd_surf, nbgp ) |
---|
[1819] | 1116 | |
---|
[1839] | 1117 | ! |
---|
[1864] | 1118 | !-- Calculation of the nacelle area |
---|
| 1119 | !-- CAUTION: Currently disabled due to segmentation faults on the FLOW HPC |
---|
| 1120 | !-- cluster (Oldenburg) |
---|
| 1121 | !! |
---|
| 1122 | !!-- Tabulate the points on the circle that are required in the following for |
---|
| 1123 | !!-- the calculation of the Riemann integral (node points; they are called |
---|
| 1124 | !!-- circle_points in the following): |
---|
| 1125 | ! |
---|
| 1126 | ! dy_int = dy / 10000.0_wp |
---|
| 1127 | ! |
---|
| 1128 | ! DO i_ip = 1, upper_end |
---|
| 1129 | ! yvalue = dy_int * ( i_ip - 0.5_wp ) + 0.5_wp * dy !<--- segmentation fault |
---|
| 1130 | ! sqrt_arg = rnac(inot)**2 - ( yvalue - rcy(inot) )**2 !<--- segmentation fault |
---|
| 1131 | ! IF ( sqrt_arg >= 0.0_wp ) THEN |
---|
| 1132 | !! |
---|
| 1133 | !!-- bottom intersection point |
---|
| 1134 | ! circle_points(1,i_ip) = rcz(inot) - SQRT( sqrt_arg ) |
---|
| 1135 | !! |
---|
| 1136 | !!-- top intersection point |
---|
| 1137 | ! circle_points(2,i_ip) = rcz(inot) + SQRT( sqrt_arg ) !<--- segmentation fault |
---|
| 1138 | ! ELSE |
---|
| 1139 | ! circle_points(:,i_ip) = -111111 !<--- segmentation fault |
---|
| 1140 | ! ENDIF |
---|
| 1141 | ! ENDDO |
---|
| 1142 | ! |
---|
| 1143 | ! |
---|
| 1144 | ! DO j = nys, nyn |
---|
| 1145 | !! |
---|
| 1146 | !!-- In case that the grid box is located completely outside the nacelle |
---|
| 1147 | !!-- (y) it can automatically be stated that there is no overlap between |
---|
| 1148 | !!-- the grid box and the nacelle and consequently we can set |
---|
[1912] | 1149 | !!-- nac_cd_surf(:,j,i) = 0.0: |
---|
[1864] | 1150 | ! IF ( ( j >= index_nacl(inot) ) .AND. ( j <= index_nacr(inot) ) ) THEN |
---|
| 1151 | ! DO k = nzb+1, nzt |
---|
| 1152 | !! |
---|
| 1153 | !!-- In case that the grid box is located completely outside the |
---|
| 1154 | !!-- nacelle (z) it can automatically be stated that there is no |
---|
| 1155 | !!-- overlap between the grid box and the nacelle and consequently |
---|
[1912] | 1156 | !!-- we can set nac_cd_surf(k,j,i) = 0.0: |
---|
[1864] | 1157 | ! IF ( ( k >= index_nacb(inot) ) .OR. & |
---|
| 1158 | ! ( k <= index_nact(inot) ) ) THEN |
---|
| 1159 | !! |
---|
| 1160 | !!-- For all other cases Riemann integrals are calculated. |
---|
| 1161 | !!-- Here, the points on the circle that have been determined |
---|
| 1162 | !!-- above are used in order to calculate the overlap between the |
---|
| 1163 | !!-- gridbox and the nacelle area (area approached by 10000 |
---|
| 1164 | !!-- rectangulars dz_int * dy_int): |
---|
| 1165 | ! DO i_ipg = 1, 10000 |
---|
| 1166 | ! dz_int = dz |
---|
| 1167 | ! i_ip = j * 10000 + i_ipg |
---|
| 1168 | !! |
---|
| 1169 | !!-- Determine the vertical extension dz_int of the circle |
---|
| 1170 | !!-- within the current grid box: |
---|
| 1171 | ! IF ( ( circle_points(2,i_ip) < zw(k) ) .AND. & !<--- segmentation fault |
---|
| 1172 | ! ( circle_points(2,i_ip) >= zw(k-1) ) ) THEN |
---|
| 1173 | ! dz_int = dz_int - & !<--- segmentation fault |
---|
| 1174 | ! ( zw(k) - circle_points(2,i_ip) ) |
---|
| 1175 | ! ENDIF |
---|
| 1176 | ! IF ( ( circle_points(1,i_ip) <= zw(k) ) .AND. & !<--- segmentation fault |
---|
| 1177 | ! ( circle_points(1,i_ip) > zw(k-1) ) ) THEN |
---|
| 1178 | ! dz_int = dz_int - & |
---|
| 1179 | ! ( circle_points(1,i_ip) - zw(k-1) ) |
---|
| 1180 | ! ENDIF |
---|
| 1181 | ! IF ( zw(k-1) > circle_points(2,i_ip) ) THEN |
---|
| 1182 | ! dz_int = 0.0_wp |
---|
| 1183 | ! ENDIF |
---|
| 1184 | ! IF ( zw(k) < circle_points(1,i_ip) ) THEN |
---|
| 1185 | ! dz_int = 0.0_wp |
---|
| 1186 | ! ENDIF |
---|
| 1187 | ! IF ( ( nxlg <= i ) .AND. ( nxrg >= i ) ) THEN |
---|
[1912] | 1188 | ! nac_cd_surf(k,j,i) = nac_cd_surf(k,j,i) + & !<--- segmentation fault |
---|
| 1189 | ! dy_int * dz_int * turb_cd_nacelle(inot) |
---|
[1864] | 1190 | ! ENDIF |
---|
| 1191 | ! ENDDO |
---|
| 1192 | ! ENDIF |
---|
| 1193 | ! ENDDO |
---|
| 1194 | ! ENDIF |
---|
[1819] | 1195 | ! |
---|
[1864] | 1196 | ! ENDDO |
---|
[1819] | 1197 | ! |
---|
[1912] | 1198 | ! CALL exchange_horiz( nac_cd_surf, nbgp ) !<--- segmentation fault |
---|
[1819] | 1199 | |
---|
[1864] | 1200 | ENDDO ! end of loop over turbines |
---|
[1819] | 1201 | |
---|
[1912] | 1202 | tow_cd_surf = tow_cd_surf / ( dx * dy * dz ) ! Normalize tower drag |
---|
| 1203 | nac_cd_surf = nac_cd_surf / ( dx * dy * dz ) ! Normalize nacelle drag |
---|
[1819] | 1204 | |
---|
[1912] | 1205 | CALL wtm_read_blade_tables |
---|
[1839] | 1206 | |
---|
[1864] | 1207 | END SUBROUTINE wtm_init |
---|
| 1208 | |
---|
| 1209 | |
---|
[1819] | 1210 | !------------------------------------------------------------------------------! |
---|
[1864] | 1211 | ! Description: |
---|
| 1212 | ! ------------ |
---|
| 1213 | !> Read in layout of the rotor blade , the lift and drag tables |
---|
| 1214 | !> and the distribution of lift and drag tables along the blade |
---|
[1819] | 1215 | !------------------------------------------------------------------------------! |
---|
[1912] | 1216 | ! |
---|
| 1217 | SUBROUTINE wtm_read_blade_tables |
---|
[1819] | 1218 | |
---|
[1839] | 1219 | |
---|
[1864] | 1220 | IMPLICIT NONE |
---|
| 1221 | |
---|
| 1222 | INTEGER(iwp) :: ii !< running index |
---|
| 1223 | INTEGER(iwp) :: jj !< running index |
---|
[1843] | 1224 | |
---|
[1864] | 1225 | INTEGER(iwp) :: ierrn !< |
---|
| 1226 | |
---|
| 1227 | CHARACTER(200) :: chmess !< Read in string |
---|
[1839] | 1228 | |
---|
[1864] | 1229 | INTEGER(iwp) :: dlen !< no. rows of local table |
---|
| 1230 | INTEGER(iwp) :: dlenbl !< no. rows of cd, cl table |
---|
| 1231 | INTEGER(iwp) :: ialpha !< table position of current alpha value |
---|
| 1232 | INTEGER(iwp) :: iialpha !< |
---|
| 1233 | INTEGER(iwp) :: iir !< |
---|
| 1234 | INTEGER(iwp) :: radres !< radial resolution |
---|
| 1235 | INTEGER(iwp) :: t1 !< no. of airfoil |
---|
| 1236 | INTEGER(iwp) :: t2 !< no. of airfoil |
---|
| 1237 | INTEGER(iwp) :: trow !< |
---|
| 1238 | INTEGER(iwp) :: dlenbl_int !< no. rows of interpolated cd, cl tables |
---|
[1839] | 1239 | |
---|
[1864] | 1240 | REAL(wp) :: alpha_attack_i !< |
---|
| 1241 | REAL(wp) :: weight_a !< |
---|
| 1242 | REAL(wp) :: weight_b !< |
---|
[1839] | 1243 | |
---|
[1864] | 1244 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ttoint1 !< |
---|
| 1245 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ttoint2 !< |
---|
[1839] | 1246 | |
---|
[1864] | 1247 | REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cd_sel1 !< |
---|
| 1248 | REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cd_sel2 !< |
---|
| 1249 | REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cl_sel1 !< |
---|
| 1250 | REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cl_sel2 !< |
---|
| 1251 | REAL(wp), DIMENSION(:), ALLOCATABLE :: read_cl_cd !< read in var array |
---|
[1839] | 1252 | |
---|
[1864] | 1253 | REAL(wp), DIMENSION(:), ALLOCATABLE :: alpha_attack_tab !< |
---|
| 1254 | REAL(wp), DIMENSION(:), ALLOCATABLE :: trad1 !< |
---|
| 1255 | REAL(wp), DIMENSION(:), ALLOCATABLE :: trad2 !< |
---|
| 1256 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cd_table !< |
---|
| 1257 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cl_table !< |
---|
[1839] | 1258 | |
---|
[1864] | 1259 | ALLOCATE ( read_cl_cd(1:2*nairfoils+1) ) |
---|
[1839] | 1260 | |
---|
[1864] | 1261 | ! |
---|
| 1262 | !-- Read in the distribution of lift and drag tables along the blade, the |
---|
| 1263 | !-- layout of the rotor blade and the lift and drag tables: |
---|
| 1264 | |
---|
[2323] | 1265 | OPEN ( 90, FILE='WTM_DATA', STATUS='OLD', FORM='FORMATTED', IOSTAT=ierrn ) |
---|
[1864] | 1266 | |
---|
| 1267 | IF ( ierrn /= 0 ) THEN |
---|
| 1268 | message_string = 'file WTM_DATA does not exist' |
---|
[2322] | 1269 | CALL message( 'wtm_init', 'PA0464', 1, 2, 0, 6, 0 ) |
---|
[1864] | 1270 | ENDIF |
---|
| 1271 | ! |
---|
| 1272 | !-- Read distribution table: |
---|
| 1273 | |
---|
| 1274 | dlen = 0 |
---|
| 1275 | |
---|
[2323] | 1276 | READ ( 90, '(3/)' ) |
---|
[1864] | 1277 | |
---|
| 1278 | rloop3: DO |
---|
[2323] | 1279 | READ ( 90, *, IOSTAT=ierrn ) chmess |
---|
[1864] | 1280 | IF ( ierrn < 0 .OR. chmess == '#' .OR. chmess == '') EXIT rloop3 |
---|
| 1281 | dlen = dlen + 1 |
---|
| 1282 | ENDDO rloop3 |
---|
| 1283 | |
---|
| 1284 | ALLOCATE( trad1(1:dlen), trad2(1:dlen), ttoint1(1:dlen), ttoint2(1:dlen)) |
---|
| 1285 | |
---|
| 1286 | DO jj = 1,dlen+1 |
---|
[2323] | 1287 | BACKSPACE ( 90, IOSTAT=ierrn ) |
---|
[1864] | 1288 | ENDDO |
---|
| 1289 | |
---|
| 1290 | DO jj = 1,dlen |
---|
[2323] | 1291 | READ ( 90, * ) trad1(jj), trad2(jj), ttoint1(jj), ttoint2(jj) |
---|
[1864] | 1292 | ENDDO |
---|
| 1293 | |
---|
| 1294 | ! |
---|
| 1295 | !-- Read layout table: |
---|
| 1296 | |
---|
[1839] | 1297 | dlen = 0 |
---|
[1864] | 1298 | |
---|
[2323] | 1299 | READ ( 90, '(3/)') |
---|
[1864] | 1300 | |
---|
[1843] | 1301 | rloop1: DO |
---|
[2323] | 1302 | READ ( 90, *, IOSTAT=ierrn ) chmess |
---|
[1864] | 1303 | IF ( ierrn < 0 .OR. chmess == '#' .OR. chmess == '') EXIT rloop1 |
---|
| 1304 | dlen = dlen + 1 |
---|
[1843] | 1305 | ENDDO rloop1 |
---|
[1864] | 1306 | |
---|
[1839] | 1307 | ALLOCATE( lrd(1:dlen), ard(1:dlen), crd(1:dlen) ) |
---|
[1864] | 1308 | DO jj = 1, dlen+1 |
---|
[2323] | 1309 | BACKSPACE ( 90, IOSTAT=ierrn ) |
---|
[1864] | 1310 | ENDDO |
---|
| 1311 | DO jj = 1, dlen |
---|
[2323] | 1312 | READ ( 90, * ) lrd(jj), ard(jj), crd(jj) |
---|
[1839] | 1313 | ENDDO |
---|
[1819] | 1314 | |
---|
[1864] | 1315 | ! |
---|
| 1316 | !-- Read tables (turb_cl(alpha),turb_cd(alpha) for the different profiles: |
---|
[1819] | 1317 | |
---|
[1864] | 1318 | dlen = 0 |
---|
| 1319 | |
---|
[2323] | 1320 | READ ( 90, '(3/)' ) |
---|
[1864] | 1321 | |
---|
[1843] | 1322 | rloop2: DO |
---|
[2323] | 1323 | READ ( 90, *, IOSTAT=ierrn ) chmess |
---|
[1864] | 1324 | IF ( ierrn < 0 .OR. chmess == '#' .OR. chmess == '') EXIT rloop2 |
---|
| 1325 | dlen = dlen + 1 |
---|
| 1326 | ENDDO rloop2 |
---|
| 1327 | |
---|
[1912] | 1328 | ALLOCATE( alpha_attack_tab(1:dlen), turb_cl_table(1:dlen,1:nairfoils), & |
---|
[1864] | 1329 | turb_cd_table(1:dlen,1:nairfoils) ) |
---|
| 1330 | |
---|
| 1331 | DO jj = 1,dlen+1 |
---|
[2323] | 1332 | BACKSPACE ( 90, IOSTAT=ierrn ) |
---|
[1864] | 1333 | ENDDO |
---|
| 1334 | |
---|
[1839] | 1335 | DO jj = 1,dlen |
---|
[2323] | 1336 | READ ( 90, * ) read_cl_cd |
---|
[1864] | 1337 | alpha_attack_tab(jj) = read_cl_cd(1) |
---|
| 1338 | DO ii= 1, nairfoils |
---|
| 1339 | turb_cl_table(jj,ii) = read_cl_cd(ii*2) |
---|
| 1340 | turb_cd_table(jj,ii) = read_cl_cd(ii*2+1) |
---|
[1819] | 1341 | ENDDO |
---|
[1864] | 1342 | |
---|
[1839] | 1343 | ENDDO |
---|
[1864] | 1344 | |
---|
| 1345 | dlenbl = dlen |
---|
| 1346 | |
---|
[2323] | 1347 | CLOSE ( 90 ) |
---|
[1819] | 1348 | |
---|
[1864] | 1349 | ! |
---|
| 1350 | !-- For each possible radial position (resolution: 0.1 m --> 630 values) and |
---|
| 1351 | !-- each possible angle of attack (resolution: 0.01 degrees --> 36000 values!) |
---|
| 1352 | !-- determine the lift and drag coefficient by interpolating between the |
---|
| 1353 | !-- tabulated values of each table (interpolate to current angle of attack) |
---|
| 1354 | !-- and between the tables (interpolate to current radial position): |
---|
[1839] | 1355 | |
---|
[1864] | 1356 | ALLOCATE( turb_cl_sel1(0:dlenbl) ) |
---|
| 1357 | ALLOCATE( turb_cl_sel2(0:dlenbl) ) |
---|
| 1358 | ALLOCATE( turb_cd_sel1(0:dlenbl) ) |
---|
| 1359 | ALLOCATE( turb_cd_sel2(0:dlenbl) ) |
---|
[1819] | 1360 | |
---|
[1864] | 1361 | radres = INT( rr(1) * 10.0_wp ) + 1_iwp |
---|
| 1362 | dlenbl_int = INT( 360.0_wp / accu_cl_cd_tab ) + 1_iwp |
---|
[1839] | 1363 | |
---|
[1864] | 1364 | |
---|
| 1365 | ALLOCATE( turb_cl_tab(0:dlenbl_int,1:radres) ) |
---|
| 1366 | ALLOCATE( turb_cd_tab(0:dlenbl_int,1:radres) ) |
---|
| 1367 | |
---|
| 1368 | |
---|
| 1369 | DO iir = 1, radres ! loop over radius |
---|
| 1370 | |
---|
| 1371 | DO iialpha = 1, dlenbl_int ! loop over angles |
---|
| 1372 | |
---|
| 1373 | cur_r = ( iir - 1_iwp ) * 0.1_wp |
---|
| 1374 | alpha_attack_i = -180.0_wp + REAL( iialpha-1 ) * accu_cl_cd_tab |
---|
| 1375 | ialpha = 1 |
---|
| 1376 | |
---|
[2152] | 1377 | DO WHILE ( ( alpha_attack_i > alpha_attack_tab(ialpha) ) .AND. (ialpha <= dlen ) ) |
---|
[1864] | 1378 | ialpha = ialpha + 1 |
---|
| 1379 | ENDDO |
---|
[1819] | 1380 | ! |
---|
[1864] | 1381 | !-- Find position in table |
---|
| 1382 | lct = MINLOC( ABS( trad1 - cur_r ) ) |
---|
| 1383 | ! lct(1) = lct(1) |
---|
[1819] | 1384 | |
---|
[1864] | 1385 | IF ( ( trad1(lct(1)) - cur_r ) .GT. 0.0 ) THEN |
---|
| 1386 | lct(1) = lct(1) - 1 |
---|
| 1387 | ENDIF |
---|
[1819] | 1388 | |
---|
[1864] | 1389 | trow = lct(1) |
---|
| 1390 | ! |
---|
| 1391 | !-- Calculate weights for interpolation |
---|
| 1392 | weight_a = ( trad2(trow) - cur_r ) / ( trad2(trow) - trad1(trow) ) |
---|
| 1393 | weight_b = ( cur_r - trad1(trow) ) / ( trad2(trow) - trad1(trow) ) |
---|
| 1394 | t1 = ttoint1(trow) |
---|
| 1395 | t2 = ttoint2(trow) |
---|
[1819] | 1396 | |
---|
[1912] | 1397 | IF ( t1 .EQ. t2 ) THEN ! if both are the same, the weights are NaN |
---|
| 1398 | weight_a = 0.5_wp ! then do interpolate in between same twice |
---|
| 1399 | weight_b = 0.5_wp ! using 0.5 as weight |
---|
[1864] | 1400 | ENDIF |
---|
[1819] | 1401 | |
---|
[1864] | 1402 | IF ( t1 == 0 .AND. t2 == 0 ) THEN |
---|
| 1403 | turb_cd_sel1 = 0.0_wp |
---|
| 1404 | turb_cd_sel2 = 0.0_wp |
---|
| 1405 | turb_cl_sel1 = 0.0_wp |
---|
| 1406 | turb_cl_sel2 = 0.0_wp |
---|
| 1407 | ELSE |
---|
| 1408 | turb_cd_sel1 = turb_cd_table(:,t1) |
---|
| 1409 | turb_cd_sel2 = turb_cd_table(:,t2) |
---|
| 1410 | turb_cl_sel1 = turb_cl_table(:,t1) |
---|
| 1411 | turb_cl_sel2 = turb_cl_table(:,t2) |
---|
| 1412 | ENDIF |
---|
| 1413 | |
---|
[1912] | 1414 | ! |
---|
| 1415 | !-- Interpolation of lift and drag coefficiencts on fine grid of radius |
---|
| 1416 | !-- segments and angles of attack |
---|
| 1417 | |
---|
| 1418 | turb_cl_tab(iialpha,iir) = ( alpha_attack_tab(ialpha) - & |
---|
| 1419 | alpha_attack_i ) / & |
---|
[1864] | 1420 | ( alpha_attack_tab(ialpha) - & |
---|
| 1421 | alpha_attack_tab(ialpha-1) ) * & |
---|
| 1422 | ( weight_a * turb_cl_sel1(ialpha-1) + & |
---|
| 1423 | weight_b * turb_cl_sel2(ialpha-1) ) +& |
---|
| 1424 | ( alpha_attack_i - & |
---|
| 1425 | alpha_attack_tab(ialpha-1) ) / & |
---|
| 1426 | ( alpha_attack_tab(ialpha) - & |
---|
| 1427 | alpha_attack_tab(ialpha-1) ) * & |
---|
| 1428 | ( weight_a * turb_cl_sel1(ialpha) + & |
---|
| 1429 | weight_b * turb_cl_sel2(ialpha) ) |
---|
[1912] | 1430 | turb_cd_tab(iialpha,iir) = ( alpha_attack_tab(ialpha) - & |
---|
| 1431 | alpha_attack_i ) / & |
---|
[1864] | 1432 | ( alpha_attack_tab(ialpha) - & |
---|
| 1433 | alpha_attack_tab(ialpha-1) ) * & |
---|
| 1434 | ( weight_a * turb_cd_sel1(ialpha-1) + & |
---|
| 1435 | weight_b * turb_cd_sel2(ialpha-1) ) +& |
---|
| 1436 | ( alpha_attack_i - & |
---|
| 1437 | alpha_attack_tab(ialpha-1) ) / & |
---|
| 1438 | ( alpha_attack_tab(ialpha) - & |
---|
| 1439 | alpha_attack_tab(ialpha-1) ) * & |
---|
| 1440 | ( weight_a * turb_cd_sel1(ialpha) + & |
---|
| 1441 | weight_b * turb_cd_sel2(ialpha) ) |
---|
[1819] | 1442 | |
---|
[1912] | 1443 | ENDDO ! end loop over angles of attack |
---|
[1819] | 1444 | |
---|
[1864] | 1445 | ENDDO ! end loop over radius |
---|
| 1446 | |
---|
[1912] | 1447 | END SUBROUTINE wtm_read_blade_tables |
---|
[1819] | 1448 | |
---|
| 1449 | |
---|
[1864] | 1450 | !------------------------------------------------------------------------------! |
---|
| 1451 | ! Description: |
---|
| 1452 | ! ------------ |
---|
| 1453 | !> The projection matrix for the coordinate system of therotor disc in respect |
---|
| 1454 | !> to the yaw and tilt angle of the rotor is calculated |
---|
| 1455 | !------------------------------------------------------------------------------! |
---|
| 1456 | SUBROUTINE wtm_rotate_rotor( inot ) |
---|
[1819] | 1457 | |
---|
[1864] | 1458 | |
---|
| 1459 | IMPLICIT NONE |
---|
| 1460 | |
---|
| 1461 | INTEGER(iwp) :: inot |
---|
| 1462 | ! |
---|
| 1463 | !-- Calculation of the rotation matrix for the application of the tilt to |
---|
| 1464 | !-- the rotors |
---|
| 1465 | rot_eigen_rad(1) = SIN( phi_yaw(inot) ) ! x-component of the radial eigenvector |
---|
| 1466 | rot_eigen_rad(2) = COS( phi_yaw(inot) ) ! y-component of the radial eigenvector |
---|
| 1467 | rot_eigen_rad(3) = 0.0_wp ! z-component of the radial eigenvector |
---|
| 1468 | |
---|
| 1469 | rot_eigen_azi(1) = 0.0_wp ! x-component of the azimuth eigenvector |
---|
| 1470 | rot_eigen_azi(2) = 0.0_wp ! y-component of the azimuth eigenvector |
---|
| 1471 | rot_eigen_azi(3) = 1.0_wp ! z-component of the azimuth eigenvector |
---|
| 1472 | |
---|
| 1473 | rot_eigen_nor(1) = COS( phi_yaw(inot) ) ! x-component of the normal eigenvector |
---|
| 1474 | rot_eigen_nor(2) = -SIN( phi_yaw(inot) ) ! y-component of the normal eigenvector |
---|
| 1475 | rot_eigen_nor(3) = 0.0_wp ! z-component of the normal eigenvector |
---|
[1839] | 1476 | |
---|
[1864] | 1477 | ! |
---|
| 1478 | !-- Calculation of the coordinate transformation matrix to apply a tilt to |
---|
| 1479 | !-- the rotor. If tilt = 0, rot_coord_trans is a unit matrix. |
---|
[1819] | 1480 | |
---|
[1912] | 1481 | rot_coord_trans(inot,1,1) = rot_eigen_rad(1)**2 * & |
---|
[1864] | 1482 | ( 1.0_wp - COS( tilt ) ) + COS( tilt ) |
---|
[1912] | 1483 | rot_coord_trans(inot,1,2) = rot_eigen_rad(1) * rot_eigen_rad(2) * & |
---|
| 1484 | ( 1.0_wp - COS( tilt ) ) - & |
---|
[1864] | 1485 | rot_eigen_rad(3) * SIN( tilt ) |
---|
[1912] | 1486 | rot_coord_trans(inot,1,3) = rot_eigen_rad(1) * rot_eigen_rad(3) * & |
---|
| 1487 | ( 1.0_wp - COS( tilt ) ) + & |
---|
[1864] | 1488 | rot_eigen_rad(2) * SIN( tilt ) |
---|
[1912] | 1489 | rot_coord_trans(inot,2,1) = rot_eigen_rad(2) * rot_eigen_rad(1) * & |
---|
| 1490 | ( 1.0_wp - COS( tilt ) ) + & |
---|
[1864] | 1491 | rot_eigen_rad(3) * SIN( tilt ) |
---|
[1912] | 1492 | rot_coord_trans(inot,2,2) = rot_eigen_rad(2)**2 * & |
---|
[1864] | 1493 | ( 1.0_wp - COS( tilt ) ) + COS( tilt ) |
---|
[1912] | 1494 | rot_coord_trans(inot,2,3) = rot_eigen_rad(2) * rot_eigen_rad(3) * & |
---|
| 1495 | ( 1.0_wp - COS( tilt ) ) - & |
---|
[1864] | 1496 | rot_eigen_rad(1) * SIN( tilt ) |
---|
[1912] | 1497 | rot_coord_trans(inot,3,1) = rot_eigen_rad(3) * rot_eigen_rad(1) * & |
---|
| 1498 | ( 1.0_wp - COS( tilt ) ) - & |
---|
[1864] | 1499 | rot_eigen_rad(2) * SIN( tilt ) |
---|
[1912] | 1500 | rot_coord_trans(inot,3,2) = rot_eigen_rad(3) * rot_eigen_rad(2) * & |
---|
| 1501 | ( 1.0_wp - COS( tilt ) ) + & |
---|
[1864] | 1502 | rot_eigen_rad(1) * SIN( tilt ) |
---|
[1912] | 1503 | rot_coord_trans(inot,3,3) = rot_eigen_rad(3)**2 * & |
---|
[1864] | 1504 | ( 1.0_wp - COS( tilt ) ) + COS( tilt ) |
---|
[1839] | 1505 | |
---|
[1864] | 1506 | ! |
---|
| 1507 | !-- Vectors for the Transformation of forces from the rotor's spheric |
---|
| 1508 | !-- coordinate system to the cartesian coordinate system |
---|
| 1509 | rotx(inot,:) = MATMUL( rot_coord_trans(inot,:,:), rot_eigen_nor ) |
---|
| 1510 | roty(inot,:) = MATMUL( rot_coord_trans(inot,:,:), rot_eigen_rad ) |
---|
| 1511 | rotz(inot,:) = MATMUL( rot_coord_trans(inot,:,:), rot_eigen_azi ) |
---|
| 1512 | |
---|
| 1513 | END SUBROUTINE wtm_rotate_rotor |
---|
[1839] | 1514 | |
---|
| 1515 | |
---|
[1819] | 1516 | !------------------------------------------------------------------------------! |
---|
| 1517 | ! Description: |
---|
| 1518 | ! ------------ |
---|
[1839] | 1519 | !> Calculation of the forces generated by the wind turbine |
---|
[1819] | 1520 | !------------------------------------------------------------------------------! |
---|
| 1521 | SUBROUTINE wtm_forces |
---|
| 1522 | |
---|
[1864] | 1523 | |
---|
[1819] | 1524 | IMPLICIT NONE |
---|
| 1525 | |
---|
| 1526 | CHARACTER (LEN=2) :: turbine_id |
---|
| 1527 | |
---|
[1839] | 1528 | INTEGER(iwp) :: i, j, k !< loop indices |
---|
| 1529 | INTEGER(iwp) :: inot !< turbine loop index (turbine id) |
---|
| 1530 | INTEGER(iwp) :: iialpha, iir !< |
---|
| 1531 | INTEGER(iwp) :: rseg, rseg_int !< |
---|
| 1532 | INTEGER(iwp) :: ring, ring_int !< |
---|
| 1533 | INTEGER(iwp) :: ii, jj, kk !< |
---|
[1912] | 1534 | |
---|
[2232] | 1535 | REAL(wp) :: flag !< flag to mask topography grid points |
---|
[1839] | 1536 | REAL(wp) :: sin_rot, cos_rot !< |
---|
| 1537 | REAL(wp) :: sin_yaw, cos_yaw !< |
---|
[1912] | 1538 | |
---|
| 1539 | REAL(wp) :: aa, bb, cc, dd !< interpolation distances |
---|
| 1540 | REAL(wp) :: gg !< interpolation volume var |
---|
| 1541 | |
---|
| 1542 | REAL(wp) :: dist_u_3d, dist_v_3d, dist_w_3d !< smearing distances |
---|
| 1543 | |
---|
| 1544 | |
---|
[1839] | 1545 | ! |
---|
[1819] | 1546 | ! Variables for pitch control |
---|
| 1547 | REAL(wp) :: torque_max=0.0_wp |
---|
[1839] | 1548 | LOGICAL :: pitch_sw=.FALSE. |
---|
| 1549 | |
---|
[1819] | 1550 | INTEGER(iwp), DIMENSION(1) :: lct=0 |
---|
[1839] | 1551 | REAL(wp), DIMENSION(1) :: rad_d=0.0_wp |
---|
[2152] | 1552 | |
---|
| 1553 | REAL(wp) :: tl_factor !< factor for tip loss correction |
---|
[1819] | 1554 | |
---|
| 1555 | |
---|
[1864] | 1556 | CALL cpu_log( log_point_s(61), 'wtm_forces', 'start' ) |
---|
[1819] | 1557 | |
---|
| 1558 | |
---|
| 1559 | IF ( simulated_time >= time_turbine_on ) THEN |
---|
| 1560 | |
---|
[1864] | 1561 | ! |
---|
[1819] | 1562 | !-- Set forces to zero for each new time step: |
---|
| 1563 | thrust(:,:,:) = 0.0_wp |
---|
| 1564 | torque_y(:,:,:) = 0.0_wp |
---|
| 1565 | torque_z(:,:,:) = 0.0_wp |
---|
| 1566 | torque_total(:) = 0.0_wp |
---|
[1912] | 1567 | rot_tend_x(:,:,:) = 0.0_wp |
---|
| 1568 | rot_tend_y(:,:,:) = 0.0_wp |
---|
| 1569 | rot_tend_z(:,:,:) = 0.0_wp |
---|
| 1570 | thrust_rotor(:) = 0.0_wp |
---|
[1819] | 1571 | ! |
---|
| 1572 | !-- Loop over number of turbines: |
---|
| 1573 | DO inot = 1, nturbines |
---|
| 1574 | |
---|
| 1575 | cos_yaw = COS(phi_yaw(inot)) |
---|
| 1576 | sin_yaw = SIN(phi_yaw(inot)) |
---|
| 1577 | ! |
---|
[1839] | 1578 | !-- Loop over rings of each turbine: |
---|
[1819] | 1579 | DO ring = 1, nrings(inot) |
---|
| 1580 | |
---|
| 1581 | thrust_seg(:) = 0.0_wp |
---|
| 1582 | torque_seg_y(:) = 0.0_wp |
---|
| 1583 | torque_seg_z(:) = 0.0_wp |
---|
| 1584 | ! |
---|
| 1585 | !-- Determine distance between each ring (center) and the hub: |
---|
| 1586 | cur_r = (ring - 0.5_wp) * delta_r(inot) |
---|
| 1587 | |
---|
| 1588 | ! |
---|
[1839] | 1589 | !-- Loop over segments of each ring of each turbine: |
---|
[1819] | 1590 | DO rseg = 1, nsegs(ring,inot) |
---|
| 1591 | ! |
---|
[1864] | 1592 | !-- !-----------------------------------------------------------! |
---|
| 1593 | !-- !-- Determine coordinates of the ring segments --! |
---|
| 1594 | !-- !-----------------------------------------------------------! |
---|
[1819] | 1595 | ! |
---|
[1864] | 1596 | !-- Determine angle of ring segment towards zero degree angle of |
---|
| 1597 | !-- rotor system (at zero degree rotor direction vectors aligned |
---|
| 1598 | !-- with y-axis): |
---|
[1819] | 1599 | phi_rotor = rseg * 2.0_wp * pi / nsegs(ring,inot) |
---|
| 1600 | cos_rot = COS( phi_rotor ) |
---|
| 1601 | sin_rot = SIN( phi_rotor ) |
---|
| 1602 | |
---|
[1864] | 1603 | !-- Now the direction vectors can be determined with respect to |
---|
| 1604 | !-- the yaw and tilt angle: |
---|
[1819] | 1605 | re(1) = cos_rot * sin_yaw |
---|
[1839] | 1606 | re(2) = cos_rot * cos_yaw |
---|
[1819] | 1607 | re(3) = sin_rot |
---|
| 1608 | |
---|
| 1609 | rote = MATMUL( rot_coord_trans(inot,:,:), re ) |
---|
| 1610 | ! |
---|
| 1611 | !-- Coordinates of the single segments (center points): |
---|
| 1612 | rbx(ring,rseg) = rcx(inot) + cur_r * rote(1) |
---|
| 1613 | rby(ring,rseg) = rcy(inot) + cur_r * rote(2) |
---|
| 1614 | rbz(ring,rseg) = rcz(inot) + cur_r * rote(3) |
---|
| 1615 | |
---|
[1864] | 1616 | !-- !-----------------------------------------------------------! |
---|
| 1617 | !-- !-- Interpolation of the velocity components from the --! |
---|
| 1618 | !-- !-- cartesian grid point to the coordinates of each ring --! |
---|
| 1619 | !-- !-- segment (follows a method used in the particle model) --! |
---|
| 1620 | !-- !-----------------------------------------------------------! |
---|
[1819] | 1621 | |
---|
| 1622 | u_int(inot,ring,rseg) = 0.0_wp |
---|
| 1623 | u_int_1_l(inot,ring,rseg) = 0.0_wp |
---|
| 1624 | |
---|
| 1625 | v_int(inot,ring,rseg) = 0.0_wp |
---|
| 1626 | v_int_1_l(inot,ring,rseg) = 0.0_wp |
---|
| 1627 | |
---|
| 1628 | w_int(inot,ring,rseg) = 0.0_wp |
---|
| 1629 | w_int_1_l(inot,ring,rseg) = 0.0_wp |
---|
| 1630 | |
---|
| 1631 | ! |
---|
| 1632 | !-- Interpolation of the u-component: |
---|
| 1633 | |
---|
| 1634 | ii = rbx(ring,rseg) * ddx |
---|
| 1635 | jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy |
---|
[2553] | 1636 | kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz |
---|
[1819] | 1637 | ! |
---|
[1864] | 1638 | !-- Interpolate only if all required information is available on |
---|
| 1639 | !-- the current PE: |
---|
[1839] | 1640 | IF ( ( ii >= nxl ) .AND. ( ii <= nxr ) ) THEN |
---|
| 1641 | IF ( ( jj >= nys ) .AND. ( jj <= nyn ) ) THEN |
---|
[1819] | 1642 | |
---|
[1839] | 1643 | aa = ( ( ii + 1 ) * dx - rbx(ring,rseg) ) * & |
---|
| 1644 | ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) |
---|
| 1645 | bb = ( rbx(ring,rseg) - ii * dx ) * & |
---|
| 1646 | ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) |
---|
| 1647 | cc = ( ( ii+1 ) * dx - rbx(ring,rseg) ) * & |
---|
| 1648 | ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) |
---|
| 1649 | dd = ( rbx(ring,rseg) - ii * dx ) * & |
---|
| 1650 | ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) |
---|
[1819] | 1651 | gg = dx * dy |
---|
| 1652 | |
---|
[1864] | 1653 | u_int_l = ( aa * u(kk,jj,ii) + & |
---|
| 1654 | bb * u(kk,jj,ii+1) + & |
---|
| 1655 | cc * u(kk,jj+1,ii) + & |
---|
| 1656 | dd * u(kk,jj+1,ii+1) & |
---|
[1819] | 1657 | ) / gg |
---|
| 1658 | |
---|
[1864] | 1659 | u_int_u = ( aa * u(kk+1,jj,ii) + & |
---|
| 1660 | bb * u(kk+1,jj,ii+1) + & |
---|
| 1661 | cc * u(kk+1,jj+1,ii) + & |
---|
| 1662 | dd * u(kk+1,jj+1,ii+1) & |
---|
[1819] | 1663 | ) / gg |
---|
| 1664 | |
---|
[1864] | 1665 | u_int_1_l(inot,ring,rseg) = u_int_l + & |
---|
| 1666 | ( rbz(ring,rseg) - zu(kk) ) / dz * & |
---|
[1819] | 1667 | ( u_int_u - u_int_l ) |
---|
| 1668 | |
---|
| 1669 | ELSE |
---|
| 1670 | u_int_1_l(inot,ring,rseg) = 0.0_wp |
---|
| 1671 | ENDIF |
---|
| 1672 | ELSE |
---|
| 1673 | u_int_1_l(inot,ring,rseg) = 0.0_wp |
---|
| 1674 | ENDIF |
---|
| 1675 | |
---|
| 1676 | |
---|
| 1677 | ! |
---|
| 1678 | !-- Interpolation of the v-component: |
---|
| 1679 | ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx |
---|
[1839] | 1680 | jj = rby(ring,rseg) * ddy |
---|
[1819] | 1681 | kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz |
---|
| 1682 | ! |
---|
[1864] | 1683 | !-- Interpolate only if all required information is available on |
---|
| 1684 | !-- the current PE: |
---|
[1839] | 1685 | IF ( ( ii >= nxl ) .AND. ( ii <= nxr ) ) THEN |
---|
| 1686 | IF ( ( jj >= nys ) .AND. ( jj <= nyn ) ) THEN |
---|
[1819] | 1687 | |
---|
[1839] | 1688 | aa = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & |
---|
| 1689 | ( ( jj + 1 ) * dy - rby(ring,rseg) ) |
---|
| 1690 | bb = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & |
---|
| 1691 | ( ( jj + 1 ) * dy - rby(ring,rseg) ) |
---|
| 1692 | cc = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & |
---|
| 1693 | ( rby(ring,rseg) - jj * dy ) |
---|
| 1694 | dd = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & |
---|
| 1695 | ( rby(ring,rseg) - jj * dy ) |
---|
[1819] | 1696 | gg = dx * dy |
---|
| 1697 | |
---|
[1864] | 1698 | v_int_l = ( aa * v(kk,jj,ii) + & |
---|
| 1699 | bb * v(kk,jj,ii+1) + & |
---|
| 1700 | cc * v(kk,jj+1,ii) + & |
---|
| 1701 | dd * v(kk,jj+1,ii+1) & |
---|
[1819] | 1702 | ) / gg |
---|
| 1703 | |
---|
[1864] | 1704 | v_int_u = ( aa * v(kk+1,jj,ii) + & |
---|
| 1705 | bb * v(kk+1,jj,ii+1) + & |
---|
| 1706 | cc * v(kk+1,jj+1,ii) + & |
---|
| 1707 | dd * v(kk+1,jj+1,ii+1) & |
---|
[1819] | 1708 | ) / gg |
---|
| 1709 | |
---|
[1864] | 1710 | v_int_1_l(inot,ring,rseg) = v_int_l + & |
---|
| 1711 | ( rbz(ring,rseg) - zu(kk) ) / dz * & |
---|
[1819] | 1712 | ( v_int_u - v_int_l ) |
---|
| 1713 | |
---|
| 1714 | ELSE |
---|
| 1715 | v_int_1_l(inot,ring,rseg) = 0.0_wp |
---|
| 1716 | ENDIF |
---|
| 1717 | ELSE |
---|
| 1718 | v_int_1_l(inot,ring,rseg) = 0.0_wp |
---|
| 1719 | ENDIF |
---|
| 1720 | |
---|
| 1721 | |
---|
| 1722 | ! |
---|
| 1723 | !-- Interpolation of the w-component: |
---|
| 1724 | ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx |
---|
| 1725 | jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy |
---|
[1839] | 1726 | kk = rbz(ring,rseg) / dz |
---|
[1819] | 1727 | ! |
---|
[1864] | 1728 | !-- Interpolate only if all required information is available on |
---|
| 1729 | !-- the current PE: |
---|
[1839] | 1730 | IF ( ( ii >= nxl ) .AND. ( ii <= nxr ) ) THEN |
---|
| 1731 | IF ( ( jj >= nys ) .AND. ( jj <= nyn ) ) THEN |
---|
[1819] | 1732 | |
---|
[1839] | 1733 | aa = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & |
---|
| 1734 | ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) |
---|
| 1735 | bb = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & |
---|
| 1736 | ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) |
---|
| 1737 | cc = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & |
---|
| 1738 | ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) |
---|
| 1739 | dd = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & |
---|
| 1740 | ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) |
---|
[1819] | 1741 | gg = dx * dy |
---|
| 1742 | |
---|
[1864] | 1743 | w_int_l = ( aa * w(kk,jj,ii) + & |
---|
| 1744 | bb * w(kk,jj,ii+1) + & |
---|
| 1745 | cc * w(kk,jj+1,ii) + & |
---|
| 1746 | dd * w(kk,jj+1,ii+1) & |
---|
[1819] | 1747 | ) / gg |
---|
| 1748 | |
---|
[1864] | 1749 | w_int_u = ( aa * w(kk+1,jj,ii) + & |
---|
| 1750 | bb * w(kk+1,jj,ii+1) + & |
---|
| 1751 | cc * w(kk+1,jj+1,ii) + & |
---|
| 1752 | dd * w(kk+1,jj+1,ii+1) & |
---|
[1819] | 1753 | ) / gg |
---|
| 1754 | |
---|
[1864] | 1755 | w_int_1_l(inot,ring,rseg) = w_int_l + & |
---|
| 1756 | ( rbz(ring,rseg) - zw(kk) ) / dz * & |
---|
[1819] | 1757 | ( w_int_u - w_int_l ) |
---|
| 1758 | ELSE |
---|
| 1759 | w_int_1_l(inot,ring,rseg) = 0.0_wp |
---|
| 1760 | ENDIF |
---|
| 1761 | ELSE |
---|
| 1762 | w_int_1_l(inot,ring,rseg) = 0.0_wp |
---|
| 1763 | ENDIF |
---|
| 1764 | |
---|
| 1765 | ENDDO |
---|
| 1766 | ENDDO |
---|
| 1767 | |
---|
| 1768 | ENDDO |
---|
| 1769 | |
---|
| 1770 | ! |
---|
| 1771 | !-- Exchange between PEs (information required on each PE): |
---|
| 1772 | #if defined( __parallel ) |
---|
[1839] | 1773 | CALL MPI_ALLREDUCE( u_int_1_l, u_int, nturbines * MAXVAL(nrings) * & |
---|
| 1774 | MAXVAL(nsegs), MPI_REAL, MPI_SUM, comm2d, ierr ) |
---|
| 1775 | CALL MPI_ALLREDUCE( v_int_1_l, v_int, nturbines * MAXVAL(nrings) * & |
---|
| 1776 | MAXVAL(nsegs), MPI_REAL, MPI_SUM, comm2d, ierr ) |
---|
| 1777 | CALL MPI_ALLREDUCE( w_int_1_l, w_int, nturbines * MAXVAL(nrings) * & |
---|
| 1778 | MAXVAL(nsegs), MPI_REAL, MPI_SUM, comm2d, ierr ) |
---|
[1819] | 1779 | #else |
---|
| 1780 | u_int = u_int_1_l |
---|
| 1781 | v_int = v_int_1_l |
---|
| 1782 | w_int = w_int_1_l |
---|
| 1783 | #endif |
---|
| 1784 | |
---|
| 1785 | |
---|
| 1786 | ! |
---|
| 1787 | !-- Loop over number of turbines: |
---|
[1912] | 1788 | |
---|
[1819] | 1789 | DO inot = 1, nturbines |
---|
[1912] | 1790 | pit_loop: DO |
---|
[1819] | 1791 | |
---|
[1912] | 1792 | IF ( pitch_sw ) THEN |
---|
[1839] | 1793 | torque_total(inot) = 0.0_wp |
---|
[1912] | 1794 | thrust_rotor(inot) = 0.0_wp |
---|
| 1795 | pitch_add(inot) = pitch_add(inot) + 0.25_wp |
---|
| 1796 | ! IF ( myid == 0 ) PRINT*, 'Pitch', inot, pitch_add(inot) |
---|
| 1797 | ELSE |
---|
| 1798 | cos_yaw = COS(phi_yaw(inot)) |
---|
| 1799 | sin_yaw = SIN(phi_yaw(inot)) |
---|
| 1800 | IF ( pitch_control ) THEN |
---|
| 1801 | pitch_add(inot) = MAX(pitch_add_old(inot) - pitch_rate * & |
---|
| 1802 | dt_3d , 0.0_wp ) |
---|
| 1803 | ENDIF |
---|
[1819] | 1804 | ENDIF |
---|
| 1805 | |
---|
[1839] | 1806 | ! |
---|
| 1807 | !-- Loop over rings of each turbine: |
---|
[1819] | 1808 | DO ring = 1, nrings(inot) |
---|
| 1809 | ! |
---|
| 1810 | !-- Determine distance between each ring (center) and the hub: |
---|
| 1811 | cur_r = (ring - 0.5_wp) * delta_r(inot) |
---|
[1839] | 1812 | ! |
---|
| 1813 | !-- Loop over segments of each ring of each turbine: |
---|
[1819] | 1814 | DO rseg = 1, nsegs(ring,inot) |
---|
| 1815 | ! |
---|
[1839] | 1816 | !-- Determine angle of ring segment towards zero degree angle of |
---|
| 1817 | !-- rotor system (at zero degree rotor direction vectors aligned |
---|
| 1818 | !-- with y-axis): |
---|
[1819] | 1819 | phi_rotor = rseg * 2.0_wp * pi / nsegs(ring,inot) |
---|
| 1820 | cos_rot = COS(phi_rotor) |
---|
| 1821 | sin_rot = SIN(phi_rotor) |
---|
| 1822 | ! |
---|
[1839] | 1823 | !-- Now the direction vectors can be determined with respect to |
---|
| 1824 | !-- the yaw and tilt angle: |
---|
[1819] | 1825 | re(1) = cos_rot * sin_yaw |
---|
| 1826 | re(2) = cos_rot * cos_yaw |
---|
| 1827 | re(3) = sin_rot |
---|
| 1828 | |
---|
[1864] | 1829 | ! The current unit vector in azimuthal direction: |
---|
[1819] | 1830 | rea(1) = - sin_rot * sin_yaw |
---|
| 1831 | rea(2) = - sin_rot * cos_yaw |
---|
| 1832 | rea(3) = cos_rot |
---|
| 1833 | |
---|
| 1834 | ! |
---|
[1864] | 1835 | !-- To respect the yawing angle for the calculations of |
---|
| 1836 | !-- velocities and forces the unit vectors perpendicular to the |
---|
| 1837 | !-- rotor area in direction of the positive yaw angle are defined: |
---|
[1819] | 1838 | ren(1) = cos_yaw |
---|
| 1839 | ren(2) = - sin_yaw |
---|
| 1840 | ren(3) = 0.0_wp |
---|
| 1841 | ! |
---|
| 1842 | !-- Multiplication with the coordinate transformation matrix |
---|
[1864] | 1843 | !-- gives the final unit vector with consideration of the rotor |
---|
| 1844 | !-- tilt: |
---|
[1819] | 1845 | rote = MATMUL( rot_coord_trans(inot,:,:), re ) |
---|
| 1846 | rota = MATMUL( rot_coord_trans(inot,:,:), rea ) |
---|
| 1847 | rotn = MATMUL( rot_coord_trans(inot,:,:), ren ) |
---|
| 1848 | ! |
---|
| 1849 | !-- Coordinates of the single segments (center points): |
---|
| 1850 | rbx(ring,rseg) = rcx(inot) + cur_r * rote(1) |
---|
| 1851 | |
---|
| 1852 | rby(ring,rseg) = rcy(inot) + cur_r * rote(2) |
---|
| 1853 | |
---|
| 1854 | rbz(ring,rseg) = rcz(inot) + cur_r * rote(3) |
---|
| 1855 | |
---|
| 1856 | ! |
---|
[1864] | 1857 | !-- !-----------------------------------------------------------! |
---|
| 1858 | !-- !-- Calculation of various angles and relative velocities --! |
---|
| 1859 | !-- !-----------------------------------------------------------! |
---|
[1819] | 1860 | ! |
---|
[1864] | 1861 | !-- In the following the 3D-velocity field is projected its |
---|
[2553] | 1862 | !-- components perpendicular and parallel to the rotor area |
---|
[1819] | 1863 | !-- The calculation of forces will be done in the rotor- |
---|
| 1864 | !-- coordinates y' and z. |
---|
| 1865 | !-- The yaw angle will be reintroduced when the force is applied |
---|
| 1866 | !-- on the hydrodynamic equations |
---|
[1864] | 1867 | ! |
---|
| 1868 | !-- Projection of the xy-velocities relative to the rotor area |
---|
| 1869 | ! |
---|
[1819] | 1870 | !-- Velocity perpendicular to the rotor area: |
---|
[1864] | 1871 | u_rot = u_int(inot,ring,rseg)*rotn(1) + & |
---|
| 1872 | v_int(inot,ring,rseg)*rotn(2) + & |
---|
[1819] | 1873 | w_int(inot,ring,rseg)*rotn(3) |
---|
| 1874 | ! |
---|
[1864] | 1875 | !-- Projection of the 3D-velocity vector in the azimuthal |
---|
| 1876 | !-- direction: |
---|
| 1877 | vtheta(rseg) = rota(1) * u_int(inot,ring,rseg) + & |
---|
| 1878 | rota(2) * v_int(inot,ring,rseg) + & |
---|
[1819] | 1879 | rota(3) * w_int(inot,ring,rseg) |
---|
| 1880 | ! |
---|
[1864] | 1881 | !-- Determination of the angle phi_rel between the rotor plane |
---|
| 1882 | !-- and the direction of the flow relative to the rotor: |
---|
[1819] | 1883 | |
---|
[1864] | 1884 | phi_rel(rseg) = ATAN( u_rot / & |
---|
| 1885 | ( omega_rot(inot) * cur_r - & |
---|
[1819] | 1886 | vtheta(rseg) ) ) |
---|
| 1887 | |
---|
| 1888 | ! |
---|
[1864] | 1889 | !-- Interpolation of the local pitch angle from tabulated values |
---|
| 1890 | !-- to the current radial position: |
---|
[1819] | 1891 | |
---|
| 1892 | lct=minloc(ABS(cur_r-lrd)) |
---|
| 1893 | rad_d=cur_r-lrd(lct) |
---|
| 1894 | |
---|
| 1895 | IF (cur_r == 0.0_wp) THEN |
---|
| 1896 | alpha_attack(rseg) = 0.0_wp |
---|
| 1897 | ELSE IF (cur_r >= lrd(size(ard))) THEN |
---|
[1864] | 1898 | alpha_attack(rseg) = ( ard(size(ard)) + & |
---|
| 1899 | ard(size(ard)-1) ) / 2.0_wp |
---|
[1819] | 1900 | ELSE |
---|
[1864] | 1901 | alpha_attack(rseg) = ( ard(lct(1)) * & |
---|
| 1902 | ( ( lrd(lct(1)+1) - cur_r ) / & |
---|
| 1903 | ( lrd(lct(1)+1) - lrd(lct(1)) ) & |
---|
| 1904 | ) ) + ( ard(lct(1)+1) * & |
---|
| 1905 | ( ( cur_r - lrd(lct(1)) ) / & |
---|
| 1906 | ( lrd(lct(1)+1) - lrd(lct(1)) ) ) ) |
---|
[1819] | 1907 | ENDIF |
---|
| 1908 | |
---|
| 1909 | ! |
---|
[1864] | 1910 | !-- In Fortran radian instead of degree is used as unit for all |
---|
| 1911 | !-- angles. Therefore, a transformation from angles given in |
---|
| 1912 | !-- degree to angles given in radian is necessary here: |
---|
| 1913 | alpha_attack(rseg) = alpha_attack(rseg) * & |
---|
[1819] | 1914 | ( (2.0_wp*pi) / 360.0_wp ) |
---|
| 1915 | ! |
---|
[1864] | 1916 | !-- Substraction of the local pitch angle to obtain the local |
---|
| 1917 | !-- angle of attack: |
---|
[1819] | 1918 | alpha_attack(rseg) = phi_rel(rseg) - alpha_attack(rseg) |
---|
| 1919 | ! |
---|
[1864] | 1920 | !-- Preliminary transformation back from angles given in radian |
---|
| 1921 | !-- to angles given in degree: |
---|
| 1922 | alpha_attack(rseg) = alpha_attack(rseg) * & |
---|
[1819] | 1923 | ( 360.0_wp / (2.0_wp*pi) ) |
---|
| 1924 | ! |
---|
[1864] | 1925 | !-- Correct with collective pitch angle: |
---|
[2015] | 1926 | alpha_attack = alpha_attack - pitch_add(inot) |
---|
[1819] | 1927 | |
---|
| 1928 | ! |
---|
[1864] | 1929 | !-- Determination of the magnitude of the flow velocity relative |
---|
| 1930 | !-- to the rotor: |
---|
[1912] | 1931 | vrel(rseg) = SQRT( u_rot**2 + & |
---|
| 1932 | ( omega_rot(inot) * cur_r - & |
---|
[1819] | 1933 | vtheta(rseg) )**2 ) |
---|
| 1934 | |
---|
| 1935 | ! |
---|
[1864] | 1936 | !-- !-----------------------------------------------------------! |
---|
| 1937 | !-- !-- Interpolation of chord as well as lift and drag --! |
---|
| 1938 | !-- !-- coefficients from tabulated values --! |
---|
| 1939 | !-- !-----------------------------------------------------------! |
---|
[1819] | 1940 | |
---|
| 1941 | ! |
---|
[1864] | 1942 | !-- Interpolation of the chord_length from tabulated values to |
---|
| 1943 | !-- the current radial position: |
---|
[1819] | 1944 | |
---|
| 1945 | IF (cur_r == 0.0_wp) THEN |
---|
| 1946 | chord(rseg) = 0.0_wp |
---|
| 1947 | ELSE IF (cur_r >= lrd(size(crd))) THEN |
---|
[1864] | 1948 | chord(rseg) = (crd(size(crd)) + ard(size(crd)-1)) / 2.0_wp |
---|
[1819] | 1949 | ELSE |
---|
[1864] | 1950 | chord(rseg) = ( crd(lct(1)) * & |
---|
| 1951 | ( ( lrd(lct(1)+1) - cur_r ) / & |
---|
| 1952 | ( lrd(lct(1)+1) - lrd(lct(1)) ) ) ) + & |
---|
| 1953 | ( crd(lct(1)+1) * & |
---|
| 1954 | ( ( cur_r-lrd(lct(1)) ) / & |
---|
| 1955 | ( lrd(lct(1)+1) - lrd(lct(1)) ) ) ) |
---|
[1819] | 1956 | ENDIF |
---|
| 1957 | |
---|
| 1958 | ! |
---|
| 1959 | !-- Determine index of current angle of attack, needed for |
---|
[1864] | 1960 | !-- finding the appropriate interpolated values of the lift and |
---|
| 1961 | !-- drag coefficients (-180.0 degrees = 0, +180.0 degrees = 36000, |
---|
| 1962 | !-- so one index every 0.01 degrees): |
---|
| 1963 | iialpha = CEILING( ( alpha_attack(rseg) + 180.0_wp ) & |
---|
[1839] | 1964 | * ( 1.0_wp / accu_cl_cd_tab ) ) |
---|
[1819] | 1965 | ! |
---|
| 1966 | !-- Determine index of current radial position, needed for |
---|
[1864] | 1967 | !-- finding the appropriate interpolated values of the lift and |
---|
| 1968 | !-- drag coefficients (one index every 0.1 m): |
---|
[1819] | 1969 | iir = CEILING( cur_r * 10.0_wp ) |
---|
| 1970 | ! |
---|
[1864] | 1971 | !-- Read in interpolated values of the lift and drag coefficients |
---|
| 1972 | !-- for the current radial position and angle of attack: |
---|
[1839] | 1973 | turb_cl(rseg) = turb_cl_tab(iialpha,iir) |
---|
| 1974 | turb_cd(rseg) = turb_cd_tab(iialpha,iir) |
---|
[1819] | 1975 | |
---|
| 1976 | ! |
---|
[1864] | 1977 | !-- Final transformation back from angles given in degree to |
---|
| 1978 | !-- angles given in radian: |
---|
| 1979 | alpha_attack(rseg) = alpha_attack(rseg) * & |
---|
[1819] | 1980 | ( (2.0_wp*pi) / 360.0_wp ) |
---|
| 1981 | |
---|
[2152] | 1982 | IF ( tl_cor ) THEN |
---|
| 1983 | |
---|
| 1984 | !-- Tip loss correction following Schito |
---|
| 1985 | !-- Schito applies the tip loss correction only to the lift force |
---|
| 1986 | !-- Therefore, the tip loss correction is only applied to the lift |
---|
| 1987 | !-- coefficient and not to the drag coefficient in our case |
---|
| 1988 | !-- |
---|
| 1989 | tl_factor = ( 2.0 / pi ) * & |
---|
| 1990 | ACOS( EXP( -1.0 * ( 3.0 * ( rr(inot) - cur_r ) / & |
---|
| 1991 | ( 2.0 * cur_r * abs( sin( phi_rel(rseg) ) ) ) ) ) ) |
---|
| 1992 | |
---|
| 1993 | turb_cl(rseg) = tl_factor * turb_cl(rseg) |
---|
| 1994 | |
---|
| 1995 | END IF |
---|
[1819] | 1996 | ! |
---|
| 1997 | !-- !-----------------------------------------------------! |
---|
| 1998 | !-- !-- Calculation of the forces --! |
---|
| 1999 | !-- !-----------------------------------------------------! |
---|
| 2000 | |
---|
| 2001 | ! |
---|
[1864] | 2002 | !-- Calculate the pre_factor for the thrust and torque forces: |
---|
[1819] | 2003 | |
---|
| 2004 | pre_factor = 0.5_wp * (vrel(rseg)**2) * 3.0_wp * & |
---|
| 2005 | chord(rseg) * delta_r(inot) / nsegs(ring,inot) |
---|
| 2006 | |
---|
| 2007 | ! |
---|
[1864] | 2008 | !-- Calculate the thrust force (x-component of the total force) |
---|
| 2009 | !-- for each ring segment: |
---|
| 2010 | thrust_seg(rseg) = pre_factor * & |
---|
| 2011 | ( turb_cl(rseg) * COS(phi_rel(rseg)) + & |
---|
| 2012 | turb_cd(rseg) * SIN(phi_rel(rseg)) ) |
---|
[1819] | 2013 | |
---|
| 2014 | ! |
---|
[1864] | 2015 | !-- Determination of the second of the additional forces acting |
---|
| 2016 | !-- on the flow in the azimuthal direction: force vector as basis |
---|
| 2017 | !-- for torque (torque itself would be the vector product of the |
---|
| 2018 | !-- radius vector and the force vector): |
---|
| 2019 | torque_seg = pre_factor * & |
---|
| 2020 | ( turb_cl(rseg) * SIN(phi_rel(rseg)) - & |
---|
[1839] | 2021 | turb_cd(rseg) * COS(phi_rel(rseg)) ) |
---|
[1819] | 2022 | ! |
---|
| 2023 | !-- Decomposition of the force vector into two parts: |
---|
[1864] | 2024 | !-- One acting along the y-direction and one acting along the |
---|
| 2025 | !-- z-direction of the rotor coordinate system: |
---|
[1819] | 2026 | |
---|
| 2027 | torque_seg_y(rseg) = -torque_seg * sin_rot |
---|
| 2028 | torque_seg_z(rseg) = torque_seg * cos_rot |
---|
| 2029 | |
---|
[1912] | 2030 | ! |
---|
| 2031 | !-- Add the segment thrust to the thrust of the whole rotor |
---|
| 2032 | thrust_rotor(inot) = thrust_rotor(inot) + & |
---|
| 2033 | thrust_seg(rseg) |
---|
| 2034 | |
---|
| 2035 | |
---|
[1819] | 2036 | torque_total(inot) = torque_total(inot) + (torque_seg * cur_r) |
---|
| 2037 | |
---|
[1864] | 2038 | ENDDO !-- end of loop over ring segments |
---|
[1819] | 2039 | |
---|
| 2040 | ! |
---|
[1864] | 2041 | !-- Restore the forces into arrays containing all the segments of |
---|
| 2042 | !-- each ring: |
---|
[1819] | 2043 | thrust_ring(ring,:) = thrust_seg(:) |
---|
| 2044 | torque_ring_y(ring,:) = torque_seg_y(:) |
---|
| 2045 | torque_ring_z(ring,:) = torque_seg_z(:) |
---|
| 2046 | |
---|
| 2047 | |
---|
[1864] | 2048 | ENDDO !-- end of loop over rings |
---|
[1819] | 2049 | |
---|
| 2050 | |
---|
[1864] | 2051 | CALL cpu_log( log_point_s(62), 'wtm_controller', 'start' ) |
---|
[1819] | 2052 | |
---|
[1912] | 2053 | |
---|
| 2054 | IF ( speed_control ) THEN |
---|
| 2055 | ! |
---|
| 2056 | !-- Calculation of the current generator speed for rotor speed control |
---|
| 2057 | |
---|
| 2058 | ! |
---|
| 2059 | !-- The acceleration of the rotor speed is calculated from |
---|
| 2060 | !-- the force balance of the accelerating torque |
---|
| 2061 | !-- and the torque of the rotating rotor and generator |
---|
| 2062 | om_rate = ( torque_total(inot) * air_dens * gear_eff - & |
---|
| 2063 | gear_ratio * torque_gen_old(inot) ) / & |
---|
| 2064 | ( inertia_rot + & |
---|
| 2065 | gear_ratio * gear_ratio * inertia_gen ) * dt_3d |
---|
[1819] | 2066 | |
---|
[1912] | 2067 | ! |
---|
| 2068 | !-- The generator speed is given by the product of gear gear_ratio |
---|
| 2069 | !-- and rotor speed |
---|
| 2070 | omega_gen(inot) = gear_ratio * ( omega_rot(inot) + om_rate ) |
---|
| 2071 | |
---|
| 2072 | ENDIF |
---|
| 2073 | |
---|
[1864] | 2074 | IF ( pitch_control ) THEN |
---|
[1819] | 2075 | |
---|
[1912] | 2076 | ! |
---|
| 2077 | !-- If the current generator speed is above rated, the pitch is not |
---|
| 2078 | !-- saturated and the change from the last time step is within the |
---|
| 2079 | !-- maximum pitch rate, then the pitch loop is repeated with a pitch |
---|
| 2080 | !-- gain |
---|
| 2081 | IF ( ( omega_gen(inot) > rated_genspeed ) .AND. & |
---|
| 2082 | ( pitch_add(inot) < 25.0_wp ) .AND. & |
---|
| 2083 | ( pitch_add(inot) < pitch_add_old(inot) + & |
---|
| 2084 | pitch_rate * dt_3d ) ) THEN |
---|
[1864] | 2085 | pitch_sw = .TRUE. |
---|
[1912] | 2086 | ! |
---|
| 2087 | !-- Go back to beginning of pit_loop |
---|
| 2088 | CYCLE pit_loop |
---|
[1819] | 2089 | ENDIF |
---|
[1912] | 2090 | |
---|
| 2091 | ! |
---|
| 2092 | !-- The current pitch is saved for the next time step |
---|
| 2093 | pitch_add_old(inot) = pitch_add(inot) |
---|
[1864] | 2094 | pitch_sw = .FALSE. |
---|
[1819] | 2095 | ENDIF |
---|
[1912] | 2096 | EXIT pit_loop |
---|
| 2097 | ENDDO pit_loop ! Recursive pitch control loop |
---|
[1819] | 2098 | |
---|
[1864] | 2099 | |
---|
[1819] | 2100 | ! |
---|
[1864] | 2101 | !-- Call the rotor speed controller |
---|
| 2102 | |
---|
[1819] | 2103 | IF ( speed_control ) THEN |
---|
| 2104 | ! |
---|
[1864] | 2105 | !-- Find processor at i_hub, j_hub |
---|
[1912] | 2106 | IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) & |
---|
| 2107 | THEN |
---|
| 2108 | IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) )& |
---|
| 2109 | THEN |
---|
[1864] | 2110 | CALL wtm_speed_control( inot ) |
---|
| 2111 | ENDIF |
---|
[1912] | 2112 | ENDIF |
---|
[1864] | 2113 | |
---|
[1819] | 2114 | ENDIF |
---|
| 2115 | |
---|
| 2116 | |
---|
[1864] | 2117 | CALL cpu_log( log_point_s(62), 'wtm_controller', 'stop' ) |
---|
[1819] | 2118 | |
---|
[1864] | 2119 | CALL cpu_log( log_point_s(63), 'wtm_smearing', 'start' ) |
---|
[1819] | 2120 | |
---|
| 2121 | |
---|
[1864] | 2122 | !-- !-----------------------------------------------------------------! |
---|
| 2123 | !-- !-- Regularization kernel --! |
---|
| 2124 | !-- !-- Smearing of the forces and interpolation to cartesian grid --! |
---|
| 2125 | !-- !-----------------------------------------------------------------! |
---|
[1819] | 2126 | ! |
---|
[1864] | 2127 | !-- The aerodynamic blade forces need to be distributed smoothly on |
---|
| 2128 | !-- several mesh points in order to avoid singular behaviour |
---|
[1819] | 2129 | ! |
---|
| 2130 | !-- Summation over sum of weighted forces. The weighting factor |
---|
[1864] | 2131 | !-- (calculated in user_init) includes information on the distance |
---|
| 2132 | !-- between the center of the grid cell and the rotor segment under |
---|
| 2133 | !-- consideration |
---|
[1819] | 2134 | ! |
---|
[1864] | 2135 | !-- To save computing time, apply smearing only for the relevant part |
---|
| 2136 | !-- of the model domain: |
---|
[1819] | 2137 | ! |
---|
| 2138 | !-- |
---|
| 2139 | !-- Calculation of the boundaries: |
---|
[1864] | 2140 | i_smear(inot) = CEILING( ( rr(inot) * ABS( roty(inot,1) ) + & |
---|
| 2141 | eps_min ) / dx ) |
---|
| 2142 | j_smear(inot) = CEILING( ( rr(inot) * ABS( roty(inot,2) ) + & |
---|
| 2143 | eps_min ) / dy ) |
---|
[1819] | 2144 | |
---|
[1864] | 2145 | DO i = MAX( nxl, i_hub(inot) - i_smear(inot) ), & |
---|
[1819] | 2146 | MIN( nxr, i_hub(inot) + i_smear(inot) ) |
---|
[1864] | 2147 | DO j = MAX( nys, j_hub(inot) - j_smear(inot) ), & |
---|
[1819] | 2148 | MIN( nyn, j_hub(inot) + j_smear(inot) ) |
---|
[2232] | 2149 | ! DO k = MAX( nzb_u_inner(j,i)+1, k_hub(inot) - k_smear(inot) ), & |
---|
| 2150 | ! k_hub(inot) + k_smear(inot) |
---|
| 2151 | DO k = nzb+1, k_hub(inot) + k_smear(inot) |
---|
[1819] | 2152 | DO ring = 1, nrings(inot) |
---|
| 2153 | DO rseg = 1, nsegs(ring,inot) |
---|
| 2154 | ! |
---|
[2232] | 2155 | !-- Predetermine flag to mask topography |
---|
| 2156 | flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1 ) ) |
---|
| 2157 | |
---|
| 2158 | ! |
---|
[1819] | 2159 | !-- Determine the square of the distance between the |
---|
| 2160 | !-- current grid point and each rotor area segment: |
---|
[1864] | 2161 | dist_u_3d = ( i * dx - rbx(ring,rseg) )**2 + & |
---|
| 2162 | ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & |
---|
| 2163 | ( k * dz - 0.5_wp * dz - rbz(ring,rseg) )**2 |
---|
| 2164 | dist_v_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & |
---|
| 2165 | ( j * dy - rby(ring,rseg) )**2 + & |
---|
| 2166 | ( k * dz - 0.5_wp * dz - rbz(ring,rseg) )**2 |
---|
| 2167 | dist_w_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & |
---|
| 2168 | ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & |
---|
| 2169 | ( k * dz - rbz(ring,rseg) )**2 |
---|
[1819] | 2170 | |
---|
| 2171 | ! |
---|
| 2172 | !-- 3D-smearing of the forces with a polynomial function |
---|
| 2173 | !-- (much faster than the old Gaussian function), using |
---|
| 2174 | !-- some parameters that have been calculated in user_init. |
---|
| 2175 | !-- The function is only similar to Gaussian function for |
---|
| 2176 | !-- squared distances <= eps_min2: |
---|
| 2177 | IF ( dist_u_3d <= eps_min2 ) THEN |
---|
[2232] | 2178 | thrust(k,j,i) = thrust(k,j,i) + & |
---|
| 2179 | thrust_ring(ring,rseg) * & |
---|
| 2180 | ( ( pol_a * dist_u_3d - pol_b ) * & |
---|
| 2181 | dist_u_3d + 1.0_wp ) * eps_factor *& |
---|
| 2182 | flag |
---|
[1819] | 2183 | ENDIF |
---|
| 2184 | IF ( dist_v_3d <= eps_min2 ) THEN |
---|
[2232] | 2185 | torque_y(k,j,i) = torque_y(k,j,i) + & |
---|
| 2186 | torque_ring_y(ring,rseg) * & |
---|
| 2187 | ( ( pol_a * dist_v_3d - pol_b ) * & |
---|
| 2188 | dist_v_3d + 1.0_wp ) * eps_factor * & |
---|
| 2189 | flag |
---|
[1819] | 2190 | ENDIF |
---|
| 2191 | IF ( dist_w_3d <= eps_min2 ) THEN |
---|
[2232] | 2192 | torque_z(k,j,i) = torque_z(k,j,i) + & |
---|
| 2193 | torque_ring_z(ring,rseg) * & |
---|
| 2194 | ( ( pol_a * dist_w_3d - pol_b ) * & |
---|
| 2195 | dist_w_3d + 1.0_wp ) * eps_factor * & |
---|
| 2196 | flag |
---|
[1819] | 2197 | ENDIF |
---|
| 2198 | |
---|
| 2199 | ENDDO ! End of loop over rseg |
---|
| 2200 | ENDDO ! End of loop over ring |
---|
| 2201 | |
---|
| 2202 | ! |
---|
| 2203 | !-- Rotation of force components: |
---|
[2232] | 2204 | rot_tend_x(k,j,i) = rot_tend_x(k,j,i) + ( & |
---|
[1864] | 2205 | thrust(k,j,i)*rotx(inot,1) + & |
---|
| 2206 | torque_y(k,j,i)*roty(inot,1) + & |
---|
[2232] | 2207 | torque_z(k,j,i)*rotz(inot,1) & |
---|
| 2208 | ) * flag |
---|
[1819] | 2209 | |
---|
[2232] | 2210 | rot_tend_y(k,j,i) = rot_tend_y(k,j,i) + ( & |
---|
[1864] | 2211 | thrust(k,j,i)*rotx(inot,2) + & |
---|
| 2212 | torque_y(k,j,i)*roty(inot,2) + & |
---|
[2232] | 2213 | torque_z(k,j,i)*rotz(inot,2) & |
---|
| 2214 | ) * flag |
---|
[1819] | 2215 | |
---|
[2232] | 2216 | rot_tend_z(k,j,i) = rot_tend_z(k,j,i) + ( & |
---|
[1864] | 2217 | thrust(k,j,i)*rotx(inot,3) + & |
---|
| 2218 | torque_y(k,j,i)*roty(inot,3) + & |
---|
[2232] | 2219 | torque_z(k,j,i)*rotz(inot,3) & |
---|
| 2220 | ) * flag |
---|
[1819] | 2221 | |
---|
| 2222 | ENDDO ! End of loop over k |
---|
| 2223 | ENDDO ! End of loop over j |
---|
| 2224 | ENDDO ! End of loop over i |
---|
| 2225 | |
---|
[1864] | 2226 | CALL cpu_log( log_point_s(63), 'wtm_smearing', 'stop' ) |
---|
[1819] | 2227 | |
---|
[1912] | 2228 | ENDDO !-- end of loop over turbines |
---|
[1819] | 2229 | |
---|
| 2230 | |
---|
| 2231 | IF ( yaw_control ) THEN |
---|
| 2232 | ! |
---|
| 2233 | !-- Allocate arrays for yaw control at first call |
---|
| 2234 | !-- Can't be allocated before dt_3d is set |
---|
| 2235 | IF ( start_up ) THEN |
---|
[1864] | 2236 | WDLON = NINT( 30.0_wp / dt_3d ) ! 30s running mean array |
---|
[1819] | 2237 | ALLOCATE( wd30(1:nturbines,1:WDLON) ) |
---|
[1912] | 2238 | wd30 = 999.0_wp ! Set to dummy value |
---|
[1819] | 2239 | ALLOCATE( wd30_l(1:WDLON) ) |
---|
| 2240 | |
---|
[1864] | 2241 | WDSHO = NINT( 2.0_wp / dt_3d ) ! 2s running mean array |
---|
[1819] | 2242 | ALLOCATE( wd2(1:nturbines,1:WDSHO) ) |
---|
[1912] | 2243 | wd2 = 999.0_wp ! Set to dummy value |
---|
[1819] | 2244 | ALLOCATE( wd2_l(1:WDSHO) ) |
---|
| 2245 | start_up = .FALSE. |
---|
| 2246 | ENDIF |
---|
| 2247 | |
---|
| 2248 | ! |
---|
| 2249 | !-- Calculate the inflow wind speed |
---|
| 2250 | !-- |
---|
| 2251 | !-- Loop over number of turbines: |
---|
| 2252 | DO inot = 1, nturbines |
---|
| 2253 | ! |
---|
| 2254 | !-- Find processor at i_hub, j_hub |
---|
[1912] | 2255 | IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) & |
---|
| 2256 | THEN |
---|
| 2257 | IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) )& |
---|
| 2258 | THEN |
---|
[1864] | 2259 | |
---|
[1819] | 2260 | u_inflow_l(inot) = u(k_hub(inot),j_hub(inot),i_hub(inot)) |
---|
[1864] | 2261 | |
---|
[1912] | 2262 | wdir_l(inot) = -1.0_wp * ATAN2( & |
---|
| 2263 | 0.5_wp * ( v(k_hub(inot),j_hub(inot),i_hub(inot)+1) + & |
---|
| 2264 | v(k_hub(inot),j_hub(inot),i_hub(inot)) ) , & |
---|
| 2265 | 0.5_wp * ( u(k_hub(inot),j_hub(inot)+1,i_hub(inot)) + & |
---|
| 2266 | u(k_hub(inot),j_hub(inot),i_hub(inot)) ) ) |
---|
[1864] | 2267 | |
---|
| 2268 | CALL wtm_yawcontrol( inot ) |
---|
| 2269 | |
---|
[1819] | 2270 | phi_yaw_l(inot) = phi_yaw(inot) |
---|
[2015] | 2271 | |
---|
[1819] | 2272 | ENDIF |
---|
| 2273 | ENDIF |
---|
| 2274 | |
---|
[1864] | 2275 | ENDDO !-- end of loop over turbines |
---|
| 2276 | |
---|
| 2277 | ! |
---|
[1929] | 2278 | !-- Transfer of information to the other cpus |
---|
| 2279 | #if defined( __parallel ) |
---|
[1864] | 2280 | CALL MPI_ALLREDUCE( u_inflow_l, u_inflow, nturbines, MPI_REAL, & |
---|
| 2281 | MPI_SUM, comm2d, ierr ) |
---|
| 2282 | CALL MPI_ALLREDUCE( wdir_l, wdir, nturbines, MPI_REAL, MPI_SUM, & |
---|
| 2283 | comm2d, ierr ) |
---|
| 2284 | CALL MPI_ALLREDUCE( phi_yaw_l, phi_yaw, nturbines, MPI_REAL, & |
---|
| 2285 | MPI_SUM, comm2d, ierr ) |
---|
[1929] | 2286 | #else |
---|
| 2287 | u_inflow = u_inflow_l |
---|
| 2288 | wdir = wdir_l |
---|
| 2289 | phi_yaw = phi_yaw_l |
---|
[2015] | 2290 | |
---|
| 2291 | |
---|
[1929] | 2292 | #endif |
---|
[1819] | 2293 | DO inot = 1, nturbines |
---|
[1864] | 2294 | ! |
---|
| 2295 | !-- Update rotor orientation |
---|
| 2296 | CALL wtm_rotate_rotor( inot ) |
---|
[1819] | 2297 | |
---|
| 2298 | ENDDO ! End of loop over turbines |
---|
| 2299 | |
---|
[2152] | 2300 | END IF ! end of yaw control |
---|
[1864] | 2301 | |
---|
| 2302 | IF ( speed_control ) THEN |
---|
| 2303 | ! |
---|
| 2304 | !-- Transfer of information to the other cpus |
---|
[1912] | 2305 | ! CALL MPI_ALLREDUCE( omega_gen, omega_gen_old, nturbines, & |
---|
| 2306 | ! MPI_REAL,MPI_SUM, comm2d, ierr ) |
---|
[1929] | 2307 | #if defined( __parallel ) |
---|
[1864] | 2308 | CALL MPI_ALLREDUCE( torque_gen, torque_gen_old, nturbines, & |
---|
| 2309 | MPI_REAL, MPI_SUM, comm2d, ierr ) |
---|
| 2310 | CALL MPI_ALLREDUCE( omega_rot_l, omega_rot, nturbines, & |
---|
| 2311 | MPI_REAL, MPI_SUM, comm2d, ierr ) |
---|
[1912] | 2312 | CALL MPI_ALLREDUCE( omega_gen_f, omega_gen_f_old, nturbines, & |
---|
| 2313 | MPI_REAL, MPI_SUM, comm2d, ierr ) |
---|
[1929] | 2314 | #else |
---|
| 2315 | torque_gen_old = torque_gen |
---|
| 2316 | omega_rot = omega_rot_l |
---|
| 2317 | omega_gen_f_old = omega_gen_f |
---|
| 2318 | #endif |
---|
[1864] | 2319 | |
---|
| 2320 | ENDIF |
---|
[1819] | 2321 | |
---|
| 2322 | DO inot = 1, nturbines |
---|
| 2323 | |
---|
| 2324 | IF ( myid == 0 ) THEN |
---|
| 2325 | IF ( openfile_turb_mod(400+inot)%opened ) THEN |
---|
| 2326 | WRITE ( 400+inot, 106 ) simulated_time, omega_rot(inot), & |
---|
[1912] | 2327 | omega_gen(inot), torque_gen_old(inot), & |
---|
[1864] | 2328 | torque_total(inot), pitch_add(inot), & |
---|
[1912] | 2329 | torque_gen_old(inot)*omega_gen(inot)*gen_eff, & |
---|
[1864] | 2330 | torque_total(inot)*omega_rot(inot)*air_dens, & |
---|
[1912] | 2331 | thrust_rotor(inot), & |
---|
[1864] | 2332 | wdir(inot)*180.0_wp/pi, & |
---|
| 2333 | (phi_yaw(inot))*180.0_wp/pi |
---|
[1912] | 2334 | |
---|
[1819] | 2335 | ELSE |
---|
| 2336 | |
---|
| 2337 | WRITE ( turbine_id,'(I2.2)') inot |
---|
[1864] | 2338 | OPEN ( 400+inot, FILE=( 'TURBINE_PARAMETERS'//turbine_id ), & |
---|
| 2339 | FORM='FORMATTED' ) |
---|
[1819] | 2340 | WRITE ( 400+inot, 105 ) inot |
---|
| 2341 | WRITE ( 400+inot, 106 ) simulated_time, omega_rot(inot), & |
---|
[1912] | 2342 | omega_gen(inot), torque_gen_old(inot), & |
---|
[1864] | 2343 | torque_total(inot), pitch_add(inot), & |
---|
[1912] | 2344 | torque_gen_old(inot)*omega_gen(inot)*gen_eff, & |
---|
[1864] | 2345 | torque_total(inot)*omega_rot(inot)*air_dens, & |
---|
[1912] | 2346 | thrust_rotor(inot), & |
---|
[1864] | 2347 | wdir(inot)*180.0_wp/pi, & |
---|
| 2348 | (phi_yaw(inot))*180.0_wp/pi |
---|
[1819] | 2349 | ENDIF |
---|
| 2350 | ENDIF |
---|
| 2351 | |
---|
[1912] | 2352 | !-- Set open flag |
---|
[1819] | 2353 | openfile_turb_mod(400+inot)%opened = .TRUE. |
---|
[1864] | 2354 | ENDDO !-- end of loop over turbines |
---|
[1819] | 2355 | |
---|
| 2356 | ENDIF |
---|
| 2357 | |
---|
[1864] | 2358 | CALL cpu_log( log_point_s(61), 'wtm_forces', 'stop' ) |
---|
[1819] | 2359 | |
---|
[1912] | 2360 | ! |
---|
| 2361 | !-- Formats |
---|
[1819] | 2362 | 105 FORMAT ('Turbine control data for turbine ',I2,1X,':'/ & |
---|
| 2363 | &'----------------------------------------'/ & |
---|
| 2364 | &' Time RSpeed GSpeed ', & |
---|
[1912] | 2365 | 'GenTorque AeroTorque Pitch Power(Gen) Power(Rot) ', & |
---|
| 2366 | 'RotThrust WDirection YawOrient') |
---|
[1819] | 2367 | |
---|
[1864] | 2368 | 106 FORMAT (F9.3,2X,F7.3,2X,F7.2,2X,F9.1,3X,F9.1,1X,F6.2,2X,F10.1,2X, & |
---|
| 2369 | F10.1,1X,F9.1,2X,F7.2,1X,F7.2) |
---|
[1819] | 2370 | |
---|
| 2371 | |
---|
| 2372 | END SUBROUTINE wtm_forces |
---|
| 2373 | |
---|
| 2374 | |
---|
[1839] | 2375 | !------------------------------------------------------------------------------! |
---|
| 2376 | ! Description: |
---|
| 2377 | ! ------------ |
---|
| 2378 | !> Yaw controller for the wind turbine model |
---|
| 2379 | !------------------------------------------------------------------------------! |
---|
| 2380 | SUBROUTINE wtm_yawcontrol( inot ) |
---|
| 2381 | |
---|
[1819] | 2382 | USE kinds |
---|
[1839] | 2383 | |
---|
| 2384 | IMPLICIT NONE |
---|
[1819] | 2385 | |
---|
| 2386 | INTEGER(iwp) :: inot |
---|
| 2387 | INTEGER(iwp) :: i_wd_30 |
---|
| 2388 | REAL(wp) :: missal |
---|
| 2389 | |
---|
| 2390 | i_wd_30 = 0_iwp |
---|
| 2391 | |
---|
| 2392 | ! |
---|
[1864] | 2393 | !-- The yaw controller computes a 30s running mean of the wind direction. |
---|
| 2394 | !-- If the difference between turbine alignment and wind direction exceeds |
---|
| 2395 | !-- 5°, the turbine is yawed. The mechanism stops as soon as the 2s-running |
---|
| 2396 | !-- mean of the missalignment is smaller than 0.5°. |
---|
| 2397 | !-- Attention: If the timestep during the simulation changes significantly |
---|
| 2398 | !-- the lengths of the running means change and it does not correspond to |
---|
| 2399 | !-- 30s/2s anymore. |
---|
| 2400 | !-- ! Needs to be modified for these situations ! |
---|
| 2401 | !-- For wind from the east, the averaging of the wind direction could cause |
---|
| 2402 | !-- problems and the yaw controller is probably flawed. -> Routine for |
---|
| 2403 | !-- averaging needs to be improved! |
---|
| 2404 | ! |
---|
| 2405 | !-- Check if turbine is not yawing |
---|
[1819] | 2406 | IF ( .NOT. doyaw(inot) ) THEN |
---|
[1843] | 2407 | ! |
---|
[1864] | 2408 | !-- Write current wind direction into array |
---|
[1843] | 2409 | wd30_l = wd30(inot,:) |
---|
| 2410 | wd30_l = CSHIFT( wd30_l, SHIFT=-1 ) |
---|
[1819] | 2411 | wd30_l(1) = wdir(inot) |
---|
[1843] | 2412 | ! |
---|
[1864] | 2413 | !-- Check if array is full ( no more dummies ) |
---|
[1819] | 2414 | IF ( .NOT. ANY( wd30_l == 999.) ) THEN |
---|
| 2415 | |
---|
| 2416 | missal = SUM( wd30_l ) / SIZE( wd30_l ) - phi_yaw(inot) |
---|
| 2417 | ! |
---|
[1864] | 2418 | !-- Check if turbine is missaligned by more than max_miss |
---|
[1843] | 2419 | IF ( ABS( missal ) > max_miss ) THEN |
---|
| 2420 | ! |
---|
[1864] | 2421 | !-- Check in which direction to yaw |
---|
[1843] | 2422 | yawdir(inot) = SIGN( 1.0_wp, missal ) |
---|
[1819] | 2423 | ! |
---|
[1864] | 2424 | !-- Start yawing of turbine |
---|
[1843] | 2425 | phi_yaw(inot) = phi_yaw(inot) + yawdir(inot) * yaw_speed * dt_3d |
---|
[1819] | 2426 | doyaw(inot) = .TRUE. |
---|
[1864] | 2427 | wd30_l = 999. ! fill with dummies again |
---|
[1819] | 2428 | ENDIF |
---|
| 2429 | ENDIF |
---|
| 2430 | |
---|
| 2431 | wd30(inot,:) = wd30_l |
---|
| 2432 | |
---|
| 2433 | ! |
---|
[1864] | 2434 | !-- If turbine is already yawing: |
---|
| 2435 | !-- Initialize 2 s running mean and yaw until the missalignment is smaller |
---|
| 2436 | !-- than min_miss |
---|
[1819] | 2437 | |
---|
| 2438 | ELSE |
---|
| 2439 | ! |
---|
| 2440 | !-- Initialize 2 s running mean |
---|
| 2441 | wd2_l = wd2(inot,:) |
---|
[1864] | 2442 | wd2_l = CSHIFT( wd2_l, SHIFT = -1 ) |
---|
[1819] | 2443 | wd2_l(1) = wdir(inot) |
---|
[1843] | 2444 | ! |
---|
[1864] | 2445 | !-- Check if array is full ( no more dummies ) |
---|
| 2446 | IF ( .NOT. ANY( wd2_l == 999.0_wp ) ) THEN |
---|
| 2447 | ! |
---|
| 2448 | !-- Calculate missalignment of turbine |
---|
[1819] | 2449 | missal = SUM( wd2_l - phi_yaw(inot) ) / SIZE( wd2_l ) |
---|
[1864] | 2450 | ! |
---|
| 2451 | !-- Check if missalignment is still larger than 0.5 degree and if the |
---|
| 2452 | !-- yaw direction is still right |
---|
| 2453 | IF ( ( ABS( missal ) > min_miss ) .AND. & |
---|
| 2454 | ( yawdir(inot) == SIGN( 1.0_wp, missal ) ) ) THEN |
---|
| 2455 | ! |
---|
| 2456 | !-- Continue yawing |
---|
| 2457 | phi_yaw(inot) = phi_yaw(inot) + yawdir(inot) * yaw_speed * dt_3d |
---|
[1819] | 2458 | ELSE |
---|
[1864] | 2459 | ! |
---|
| 2460 | !-- Stop yawing |
---|
[1819] | 2461 | doyaw(inot) = .FALSE. |
---|
[1864] | 2462 | wd2_l = 999.0_wp ! fill with dummies again |
---|
[1819] | 2463 | ENDIF |
---|
| 2464 | ELSE |
---|
[1864] | 2465 | ! |
---|
| 2466 | !-- Continue yawing |
---|
[1843] | 2467 | phi_yaw(inot) = phi_yaw(inot) + yawdir(inot) * yaw_speed * dt_3d |
---|
[1819] | 2468 | ENDIF |
---|
| 2469 | |
---|
| 2470 | wd2(inot,:) = wd2_l |
---|
| 2471 | |
---|
| 2472 | ENDIF |
---|
| 2473 | |
---|
[1839] | 2474 | END SUBROUTINE wtm_yawcontrol |
---|
[1819] | 2475 | |
---|
[1864] | 2476 | |
---|
[1819] | 2477 | !------------------------------------------------------------------------------! |
---|
| 2478 | ! Description: |
---|
| 2479 | ! ------------ |
---|
[1864] | 2480 | !> Initialization of the speed control |
---|
| 2481 | !------------------------------------------------------------------------------! |
---|
| 2482 | SUBROUTINE wtm_init_speed_control |
---|
| 2483 | |
---|
| 2484 | |
---|
| 2485 | IMPLICIT NONE |
---|
| 2486 | |
---|
| 2487 | ! |
---|
| 2488 | !-- If speed control is set, remaining variables and control_parameters for |
---|
| 2489 | !-- the control algorithm are calculated |
---|
| 2490 | ! |
---|
| 2491 | !-- Calculate slope constant for region 15 |
---|
| 2492 | slope15 = ( slope2 * min_reg2 * min_reg2 ) / ( min_reg2 - min_reg15 ) |
---|
| 2493 | ! |
---|
| 2494 | !-- Calculate upper limit of slipage region |
---|
| 2495 | vs_sysp = rated_genspeed / 1.1_wp |
---|
| 2496 | ! |
---|
| 2497 | !-- Calculate slope of slipage region |
---|
| 2498 | slope25 = ( rated_power / rated_genspeed ) / & |
---|
| 2499 | ( rated_genspeed - vs_sysp ) |
---|
| 2500 | ! |
---|
| 2501 | !-- Calculate lower limit of slipage region |
---|
| 2502 | min_reg25 = ( slope25 - SQRT( slope25 * ( slope25 - 4.0_wp * & |
---|
| 2503 | slope2 * vs_sysp ) ) ) / & |
---|
[2410] | 2504 | ( 2.0_wp * slope2 ) |
---|
[1864] | 2505 | ! |
---|
| 2506 | !-- Frequency for the simple low pass filter |
---|
| 2507 | Fcorner = 0.25_wp |
---|
| 2508 | ! |
---|
| 2509 | !-- At the first timestep the torque is set to its maximum to prevent |
---|
| 2510 | !-- an overspeeding of the rotor |
---|
[2563] | 2511 | IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN |
---|
| 2512 | torque_gen_old(:) = max_torque_gen |
---|
| 2513 | ENDIF |
---|
[1864] | 2514 | |
---|
| 2515 | END SUBROUTINE wtm_init_speed_control |
---|
| 2516 | |
---|
| 2517 | |
---|
| 2518 | !------------------------------------------------------------------------------! |
---|
| 2519 | ! Description: |
---|
| 2520 | ! ------------ |
---|
| 2521 | !> Simple controller for the regulation of the rotor speed |
---|
| 2522 | !------------------------------------------------------------------------------! |
---|
| 2523 | SUBROUTINE wtm_speed_control( inot ) |
---|
| 2524 | |
---|
| 2525 | |
---|
| 2526 | IMPLICIT NONE |
---|
| 2527 | |
---|
[1912] | 2528 | INTEGER(iwp) :: inot |
---|
| 2529 | |
---|
| 2530 | |
---|
[1864] | 2531 | |
---|
| 2532 | ! |
---|
| 2533 | !-- The controller is based on the fortran script from Jonkman |
---|
| 2534 | !-- et al. 2009 "Definition of a 5 MW Reference Wind Turbine for |
---|
| 2535 | !-- offshore system developement" |
---|
| 2536 | |
---|
| 2537 | ! |
---|
| 2538 | !-- The generator speed is filtered by a low pass filter |
---|
| 2539 | !-- for the control of the generator torque |
---|
| 2540 | lp_coeff = EXP( -2.0_wp * 3.14_wp * dt_3d * Fcorner ) |
---|
[1912] | 2541 | omega_gen_f(inot) = ( 1.0_wp - lp_coeff ) * omega_gen(inot) + lp_coeff *& |
---|
[1864] | 2542 | omega_gen_f_old(inot) |
---|
| 2543 | |
---|
| 2544 | IF ( omega_gen_f(inot) <= min_reg15 ) THEN |
---|
| 2545 | ! |
---|
| 2546 | !-- Region 1: Generator torque is set to zero to accelerate the rotor: |
---|
| 2547 | torque_gen(inot) = 0 |
---|
| 2548 | |
---|
| 2549 | ELSEIF ( omega_gen_f(inot) <= min_reg2 ) THEN |
---|
| 2550 | ! |
---|
| 2551 | !-- Region 1.5: Generator torque is increasing linearly with rotor speed: |
---|
| 2552 | torque_gen(inot) = slope15 * ( omega_gen_f(inot) - min_reg15 ) |
---|
| 2553 | |
---|
| 2554 | ELSEIF ( omega_gen_f(inot) <= min_reg25 ) THEN |
---|
| 2555 | ! |
---|
| 2556 | !-- Region 2: Generator torque is increased by the square of the generator |
---|
| 2557 | !-- speed to keep the TSR optimal: |
---|
| 2558 | torque_gen(inot) = slope2 * omega_gen_f(inot) * omega_gen_f(inot) |
---|
| 2559 | |
---|
| 2560 | ELSEIF ( omega_gen_f(inot) < rated_genspeed ) THEN |
---|
| 2561 | ! |
---|
| 2562 | !-- Region 2.5: Slipage region between 2 and 3: |
---|
| 2563 | torque_gen(inot) = slope25 * ( omega_gen_f(inot) - vs_sysp ) |
---|
| 2564 | |
---|
| 2565 | ELSE |
---|
| 2566 | ! |
---|
| 2567 | !-- Region 3: Generator torque is antiproportional to the rotor speed to |
---|
| 2568 | !-- keep the power constant: |
---|
| 2569 | torque_gen(inot) = rated_power / omega_gen_f(inot) |
---|
| 2570 | |
---|
| 2571 | ENDIF |
---|
| 2572 | ! |
---|
| 2573 | !-- Calculate torque rate and confine with a max |
---|
| 2574 | trq_rate = ( torque_gen(inot) - torque_gen_old(inot) ) / dt_3d |
---|
| 2575 | trq_rate = MIN( MAX( trq_rate, -1.0_wp * max_trq_rate ), max_trq_rate ) |
---|
| 2576 | ! |
---|
| 2577 | !-- Calculate new gen torque and confine with max torque |
---|
| 2578 | torque_gen(inot) = torque_gen_old(inot) + trq_rate * dt_3d |
---|
| 2579 | torque_gen(inot) = MIN( torque_gen(inot), max_torque_gen ) |
---|
| 2580 | ! |
---|
| 2581 | !-- Overwrite values for next timestep |
---|
[1912] | 2582 | omega_rot_l(inot) = omega_gen(inot) / gear_ratio |
---|
[1864] | 2583 | |
---|
| 2584 | |
---|
| 2585 | END SUBROUTINE wtm_speed_control |
---|
| 2586 | |
---|
| 2587 | |
---|
| 2588 | !------------------------------------------------------------------------------! |
---|
| 2589 | ! Description: |
---|
| 2590 | ! ------------ |
---|
[1839] | 2591 | !> Application of the additional forces generated by the wind turbine on the |
---|
| 2592 | !> flow components (tendency terms) |
---|
| 2593 | !> Call for all grid points |
---|
[1819] | 2594 | !------------------------------------------------------------------------------! |
---|
| 2595 | SUBROUTINE wtm_tendencies( component ) |
---|
| 2596 | |
---|
[1839] | 2597 | |
---|
[1819] | 2598 | IMPLICIT NONE |
---|
| 2599 | |
---|
[1839] | 2600 | INTEGER(iwp) :: component !< prognostic variable (u,v,w) |
---|
| 2601 | INTEGER(iwp) :: i !< running index |
---|
| 2602 | INTEGER(iwp) :: j !< running index |
---|
| 2603 | INTEGER(iwp) :: k !< running index |
---|
[1819] | 2604 | |
---|
| 2605 | |
---|
| 2606 | SELECT CASE ( component ) |
---|
| 2607 | |
---|
| 2608 | CASE ( 1 ) |
---|
| 2609 | ! |
---|
| 2610 | !-- Apply the x-component of the force to the u-component of the flow: |
---|
[1864] | 2611 | IF ( simulated_time >= time_turbine_on ) THEN |
---|
[1819] | 2612 | DO i = nxlg, nxrg |
---|
| 2613 | DO j = nysg, nyng |
---|
[2553] | 2614 | DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) |
---|
[1819] | 2615 | ! |
---|
[1864] | 2616 | !-- Calculate the thrust generated by the nacelle and the tower |
---|
[1912] | 2617 | tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) * & |
---|
[2553] | 2618 | SIGN( u(k,j,i)**2 , u(k,j,i) ) |
---|
[1912] | 2619 | tend_tow_x = 0.5_wp * tow_cd_surf(k,j,i) * & |
---|
[1819] | 2620 | SIGN( u(k,j,i)**2 , u(k,j,i) ) |
---|
[2553] | 2621 | |
---|
[2232] | 2622 | tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i) & |
---|
| 2623 | - tend_nac_x - tend_tow_x ) & |
---|
| 2624 | * MERGE( 1.0_wp, 0.0_wp, & |
---|
| 2625 | BTEST( wall_flags_0(k,j,i), 1 ) ) |
---|
[1819] | 2626 | ENDDO |
---|
| 2627 | ENDDO |
---|
| 2628 | ENDDO |
---|
| 2629 | ENDIF |
---|
| 2630 | |
---|
| 2631 | CASE ( 2 ) |
---|
| 2632 | ! |
---|
| 2633 | !-- Apply the y-component of the force to the v-component of the flow: |
---|
[1864] | 2634 | IF ( simulated_time >= time_turbine_on ) THEN |
---|
[1819] | 2635 | DO i = nxlg, nxrg |
---|
| 2636 | DO j = nysg, nyng |
---|
[2553] | 2637 | DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) |
---|
[1912] | 2638 | tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * & |
---|
[1819] | 2639 | SIGN( v(k,j,i)**2 , v(k,j,i) ) |
---|
[1912] | 2640 | tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * & |
---|
[1819] | 2641 | SIGN( v(k,j,i)**2 , v(k,j,i) ) |
---|
[2232] | 2642 | tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) & |
---|
| 2643 | - tend_nac_y - tend_tow_y ) & |
---|
| 2644 | * MERGE( 1.0_wp, 0.0_wp, & |
---|
| 2645 | BTEST( wall_flags_0(k,j,i), 2 ) ) |
---|
[1819] | 2646 | ENDDO |
---|
| 2647 | ENDDO |
---|
| 2648 | ENDDO |
---|
| 2649 | ENDIF |
---|
| 2650 | |
---|
| 2651 | CASE ( 3 ) |
---|
| 2652 | ! |
---|
| 2653 | !-- Apply the z-component of the force to the w-component of the flow: |
---|
[1864] | 2654 | IF ( simulated_time >= time_turbine_on ) THEN |
---|
[1819] | 2655 | DO i = nxlg, nxrg |
---|
| 2656 | DO j = nysg, nyng |
---|
[2553] | 2657 | DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) |
---|
[2232] | 2658 | tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) & |
---|
| 2659 | * MERGE( 1.0_wp, 0.0_wp, & |
---|
| 2660 | BTEST( wall_flags_0(k,j,i), 3 ) ) |
---|
[1819] | 2661 | ENDDO |
---|
| 2662 | ENDDO |
---|
| 2663 | ENDDO |
---|
| 2664 | ENDIF |
---|
| 2665 | |
---|
| 2666 | |
---|
| 2667 | CASE DEFAULT |
---|
| 2668 | |
---|
| 2669 | WRITE( message_string, * ) 'unknown prognostic variable: ', component |
---|
| 2670 | CALL message( 'wtm_tendencies', 'PA04??', 1, 2, 0, 6, 0 ) |
---|
| 2671 | |
---|
| 2672 | END SELECT |
---|
| 2673 | |
---|
| 2674 | |
---|
| 2675 | END SUBROUTINE wtm_tendencies |
---|
| 2676 | |
---|
| 2677 | |
---|
| 2678 | !------------------------------------------------------------------------------! |
---|
| 2679 | ! Description: |
---|
| 2680 | ! ------------ |
---|
[1839] | 2681 | !> Application of the additional forces generated by the wind turbine on the |
---|
| 2682 | !> flow components (tendency terms) |
---|
| 2683 | !> Call for grid point i,j |
---|
[1819] | 2684 | !------------------------------------------------------------------------------! |
---|
| 2685 | SUBROUTINE wtm_tendencies_ij( i, j, component ) |
---|
| 2686 | |
---|
| 2687 | |
---|
| 2688 | IMPLICIT NONE |
---|
| 2689 | |
---|
[1839] | 2690 | INTEGER(iwp) :: component !< prognostic variable (u,v,w) |
---|
| 2691 | INTEGER(iwp) :: i !< running index |
---|
| 2692 | INTEGER(iwp) :: j !< running index |
---|
| 2693 | INTEGER(iwp) :: k !< running index |
---|
[1819] | 2694 | |
---|
| 2695 | SELECT CASE ( component ) |
---|
| 2696 | |
---|
| 2697 | CASE ( 1 ) |
---|
| 2698 | ! |
---|
| 2699 | !-- Apply the x-component of the force to the u-component of the flow: |
---|
[1839] | 2700 | IF ( simulated_time >= time_turbine_on ) THEN |
---|
[2553] | 2701 | DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) |
---|
[1819] | 2702 | ! |
---|
[1839] | 2703 | !-- Calculate the thrust generated by the nacelle and the tower |
---|
[1912] | 2704 | tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) * & |
---|
[1819] | 2705 | SIGN( u(k,j,i)**2 , u(k,j,i) ) |
---|
[1912] | 2706 | tend_tow_x = 0.5_wp * tow_cd_surf(k,j,i) * & |
---|
[1819] | 2707 | SIGN( u(k,j,i)**2 , u(k,j,i) ) |
---|
[2232] | 2708 | tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i) & |
---|
| 2709 | - tend_nac_x - tend_tow_x ) & |
---|
| 2710 | * MERGE( 1.0_wp, 0.0_wp, & |
---|
[2553] | 2711 | BTEST( wall_flags_0(k,j,i), 1 ) ) |
---|
[1819] | 2712 | ENDDO |
---|
| 2713 | ENDIF |
---|
| 2714 | |
---|
| 2715 | CASE ( 2 ) |
---|
| 2716 | ! |
---|
| 2717 | !-- Apply the y-component of the force to the v-component of the flow: |
---|
[1839] | 2718 | IF ( simulated_time >= time_turbine_on ) THEN |
---|
[2553] | 2719 | DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) |
---|
[1912] | 2720 | tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * & |
---|
[1819] | 2721 | SIGN( v(k,j,i)**2 , v(k,j,i) ) |
---|
[1912] | 2722 | tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * & |
---|
[1819] | 2723 | SIGN( v(k,j,i)**2 , v(k,j,i) ) |
---|
[2232] | 2724 | tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) & |
---|
| 2725 | - tend_nac_y - tend_tow_y ) & |
---|
| 2726 | * MERGE( 1.0_wp, 0.0_wp, & |
---|
| 2727 | BTEST( wall_flags_0(k,j,i), 2 ) ) |
---|
[2553] | 2728 | ENDDO |
---|
[1819] | 2729 | ENDIF |
---|
| 2730 | |
---|
| 2731 | CASE ( 3 ) |
---|
| 2732 | ! |
---|
| 2733 | !-- Apply the z-component of the force to the w-component of the flow: |
---|
[1839] | 2734 | IF ( simulated_time >= time_turbine_on ) THEN |
---|
[2553] | 2735 | DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) |
---|
[2232] | 2736 | tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) & |
---|
| 2737 | * MERGE( 1.0_wp, 0.0_wp, & |
---|
| 2738 | BTEST( wall_flags_0(k,j,i), 3 ) ) |
---|
[1819] | 2739 | ENDDO |
---|
| 2740 | ENDIF |
---|
| 2741 | |
---|
| 2742 | |
---|
| 2743 | CASE DEFAULT |
---|
| 2744 | |
---|
| 2745 | WRITE( message_string, * ) 'unknown prognostic variable: ', component |
---|
| 2746 | CALL message( 'wtm_tendencies', 'PA04??', 1, 2, 0, 6, 0 ) |
---|
| 2747 | |
---|
| 2748 | END SELECT |
---|
| 2749 | |
---|
| 2750 | |
---|
| 2751 | END SUBROUTINE wtm_tendencies_ij |
---|
| 2752 | |
---|
| 2753 | END MODULE wind_turbine_model_mod |
---|