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