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