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