Changeset 4447
- Timestamp:
- Mar 6, 2020 11:05:30 AM (5 years ago)
- Location:
- palm/trunk
- Files:
-
- 2 edited
-
SOURCE/wind_turbine_model_mod.f90 (modified) (98 diffs)
-
TESTS/cases/wind_turbine_model/INPUT/wind_turbine_model_p3d (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wind_turbine_model_mod.f90
r4438 r4447 26 26 ! ----------------- 27 27 ! $Id$ 28 ! renamed wind_turbine_parameters namelist variables 29 ! 30 ! 4438 2020-03-03 20:49:28Z suehring 28 31 ! Bugfix: shifted netcdf preprocessor directive to correct position 29 32 ! … … 179 182 !-- Variables specified in the namelist wind_turbine_par 180 183 181 INTEGER(iwp) :: n airfoils = 8 !< number of airfoils of the used turbine model (for ADM-R and ALM)182 INTEGER(iwp) :: n turbines = 1 !< number of turbines184 INTEGER(iwp) :: n_airfoils = 8 !< number of airfoils of the used turbine model (for ADM-R and ALM) 185 INTEGER(iwp) :: n_turbines = 1 !< number of turbines 183 186 184 187 … … 189 192 REAL(wp), TARGET :: output_values_0d_target !< pointer for 2d output array 190 193 191 LOGICAL :: pitch_control = .FALSE. !< switch for use of pitch controller192 LOGICAL :: speed_control = .FALSE. !< switch for use of speed controller193 LOGICAL :: yaw_control = .FALSE. !< switch for use of yaw controller194 LOGICAL :: t l_cor = .FALSE.!< switch for use of tip loss correct.194 LOGICAL :: pitch_control = .FALSE. !< switch for use of pitch controller 195 LOGICAL :: speed_control = .FALSE. !< switch for use of speed controller 196 LOGICAL :: yaw_control = .FALSE. !< switch for use of yaw controller 197 LOGICAL :: tip_loss_correction = .FALSE. !< switch for use of tip loss correct. 195 198 196 199 LOGICAL :: initial_write_coordinates = .FALSE. … … 200 203 201 204 202 REAL(wp) :: segment_length = 1.0_wp!< length of the segments, the rotor area is divided into203 !< (in tangential direction, as factor of MIN(dx,dy,dz))204 REAL(wp) :: segment_width = 0.5_wp!< width of the segments, the rotor area is divided into205 !< (in radial direction, as factor of MIN(dx,dy,dz))206 REAL(wp) :: time_turbine_on = 0.0_wp!< time at which turbines are started207 REAL(wp) :: tilt = 0.0_wp !< vertical tiltof the rotor [degree] ( positive = backwards )208 209 210 REAL(wp), DIMENSION(1:100) :: dtow = 0.0_wp!< tower diameter [m]211 REAL(wp), DIMENSION(1:100), TARGET :: omega_rot = 0.9_wp!< inital or constant rotor speed [rad/s]212 REAL(wp), DIMENSION(1:100) :: phi_yaw = 0.0_wp!< yaw angle [degree] ( clockwise, 0 = facing west )213 REAL(wp), DIMENSION(1:100) :: pitch_a dd = 0.0_wp!< constant pitch angle214 REAL(wp), DIMENSION(1:100) :: rcx = 9999999.9_wp!< position of hub in x-direction215 REAL(wp), DIMENSION(1:100) :: rcy = 9999999.9_wp!< position of hub in y-direction216 REAL(wp), DIMENSION(1:100) :: rcz = 9999999.9_wp!< position of hub in z-direction217 REAL(wp), DIMENSION(1:100) :: rnac = 0.0_wp!< nacelle diameter [m]218 REAL(wp), DIMENSION(1:100) :: r r = 63.0_wp!< rotor radius [m]219 ! REAL(wp), DIMENSION(1:100) :: turb_cd_nacelle = 0.85_wp!< drag coefficient for nacelle220 REAL(wp), DIMENSION(1:100) :: t urb_cd_tower = 1.2_wp!< drag coefficient for tower205 REAL(wp) :: segment_length_tangential = 1.0_wp !< length of the segments, the rotor area is divided into 206 !< (in tangential direction, as factor of MIN(dx,dy,dz)) 207 REAL(wp) :: segment_width_radial = 0.5_wp !< width of the segments, the rotor area is divided into 208 !< (in radial direction, as factor of MIN(dx,dy,dz)) 209 REAL(wp) :: time_turbine_on = 0.0_wp !< time at which turbines are started 210 REAL(wp) :: tilt_angle = 0.0_wp !< vertical tilt_angle of the rotor [degree] ( positive = backwards ) 211 212 213 REAL(wp), DIMENSION(1:100) :: tower_diameter = 0.0_wp !< tower diameter [m] 214 REAL(wp), DIMENSION(1:100), TARGET :: rotor_speed = 0.9_wp !< inital or constant rotor speed [rad/s] 215 REAL(wp), DIMENSION(1:100) :: yaw_angle = 0.0_wp !< yaw angle [degree] ( clockwise, 0 = facing west ) 216 REAL(wp), DIMENSION(1:100) :: pitch_angle = 0.0_wp !< constant pitch angle 217 REAL(wp), DIMENSION(1:100) :: hub_x = 9999999.9_wp !< position of hub in x-direction 218 REAL(wp), DIMENSION(1:100) :: hub_y = 9999999.9_wp !< position of hub in y-direction 219 REAL(wp), DIMENSION(1:100) :: hub_z = 9999999.9_wp !< position of hub in z-direction 220 REAL(wp), DIMENSION(1:100) :: nacelle_radius = 0.0_wp !< nacelle diameter [m] 221 REAL(wp), DIMENSION(1:100) :: rotor_radius = 63.0_wp !< rotor radius [m] 222 ! REAL(wp), DIMENSION(1:100) :: nacelle_cd = 0.85_wp !< drag coefficient for nacelle 223 REAL(wp), DIMENSION(1:100) :: tower_cd = 1.2_wp !< drag coefficient for tower 221 224 222 225 ! … … 224 227 !-- Default values are from the NREL 5MW research turbine (Jonkman, 2008) 225 228 226 REAL(wp) :: rated_power= 5296610.0_wp !< rated turbine power [W]227 REAL(wp) :: gear_ratio = 97.0_wp !< Gear ratio from rotor to generator228 REAL(wp) :: inertia_rot= 34784179.0_wp !< Inertia of the rotor [kg*m2]229 REAL(wp) :: inertia_gen= 534.116_wp !< Inertia of the generator [kg*m2]230 REAL(wp) :: gen _eff= 0.944_wp !< Electric efficiency of the generator231 REAL(wp) :: gear_eff = 1.0_wp !< Loss between rotor and generator232 REAL(wp) :: air_dens = 1.225_wp !< Air density to convert to W [kg/m3]233 REAL(wp) :: rated_genspeed= 121.6805_wp !< Rated generator speed [rad/s]234 REAL(wp) :: max_torque_gen= 47402.91_wp !< Maximum of the generator torque [Nm]235 REAL(wp) :: slope2= 2.332287_wp !< Slope constant for region 2236 REAL(wp) :: min_reg2= 91.21091_wp !< Lower generator speed boundary of region 2 [rad/s]237 REAL(wp) :: min_reg15= 70.16224_wp !< Lower generator speed boundary of region 1.5 [rad/s]238 REAL(wp) :: max_trq_rate= 15000.0_wp !< Max generator torque increase [Nm/s]239 REAL(wp) :: pitch_rate = 8.0_wp !< Max pitch rate [degree/s]229 REAL(wp) :: generator_power_rated = 5296610.0_wp !< rated turbine power [W] 230 REAL(wp) :: gear_ratio = 97.0_wp !< Gear ratio from rotor to generator 231 REAL(wp) :: rotor_inertia = 34784179.0_wp !< Inertia of the rotor [kg*m2] 232 REAL(wp) :: generator_inertia = 534.116_wp !< Inertia of the generator [kg*m2] 233 REAL(wp) :: generator_efficiency = 0.944_wp !< Electric efficiency of the generator 234 REAL(wp) :: gear_efficiency = 1.0_wp !< Loss between rotor and generator 235 REAL(wp) :: air_density = 1.225_wp !< Air density to convert to W [kg/m3] 236 REAL(wp) :: generator_speed_rated = 121.6805_wp !< Rated generator speed [rad/s] 237 REAL(wp) :: generator_torque_max = 47402.91_wp !< Maximum of the generator torque [Nm] 238 REAL(wp) :: region_2_slope = 2.332287_wp !< Slope constant for region 2 239 REAL(wp) :: region_2_min = 91.21091_wp !< Lower generator speed boundary of region 2 [rad/s] 240 REAL(wp) :: region_15_min = 70.16224_wp !< Lower generator speed boundary of region 1.5 [rad/s] 241 REAL(wp) :: generator_torque_rate_max = 15000.0_wp !< Max generator torque increase [Nm/s] 242 REAL(wp) :: pitch_rate = 8.0_wp !< Max pitch rate [degree/s] 240 243 241 244 … … 243 246 !-- Variables specified in the namelist for yaw control 244 247 245 REAL(wp) :: yaw_speed = 0.005236_wp!< speed of the yaw actuator [rad/s]246 REAL(wp) :: max_miss= 0.08726_wp !< maximum tolerated yaw missalignment [rad]247 REAL(wp) :: min_miss= 0.008726_wp !< minimum yaw missalignment for which the actuator stops [rad]248 REAL(wp) :: yaw_speed = 0.005236_wp !< speed of the yaw actuator [rad/s] 249 REAL(wp) :: yaw_misalignment_max = 0.08726_wp !< maximum tolerated yaw missalignment [rad] 250 REAL(wp) :: yaw_misalignment_min = 0.008726_wp !< minimum yaw missalignment for which the actuator stops [rad] 248 251 249 252 ! … … 392 395 393 396 REAL(wp) :: Fcorner !< corner freq for the controller low pass filter 394 REAL(wp) :: min_reg25!< min region 2.5397 REAL(wp) :: region_25_min !< min region 2.5 395 398 REAL(wp) :: om_rate !< rotor speed change 396 399 REAL(wp) :: slope15 !< slope in region 1.5 397 REAL(wp) :: slope25!< slope in region 2.5400 REAL(wp) :: region_25_slope !< slope in region 2.5 398 401 REAL(wp) :: trq_rate !< torque change 399 402 REAL(wp) :: vs_sysp !< 400 403 REAL(wp) :: lp_coeff !< coeff for the controller low pass filter 401 404 402 REAL(wp), DIMENSION(100) :: omega_rot_l = 0.0_wp !< local rot speed [rad/s]405 REAL(wp), DIMENSION(100) :: rotor_speed_l = 0.0_wp !< local rot speed [rad/s] 403 406 404 407 ! … … 406 409 407 410 REAL(wp), DIMENSION(:) , ALLOCATABLE :: yawdir !< direction to yaw 408 REAL(wp), DIMENSION(:) , ALLOCATABLE :: phi_yaw_l!< local (cpu) yaw angle411 REAL(wp), DIMENSION(:) , ALLOCATABLE :: yaw_angle_l !< local (cpu) yaw angle 409 412 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wd30_l !< local (cpu) long running avg of the wd 410 413 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wd2_l !< local (cpu) short running avg of the wd … … 421 424 ! 422 425 !-- Variables that have to be saved in the binary file for restarts 423 REAL(wp), DIMENSION(1:100) :: pitch_a dd_old = 0.0_wp !< old constant pitch angle424 REAL(wp), DIMENSION(1:100) :: omega_gen= 0.0_wp !< curr. generator speed425 REAL(wp), DIMENSION(1:100) :: omega_gen_f= 0.0_wp !< filtered generator speed426 REAL(wp), DIMENSION(1:100) :: omega_gen_old= 0.0_wp !< last generator speed427 REAL(wp), DIMENSION(1:100) :: omega_gen_f_old= 0.0_wp !< last filtered generator speed426 REAL(wp), DIMENSION(1:100) :: pitch_angle_old = 0.0_wp !< old constant pitch angle 427 REAL(wp), DIMENSION(1:100) :: generator_speed = 0.0_wp !< curr. generator speed 428 REAL(wp), DIMENSION(1:100) :: generator_speed_f = 0.0_wp !< filtered generator speed 429 REAL(wp), DIMENSION(1:100) :: generator_speed_old = 0.0_wp !< last generator speed 430 REAL(wp), DIMENSION(1:100) :: generator_speed_f_old = 0.0_wp !< last filtered generator speed 428 431 REAL(wp), DIMENSION(1:100) :: torque_gen = 0.0_wp !< generator torque 429 432 REAL(wp), DIMENSION(1:100) :: torque_gen_old = 0.0_wp !< last generator torque … … 503 506 504 507 NAMELIST /wind_turbine_parameters/ & 505 air_dens, dtow, dt_data_output_wtm, gear_eff,& 506 gear_ratio, & 507 gen_eff, inertia_gen, inertia_rot, max_miss, & 508 max_torque_gen, max_trq_rate, min_miss, & 509 min_reg15, min_reg2, nairfoils, nturbines, & 510 omega_rot, phi_yaw, pitch_add, pitch_control,& 511 rated_genspeed, rated_power, rcx, rcy, rcz, & 512 rnac, rr, segment_length, segment_width, & 513 slope2, speed_control, tilt, time_turbine_on,& 514 turb_cd_tower, pitch_rate, & 515 yaw_control, yaw_speed, tl_cor 516 ! , turb_cd_nacelle 508 air_density, tower_diameter, dt_data_output_wtm, & 509 gear_efficiency, gear_ratio, generator_efficiency, & 510 generator_inertia, rotor_inertia, yaw_misalignment_max, & 511 generator_torque_max, generator_torque_rate_max, yaw_misalignment_min, & 512 region_15_min, region_2_min, n_airfoils, n_turbines, & 513 rotor_speed, yaw_angle, pitch_angle, pitch_control, & 514 generator_speed_rated, generator_power_rated, hub_x, hub_y, hub_z,& 515 nacelle_radius, rotor_radius, segment_length_tangential, segment_width_radial, & 516 region_2_slope, speed_control, tilt_angle, time_turbine_on, & 517 tower_cd, pitch_rate, & 518 yaw_control, yaw_speed, tip_loss_correction 519 ! , nacelle_cd 520 517 521 ! 518 522 !-- Try to find wind turbine model package … … 553 557 554 558 555 CALL wrd_write_string( ' omega_gen' )556 WRITE ( 14 ) omega_gen557 558 CALL wrd_write_string( ' omega_gen_f' )559 WRITE ( 14 ) omega_gen_f560 561 CALL wrd_write_string( ' omega_gen_f_old' )562 WRITE ( 14 ) omega_gen_f_old563 564 CALL wrd_write_string( ' omega_gen_old' )565 WRITE ( 14 ) omega_gen_old566 567 CALL wrd_write_string( ' omega_rot' )568 WRITE ( 14 ) omega_rot569 570 CALL wrd_write_string( ' phi_yaw' )571 WRITE ( 14 ) phi_yaw572 573 CALL wrd_write_string( 'pitch_a dd' )574 WRITE ( 14 ) pitch_a dd575 576 CALL wrd_write_string( 'pitch_a dd_old' )577 WRITE ( 14 ) pitch_a dd_old559 CALL wrd_write_string( 'generator_speed' ) 560 WRITE ( 14 ) generator_speed 561 562 CALL wrd_write_string( 'generator_speed_f' ) 563 WRITE ( 14 ) generator_speed_f 564 565 CALL wrd_write_string( 'generator_speed_f_old' ) 566 WRITE ( 14 ) generator_speed_f_old 567 568 CALL wrd_write_string( 'generator_speed_old' ) 569 WRITE ( 14 ) generator_speed_old 570 571 CALL wrd_write_string( 'rotor_speed' ) 572 WRITE ( 14 ) rotor_speed 573 574 CALL wrd_write_string( 'yaw_angle' ) 575 WRITE ( 14 ) yaw_angle 576 577 CALL wrd_write_string( 'pitch_angle' ) 578 WRITE ( 14 ) pitch_angle 579 580 CALL wrd_write_string( 'pitch_angle_old' ) 581 WRITE ( 14 ) pitch_angle_old 578 582 579 583 CALL wrd_write_string( 'torque_gen' ) … … 609 613 SELECT CASE ( restart_string(1:length) ) 610 614 611 CASE ( ' omega_gen' )612 READ ( 13 ) omega_gen613 CASE ( ' omega_gen_f' )614 READ ( 13 ) omega_gen_f615 CASE ( ' omega_gen_f_old' )616 READ ( 13 ) omega_gen_f_old617 CASE ( ' omega_gen_old' )618 READ ( 13 ) omega_gen_old619 CASE ( ' omega_rot' )620 READ ( 13 ) omega_rot621 CASE ( ' phi_yaw' )622 READ ( 13 ) phi_yaw623 CASE ( 'pitch_a dd' )624 READ ( 13 ) pitch_a dd625 CASE ( 'pitch_a dd_old' )626 READ ( 13 ) pitch_a dd_old615 CASE ( 'generator_speed' ) 616 READ ( 13 ) generator_speed 617 CASE ( 'generator_speed_f' ) 618 READ ( 13 ) generator_speed_f 619 CASE ( 'generator_speed_f_old' ) 620 READ ( 13 ) generator_speed_f_old 621 CASE ( 'generator_speed_old' ) 622 READ ( 13 ) generator_speed_old 623 CASE ( 'rotor_speed' ) 624 READ ( 13 ) rotor_speed 625 CASE ( 'yaw_angle' ) 626 READ ( 13 ) yaw_angle 627 CASE ( 'pitch_angle' ) 628 READ ( 13 ) pitch_angle 629 CASE ( 'pitch_angle_old' ) 630 READ ( 13 ) pitch_angle_old 627 631 CASE ( 'torque_gen' ) 628 632 READ ( 13 ) torque_gen … … 656 660 ENDIF 657 661 658 IF ( ANY( omega_rot(1:nturbines) < 0.0 ) ) THEN659 message_string = ' omega_rot < 0.0, Please set omega_rotto ' // &662 IF ( ANY( rotor_speed(1:n_turbines) < 0.0 ) ) THEN 663 message_string = 'rotor_speed < 0.0, Please set rotor_speed to ' // & 660 664 'a value equal or larger than zero' 661 665 CALL message( 'wtm_check_parameters', 'PA0462', 1, 2, 0, 6, 0 ) … … 663 667 664 668 665 IF ( ANY( rcx(1:nturbines) == 9999999.9_wp ) .OR.&666 ANY( rcy(1:nturbines) == 9999999.9_wp ) .OR.&667 ANY( rcz(1:nturbines) == 9999999.9_wp ) ) THEN669 IF ( ANY( hub_x(1:n_turbines) == 9999999.9_wp ) .OR. & 670 ANY( hub_y(1:n_turbines) == 9999999.9_wp ) .OR. & 671 ANY( hub_z(1:n_turbines) == 9999999.9_wp ) ) THEN 668 672 669 message_string = ' rcx, rcy, rcz '// &673 message_string = 'hub_x, hub_y, hub_z ' // & 670 674 'have to be given for each turbine.' 671 675 CALL message( 'wtm_check_parameters', 'PA0463', 1, 2, 0, 6, 0 ) … … 707 711 ! 708 712 !-- Input of all wtm parameters 709 IF ( check_existence( vars_pids, ' dtow' ) ) THEN710 CALL get_variable( pids_id, ' dtow', dtow(1:nturbines))713 IF ( check_existence( vars_pids, 'tower_diameter' ) ) THEN 714 CALL get_variable( pids_id, 'tower_diameter', tower_diameter(1:n_turbines)) 711 715 ENDIF 712 716 713 IF ( check_existence( vars_pids, ' omega_rot' ) ) THEN714 CALL get_variable( pids_id, ' omega_rot', omega_rot(1:nturbines))717 IF ( check_existence( vars_pids, 'rotor_speed' ) ) THEN 718 CALL get_variable( pids_id, 'rotor_speed', rotor_speed(1:n_turbines)) 715 719 ENDIF 716 720 717 IF ( check_existence( vars_pids, 'pitch_a dd' ) ) THEN718 CALL get_variable( pids_id, 'pitch_a dd', pitch_add(1:nturbines))721 IF ( check_existence( vars_pids, 'pitch_angle' ) ) THEN 722 CALL get_variable( pids_id, 'pitch_angle', pitch_angle(1:n_turbines)) 719 723 ENDIF 720 724 721 IF ( check_existence( vars_pids, ' phi_yaw' ) ) THEN722 CALL get_variable( pids_id, ' phi_yaw', phi_yaw(1:nturbines))725 IF ( check_existence( vars_pids, 'yaw_angle' ) ) THEN 726 CALL get_variable( pids_id, 'yaw_angle', yaw_angle(1:n_turbines)) 723 727 ENDIF 724 728 725 IF ( check_existence( vars_pids, ' rcx' ) ) THEN726 CALL get_variable( pids_id, ' rcx', rcx(1:nturbines))729 IF ( check_existence( vars_pids, 'hub_x' ) ) THEN 730 CALL get_variable( pids_id, 'hub_x', hub_x(1:n_turbines)) 727 731 ENDIF 728 732 729 IF ( check_existence( vars_pids, ' rcy' ) ) THEN730 CALL get_variable( pids_id, ' rcy', rcy(1:nturbines))733 IF ( check_existence( vars_pids, 'hub_y' ) ) THEN 734 CALL get_variable( pids_id, 'hub_y', hub_y(1:n_turbines)) 731 735 ENDIF 732 736 733 IF ( check_existence( vars_pids, ' rcz' ) ) THEN734 CALL get_variable( pids_id, ' rcz', rcz(1:nturbines))737 IF ( check_existence( vars_pids, 'hub_z' ) ) THEN 738 CALL get_variable( pids_id, 'hub_z', hub_z(1:n_turbines)) 735 739 ENDIF 736 740 737 IF ( check_existence( vars_pids, ' rnac' ) ) THEN738 CALL get_variable( pids_id, ' rnac', rnac(1:nturbines))741 IF ( check_existence( vars_pids, 'nacelle_radius' ) ) THEN 742 CALL get_variable( pids_id, 'nacelle_radius', nacelle_radius(1:n_turbines)) 739 743 ENDIF 740 744 741 IF ( check_existence( vars_pids, 'r r' ) ) THEN742 CALL get_variable( pids_id, 'r r', rr(1:nturbines))745 IF ( check_existence( vars_pids, 'rotor_radius' ) ) THEN 746 CALL get_variable( pids_id, 'rotor_radius', rotor_radius(1:n_turbines)) 743 747 ENDIF 744 748 ! 745 ! IF ( check_existence( vars_pids, ' turb_cd_nacelle' ) ) THEN746 ! CALL get_variable( pids_id, ' turb_cd_nacelle', turb_cd_nacelle(1:nturbines))749 ! IF ( check_existence( vars_pids, 'nacelle_cd' ) ) THEN 750 ! CALL get_variable( pids_id, 'nacelle_cd', nacelle_cd(1:n_turbines)) 747 751 ! ENDIF 748 752 749 IF ( check_existence( vars_pids, 't urb_cd_tower' ) ) THEN750 CALL get_variable( pids_id, 't urb_cd_tower', turb_cd_tower(1:nturbines))753 IF ( check_existence( vars_pids, 'tower_cd' ) ) THEN 754 CALL get_variable( pids_id, 'tower_cd', tower_cd(1:n_turbines)) 751 755 ENDIF 752 756 ! … … 761 765 !-- the maximum possible numbers of rings and segments have to be calculated: 762 766 763 ALLOCATE( nrings(1:n turbines) )764 ALLOCATE( delta_r(1:n turbines) )767 ALLOCATE( nrings(1:n_turbines) ) 768 ALLOCATE( delta_r(1:n_turbines) ) 765 769 766 770 nrings(:) = 0 … … 769 773 ! 770 774 !-- Thickness (radial) of each ring and length (tangential) of each segment: 771 delta_r_factor = segment_width 772 delta_t_factor = segment_length 775 delta_r_factor = segment_width_radial 776 delta_t_factor = segment_length_tangential 773 777 delta_r_init = delta_r_factor * MIN( dx, dy, dz(1)) 774 778 775 DO inot = 1, n turbines779 DO inot = 1, n_turbines 776 780 ! 777 781 !-- Determine number of rings: 778 nrings(inot) = NINT( r r(inot) / delta_r_init )779 780 delta_r(inot) = r r(inot) / nrings(inot)782 nrings(inot) = NINT( rotor_radius(inot) / delta_r_init ) 783 784 delta_r(inot) = rotor_radius(inot) / nrings(inot) 781 785 782 786 ENDDO … … 784 788 nrings_max = MAXVAL(nrings) 785 789 786 ALLOCATE( nsegs(1:nrings_max,1:n turbines) )787 ALLOCATE( nsegs_total(1:n turbines) )790 ALLOCATE( nsegs(1:nrings_max,1:n_turbines) ) 791 ALLOCATE( nsegs_total(1:n_turbines) ) 788 792 789 793 nsegs(:,:) = 0 … … 791 795 792 796 793 DO inot = 1, n turbines797 DO inot = 1, n_turbines 794 798 DO ring = 1, nrings(inot) 795 799 ! … … 824 828 ! 825 829 !-- Allocate 1D arrays (dimension = number of turbines) 826 ALLOCATE( i_hub(1:n turbines) )827 ALLOCATE( i_smear(1:n turbines) )828 ALLOCATE( j_hub(1:n turbines) )829 ALLOCATE( j_smear(1:n turbines) )830 ALLOCATE( k_hub(1:n turbines) )831 ALLOCATE( k_smear(1:n turbines) )832 ALLOCATE( torque_total(1:n turbines) )833 ALLOCATE( thrust_rotor(1:n turbines) )830 ALLOCATE( i_hub(1:n_turbines) ) 831 ALLOCATE( i_smear(1:n_turbines) ) 832 ALLOCATE( j_hub(1:n_turbines) ) 833 ALLOCATE( j_smear(1:n_turbines) ) 834 ALLOCATE( k_hub(1:n_turbines) ) 835 ALLOCATE( k_smear(1:n_turbines) ) 836 ALLOCATE( torque_total(1:n_turbines) ) 837 ALLOCATE( thrust_rotor(1:n_turbines) ) 834 838 835 839 ! 836 840 !-- Allocation of the 1D arrays for yaw control 837 ALLOCATE( yawdir(1:n turbines) )838 ALLOCATE( u_inflow(1:n turbines) )839 ALLOCATE( wdir(1:n turbines) )840 ALLOCATE( u_inflow_l(1:n turbines) )841 ALLOCATE( wdir_l(1:n turbines) )842 ALLOCATE( phi_yaw_l(1:nturbines) )841 ALLOCATE( yawdir(1:n_turbines) ) 842 ALLOCATE( u_inflow(1:n_turbines) ) 843 ALLOCATE( wdir(1:n_turbines) ) 844 ALLOCATE( u_inflow_l(1:n_turbines) ) 845 ALLOCATE( wdir_l(1:n_turbines) ) 846 ALLOCATE( yaw_angle_l(1:n_turbines) ) 843 847 844 848 ! … … 866 870 ! 867 871 !-- Allocate additional 2D arrays 868 ALLOCATE( rotx(1:n turbines,1:3) )869 ALLOCATE( roty(1:n turbines,1:3) )870 ALLOCATE( rotz(1:n turbines,1:3) )872 ALLOCATE( rotx(1:n_turbines,1:3) ) 873 ALLOCATE( roty(1:n_turbines,1:3) ) 874 ALLOCATE( rotz(1:n_turbines,1:3) ) 871 875 872 876 ! … … 883 887 ! 884 888 !-- Allocate additional 3D arrays 885 ALLOCATE( u_int(1:n turbines,1:nrings_max,1:nsegs_max) )886 ALLOCATE( u_int_1_l(1:n turbines,1:nrings_max,1:nsegs_max) )887 ALLOCATE( v_int(1:n turbines,1:nrings_max,1:nsegs_max) )888 ALLOCATE( v_int_1_l(1:n turbines,1:nrings_max,1:nsegs_max) )889 ALLOCATE( w_int(1:n turbines,1:nrings_max,1:nsegs_max) )890 ALLOCATE( w_int_1_l(1:n turbines,1:nrings_max,1:nsegs_max) )889 ALLOCATE( u_int(1:n_turbines,1:nrings_max,1:nsegs_max) ) 890 ALLOCATE( u_int_1_l(1:n_turbines,1:nrings_max,1:nsegs_max) ) 891 ALLOCATE( v_int(1:n_turbines,1:nrings_max,1:nsegs_max) ) 892 ALLOCATE( v_int_1_l(1:n_turbines,1:nrings_max,1:nsegs_max) ) 893 ALLOCATE( w_int(1:n_turbines,1:nrings_max,1:nsegs_max) ) 894 ALLOCATE( w_int_1_l(1:n_turbines,1:nrings_max,1:nsegs_max) ) 891 895 892 896 ! … … 903 907 904 908 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 905 omega_gen(:) = 0.0_wp906 omega_gen_old(:) = 0.0_wp907 omega_gen_f(:) = 0.0_wp908 omega_gen_f_old(:) = 0.0_wp909 pitch_a dd_old(:) = 0.0_wp909 generator_speed(:) = 0.0_wp 910 generator_speed_old(:) = 0.0_wp 911 generator_speed_f(:) = 0.0_wp 912 generator_speed_f_old(:) = 0.0_wp 913 pitch_angle_old(:) = 0.0_wp 910 914 torque_gen(:) = 0.0_wp 911 915 torque_gen_old(:) = 0.0_wp … … 917 921 u_inflow(:) = 0.0_wp 918 922 u_inflow_l(:) = 0.0_wp 919 phi_yaw_l(:)= 0.0_wp923 yaw_angle_l(:) = 0.0_wp 920 924 921 925 ! … … 1001 1005 IF ( debug_output ) CALL debug_message( 'wtm_init', 'start' ) 1002 1006 1003 ALLOCATE( index_nacb(1:n turbines) )1004 ALLOCATE( index_nacl(1:n turbines) )1005 ALLOCATE( index_nacr(1:n turbines) )1006 ALLOCATE( index_nact(1:n turbines) )1007 ALLOCATE( index_nacb(1:n_turbines) ) 1008 ALLOCATE( index_nacl(1:n_turbines) ) 1009 ALLOCATE( index_nacr(1:n_turbines) ) 1010 ALLOCATE( index_nact(1:n_turbines) ) 1007 1011 1008 1012 ! … … 1030 1034 !-- the area where the wtm is active. ABS (...) is required because the 1031 1035 !-- default value of dz_stretch_level_start is -9999999.9_wp (negative). 1032 IF ( ABS( dz_stretch_level_start(1) ) <= MAXVAL( rcz(1:nturbines)) +&1033 MAXVAL(r r(1:nturbines)) +&1036 IF ( ABS( dz_stretch_level_start(1) ) <= MAXVAL(hub_z(1:n_turbines)) + & 1037 MAXVAL(rotor_radius(1:n_turbines)) + & 1034 1038 eps_min) THEN 1035 1039 WRITE( message_string, * ) 'The lowest level where vertical ', & 1036 1040 'stretching is applied &have to be ', & 1037 'greater than ',MAXVAL( rcz(1:nturbines)) +&1038 MAXVAL(rr(1:nturbines)) + eps_min1041 'greater than ',MAXVAL(hub_z(1:n_turbines)) & 1042 + MAXVAL(rotor_radius(1:n_turbines)) + eps_min 1039 1043 CALL message( 'wtm_init', 'PA0484', 1, 2, 0, 6, 0 ) 1040 1044 ENDIF … … 1052 1056 1053 1057 !-- Change tilt angle to rad: 1054 tilt = tilt* pi / 180.0_wp1058 tilt_angle = tilt_angle * pi / 180.0_wp 1055 1059 1056 1060 ! 1057 1061 !-- Change yaw angle to rad: 1058 1062 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 1059 phi_yaw(:) = phi_yaw(:) * pi / 180.0_wp1063 yaw_angle(:) = yaw_angle(:) * pi / 180.0_wp 1060 1064 ENDIF 1061 1065 1062 1066 1063 DO inot = 1, n turbines1067 DO inot = 1, n_turbines 1064 1068 ! 1065 1069 !-- Rotate the rotor coordinates in case yaw and tilt are defined … … 1068 1072 ! 1069 1073 !-- Determine the indices of the hub height 1070 i_hub(inot) = INT( rcx(inot) / dx )1071 j_hub(inot) = INT( ( rcy(inot) + 0.5_wp * dy ) / dy )1072 k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz(1) ) / dz(1) )1074 i_hub(inot) = INT( hub_x(inot) / dx ) 1075 j_hub(inot) = INT( ( hub_y(inot) + 0.5_wp * dy ) / dy ) 1076 k_hub(inot) = INT( ( hub_z(inot) + 0.5_wp * dz(1) ) / dz(1) ) 1073 1077 1074 1078 ! … … 1077 1081 !-- eps_min, the smearing area can be further limited and regarded as a 1078 1082 !-- function of eps_min: 1079 i_smear(inot) = CEILING( ( r r(inot) + eps_min ) / dx )1080 j_smear(inot) = CEILING( ( r r(inot) + eps_min ) / dy )1081 k_smear(inot) = CEILING( ( r r(inot) + eps_min ) / dz(1) )1083 i_smear(inot) = CEILING( ( rotor_radius(inot) + eps_min ) / dx ) 1084 j_smear(inot) = CEILING( ( rotor_radius(inot) + eps_min ) / dy ) 1085 k_smear(inot) = CEILING( ( rotor_radius(inot) + eps_min ) / dz(1) ) 1082 1086 1083 1087 ENDDO … … 1085 1089 ! 1086 1090 !-- Call the wtm_init_speed_control subroutine and calculate the local 1087 !-- omega_rotfor the respective processor.1091 !-- rotor_speed for the respective processor. 1088 1092 IF ( speed_control) THEN 1089 1093 … … 1092 1096 IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN 1093 1097 1094 DO inot = 1, n turbines1098 DO inot = 1, n_turbines 1095 1099 1096 1100 IF ( nxl > i_hub(inot) ) THEN 1097 1101 torque_gen(inot) = 0.0_wp 1098 omega_gen_f(inot) = 0.0_wp1099 omega_rot_l(inot) = 0.0_wp1102 generator_speed_f(inot) = 0.0_wp 1103 rotor_speed_l(inot) = 0.0_wp 1100 1104 ENDIF 1101 1105 1102 1106 IF ( nxr < i_hub(inot) ) THEN 1103 1107 torque_gen(inot) = 0.0_wp 1104 omega_gen_f(inot) = 0.0_wp1105 omega_rot_l(inot) = 0.0_wp1108 generator_speed_f(inot) = 0.0_wp 1109 rotor_speed_l(inot) = 0.0_wp 1106 1110 ENDIF 1107 1111 1108 1112 IF ( nys > j_hub(inot) ) THEN 1109 1113 torque_gen(inot) = 0.0_wp 1110 omega_gen_f(inot) = 0.0_wp1111 omega_rot_l(inot) = 0.0_wp1114 generator_speed_f(inot) = 0.0_wp 1115 rotor_speed_l(inot) = 0.0_wp 1112 1116 ENDIF 1113 1117 1114 1118 IF ( nyn < j_hub(inot) ) THEN 1115 1119 torque_gen(inot) = 0.0_wp 1116 omega_gen_f(inot) = 0.0_wp1117 omega_rot_l(inot) = 0.0_wp1120 generator_speed_f(inot) = 0.0_wp 1121 rotor_speed_l(inot) = 0.0_wp 1118 1122 ENDIF 1119 1123 … … 1121 1125 IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) ) THEN 1122 1126 1123 omega_rot_l(inot) = omega_gen(inot) / gear_ratio1127 rotor_speed_l(inot) = generator_speed(inot) / gear_ratio 1124 1128 1125 1129 ENDIF … … 1151 1155 1152 1156 1153 DO inot = 1, n turbines ! loop over number of turbines1157 DO inot = 1, n_turbines ! loop over number of turbines 1154 1158 ! 1155 1159 !-- Determine the grid index (u-grid) that corresponds to the location of … … 1161 1165 !-- Determine the left and the right edge of the nacelle (corresponding 1162 1166 !-- grid point indices): 1163 index_nacl(inot) = INT( ( rcy(inot) - rnac(inot) + 0.5_wp * dy ) / dy )1164 index_nacr(inot) = INT( ( rcy(inot) + rnac(inot) + 0.5_wp * dy ) / dy )1167 index_nacl(inot) = INT( ( hub_y(inot) - nacelle_radius(inot) + 0.5_wp * dy ) / dy ) 1168 index_nacr(inot) = INT( ( hub_y(inot) + nacelle_radius(inot) + 0.5_wp * dy ) / dy ) 1165 1169 ! 1166 1170 !-- Determine the bottom and the top edge of the nacelle (corresponding … … 1169 1173 !-- surface. All points between z=0 and z=dz/s would already be contained 1170 1174 !-- in grid box 1. 1171 index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz(1) ) + 11172 index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz(1) ) + 11175 index_nacb(inot) = INT( ( hub_z(inot) - nacelle_radius(inot) ) / dz(1) ) + 1 1176 index_nact(inot) = INT( ( hub_z(inot) + nacelle_radius(inot) ) / dz(1) ) + 1 1173 1177 1174 1178 ! 1175 1179 !-- Determine the indices of the grid boxes containing the left and 1176 1180 !-- the right boundaries of the tower: 1177 tower_n = ( rcy(inot) + 0.5_wp * dtow(inot) - 0.5_wp * dy ) / dy1178 tower_s = ( rcy(inot) - 0.5_wp * dtow(inot) - 0.5_wp * dy ) / dy1181 tower_n = ( hub_y(inot) + 0.5_wp * tower_diameter(inot) - 0.5_wp * dy ) / dy 1182 tower_s = ( hub_y(inot) - 0.5_wp * tower_diameter(inot) - 0.5_wp * dy ) / dy 1179 1183 1180 1184 ! … … 1195 1199 !-- leftmost and rightmost grid box: 1196 1200 IF ( j == tower_s ) THEN 1197 tow_cd_surf(k,j,i) = ( rcz(inot) - &1201 tow_cd_surf(k,j,i) = ( hub_z(inot) - & 1198 1202 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction 1199 1203 ( ( tower_s + 1.0_wp + 0.5_wp ) * dy - & 1200 ( rcy(inot) - 0.5_wp * dtow(inot) ) ) * & ! extension in y-direction1201 t urb_cd_tower(inot)1204 ( hub_y(inot) - 0.5_wp * tower_diameter(inot) ) ) * & ! extension in y-direction 1205 tower_cd(inot) 1202 1206 ELSEIF ( j == tower_n ) THEN 1203 tow_cd_surf(k,j,i) = ( rcz(inot) - &1207 tow_cd_surf(k,j,i) = ( hub_z(inot) - & 1204 1208 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction 1205 ( ( rcy(inot) + 0.5_wp * dtow(inot) ) - &1209 ( ( hub_y(inot) + 0.5_wp * tower_diameter(inot) ) - & 1206 1210 ( tower_n + 0.5_wp ) * dy ) * & ! extension in y-direction 1207 t urb_cd_tower(inot)1211 tower_cd(inot) 1208 1212 ! 1209 1213 !-- grid boxes inbetween 1210 1214 !-- (where tow_cd_surf = grid box area): 1211 1215 ELSE 1212 tow_cd_surf(k,j,i) = ( rcz(inot) - &1216 tow_cd_surf(k,j,i) = ( hub_z(inot) - & 1213 1217 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& 1214 dy * t urb_cd_tower(inot)1218 dy * tower_cd(inot) 1215 1219 ENDIF 1216 1220 ! 1217 1221 !-- tower lies completely within one grid box: 1218 1222 ELSE 1219 tow_cd_surf(k,j,i) = ( rcz(inot) - ( k_hub(inot) * &1223 tow_cd_surf(k,j,i) = ( hub_z(inot) - ( k_hub(inot) * & 1220 1224 dz(1) - 0.5_wp * dz(1) ) ) * & 1221 dtow(inot) * turb_cd_tower(inot)1225 tower_diameter(inot) * tower_cd(inot) 1222 1226 ENDIF 1223 1227 ! … … 1232 1236 tow_cd_surf(k,j,i) = dz(1) * ( & 1233 1237 ( tower_s + 1 + 0.5_wp ) * dy - & 1234 ( rcy(inot) - 0.5_wp * dtow(inot) ) &1235 ) * t urb_cd_tower(inot)1238 ( hub_y(inot) - 0.5_wp * tower_diameter(inot) ) & 1239 ) * tower_cd(inot) 1236 1240 ELSEIF ( j == tower_n ) THEN 1237 1241 tow_cd_surf(k,j,i) = dz(1) * ( & 1238 ( rcy(inot) + 0.5_wp * dtow(inot) ) - &1242 ( hub_y(inot) + 0.5_wp * tower_diameter(inot) ) - & 1239 1243 ( tower_n + 0.5_wp ) * dy & 1240 ) * t urb_cd_tower(inot)1244 ) * tower_cd(inot) 1241 1245 ! 1242 1246 !-- grid boxes inbetween … … 1244 1248 ELSE 1245 1249 tow_cd_surf(k,j,i) = dz(1) * dy * & 1246 t urb_cd_tower(inot)1250 tower_cd(inot) 1247 1251 ENDIF 1248 1252 ! 1249 1253 !-- tower lies completely within one grid box: 1250 1254 ELSE 1251 tow_cd_surf(k,j,i) = dz(1) * dtow(inot) * &1252 t urb_cd_tower(inot)1255 tow_cd_surf(k,j,i) = dz(1) * tower_diameter(inot) * & 1256 tower_cd(inot) 1253 1257 ENDIF ! end if larger than grid box 1254 1258 … … 1279 1283 ! DO i_ip = 1, upper_end 1280 1284 ! yvalue = dy_int * ( i_ip - 0.5_wp ) + 0.5_wp * dy !<--- segmentation fault 1281 ! sqrt_arg = rnac(inot)**2 - ( yvalue - rcy(inot) )**2 !<--- segmentation fault1285 ! sqrt_arg = nacelle_radius(inot)**2 - ( yvalue - hub_y(inot) )**2 !<--- segmentation fault 1282 1286 ! IF ( sqrt_arg >= 0.0_wp ) THEN 1283 1287 !! 1284 1288 !!-- bottom intersection point 1285 ! circle_points(1,i_ip) = rcz(inot) - SQRT( sqrt_arg )1289 ! circle_points(1,i_ip) = hub_z(inot) - SQRT( sqrt_arg ) 1286 1290 !! 1287 1291 !!-- top intersection point 1288 ! circle_points(2,i_ip) = rcz(inot) + SQRT( sqrt_arg ) !<--- segmentation fault1292 ! circle_points(2,i_ip) = hub_z(inot) + SQRT( sqrt_arg ) !<--- segmentation fault 1289 1293 ! ELSE 1290 1294 ! circle_points(:,i_ip) = -111111 !<--- segmentation fault … … 1338 1342 ! IF ( ( nxlg <= i ) .AND. ( nxrg >= i ) ) THEN 1339 1343 ! nac_cd_surf(k,j,i) = nac_cd_surf(k,j,i) + & !<--- segmentation fault 1340 ! dy_int * dz_int * turb_cd_nacelle(inot)1344 ! dy_int * dz_int * nacelle_cd(inot) 1341 1345 ! ENDIF 1342 1346 ! ENDDO … … 1382 1386 ! 1383 1387 !-- Define dimensions in output file 1384 ALLOCATE( ndim(1:n turbines) )1385 DO n = 1, n turbines1388 ALLOCATE( ndim(1:n_turbines) ) 1389 DO n = 1, n_turbines 1386 1390 ndim(n) = n 1387 1391 ENDDO … … 1389 1393 dimension_name = 'turbine', & 1390 1394 output_type = 'int32', & 1391 bounds = (/1_iwp, n turbines/),&1395 bounds = (/1_iwp, n_turbines/), & 1392 1396 values_int32 = ndim ) 1393 1397 DEALLOCATE( ndim ) … … 1649 1653 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cl_table !< 1650 1654 1651 ALLOCATE ( read_cl_cd(1:2*n airfoils+1) )1655 ALLOCATE ( read_cl_cd(1:2*n_airfoils+1) ) 1652 1656 1653 1657 ! … … 1718 1722 ENDDO rloop2 1719 1723 1720 ALLOCATE( alpha_attack_tab(1:dlen), turb_cl_table(1:dlen,1:n airfoils),&1721 turb_cd_table(1:dlen,1:n airfoils) )1724 ALLOCATE( alpha_attack_tab(1:dlen), turb_cl_table(1:dlen,1:n_airfoils), & 1725 turb_cd_table(1:dlen,1:n_airfoils) ) 1722 1726 1723 1727 DO jj = 1,dlen+1 … … 1728 1732 READ ( 90, * ) read_cl_cd 1729 1733 alpha_attack_tab(jj) = read_cl_cd(1) 1730 DO ii= 1, n airfoils1734 DO ii= 1, n_airfoils 1731 1735 turb_cl_table(jj,ii) = read_cl_cd(ii*2) 1732 1736 turb_cd_table(jj,ii) = read_cl_cd(ii*2+1) … … 1740 1744 1741 1745 ! 1742 !-- For each possible radial position (resolution: 0.1 m --> 631 values if r r(1)=63m)1746 !-- For each possible radial position (resolution: 0.1 m --> 631 values if rotor_radius(1)=63m) 1743 1747 !-- and each possible angle of attack (resolution: 0.1 degrees --> 3601 values!) 1744 1748 !-- determine the lift and drag coefficient by interpolating between the … … 1751 1755 ALLOCATE( turb_cd_sel2(1:dlenbl) ) 1752 1756 1753 radres = INT( r r(1) * 10.0_wp ) + 1_iwp1757 radres = INT( rotor_radius(1) * 10.0_wp ) + 1_iwp 1754 1758 dlenbl_int = INT( 360.0_wp / accu_cl_cd_tab ) + 1_iwp 1755 1759 … … 1866 1870 !-- Calculation of the rotation matrix for the application of the tilt to 1867 1871 !-- the rotors 1868 rot_eigen_rad(1) = SIN( phi_yaw(inot) )! x-component of the radial eigenvector1869 rot_eigen_rad(2) = COS( phi_yaw(inot) )! y-component of the radial eigenvector1872 rot_eigen_rad(1) = SIN( yaw_angle(inot) ) ! x-component of the radial eigenvector 1873 rot_eigen_rad(2) = COS( yaw_angle(inot) ) ! y-component of the radial eigenvector 1870 1874 rot_eigen_rad(3) = 0.0_wp ! z-component of the radial eigenvector 1871 1875 … … 1874 1878 rot_eigen_azi(3) = 1.0_wp ! z-component of the azimuth eigenvector 1875 1879 1876 rot_eigen_nor(1) = COS( phi_yaw(inot) )! x-component of the normal eigenvector1877 rot_eigen_nor(2) = -SIN( phi_yaw(inot) )! y-component of the normal eigenvector1880 rot_eigen_nor(1) = COS( yaw_angle(inot) ) ! x-component of the normal eigenvector 1881 rot_eigen_nor(2) = -SIN( yaw_angle(inot) ) ! y-component of the normal eigenvector 1878 1882 rot_eigen_nor(3) = 0.0_wp ! z-component of the normal eigenvector 1879 1883 … … 1883 1887 1884 1888 rot_coord_trans(inot,1,1) = rot_eigen_rad(1)**2 * & 1885 ( 1.0_wp - COS( tilt ) ) + COS( tilt)1889 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1886 1890 rot_coord_trans(inot,1,2) = rot_eigen_rad(1) * rot_eigen_rad(2) * & 1887 ( 1.0_wp - COS( tilt ) ) - &1888 rot_eigen_rad(3) * SIN( tilt )1891 ( 1.0_wp - COS( tilt_angle ) ) - & 1892 rot_eigen_rad(3) * SIN( tilt_angle ) 1889 1893 rot_coord_trans(inot,1,3) = rot_eigen_rad(1) * rot_eigen_rad(3) * & 1890 ( 1.0_wp - COS( tilt ) ) + &1891 rot_eigen_rad(2) * SIN( tilt )1894 ( 1.0_wp - COS( tilt_angle ) ) + & 1895 rot_eigen_rad(2) * SIN( tilt_angle ) 1892 1896 rot_coord_trans(inot,2,1) = rot_eigen_rad(2) * rot_eigen_rad(1) * & 1893 ( 1.0_wp - COS( tilt ) ) + &1894 rot_eigen_rad(3) * SIN( tilt )1897 ( 1.0_wp - COS( tilt_angle ) ) + & 1898 rot_eigen_rad(3) * SIN( tilt_angle ) 1895 1899 rot_coord_trans(inot,2,2) = rot_eigen_rad(2)**2 * & 1896 ( 1.0_wp - COS( tilt ) ) + COS( tilt)1900 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1897 1901 rot_coord_trans(inot,2,3) = rot_eigen_rad(2) * rot_eigen_rad(3) * & 1898 ( 1.0_wp - COS( tilt ) ) - &1899 rot_eigen_rad(1) * SIN( tilt )1902 ( 1.0_wp - COS( tilt_angle ) ) - & 1903 rot_eigen_rad(1) * SIN( tilt_angle ) 1900 1904 rot_coord_trans(inot,3,1) = rot_eigen_rad(3) * rot_eigen_rad(1) * & 1901 ( 1.0_wp - COS( tilt ) ) - &1902 rot_eigen_rad(2) * SIN( tilt )1905 ( 1.0_wp - COS( tilt_angle ) ) - & 1906 rot_eigen_rad(2) * SIN( tilt_angle ) 1903 1907 rot_coord_trans(inot,3,2) = rot_eigen_rad(3) * rot_eigen_rad(2) * & 1904 ( 1.0_wp - COS( tilt ) ) + &1905 rot_eigen_rad(1) * SIN( tilt )1908 ( 1.0_wp - COS( tilt_angle ) ) + & 1909 rot_eigen_rad(1) * SIN( tilt_angle ) 1906 1910 rot_coord_trans(inot,3,3) = rot_eigen_rad(3)**2 * & 1907 ( 1.0_wp - COS( tilt ) ) + COS( tilt)1911 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1908 1912 1909 1913 ! … … 1973 1977 ! 1974 1978 !-- Loop over number of turbines: 1975 DO inot = 1, n turbines1976 1977 cos_yaw = COS( phi_yaw(inot))1978 sin_yaw = SIN( phi_yaw(inot))1979 DO inot = 1, n_turbines 1980 1981 cos_yaw = COS(yaw_angle(inot)) 1982 sin_yaw = SIN(yaw_angle(inot)) 1979 1983 ! 1980 1984 !-- Loop over rings of each turbine: … … 2016 2020 ! 2017 2021 !-- Coordinates of the single segments (center points): 2018 rbx(ring,rseg) = rcx(inot) + cur_r * rote(1)2019 rby(ring,rseg) = rcy(inot) + cur_r * rote(2)2020 rbz(ring,rseg) = rcz(inot) + cur_r * rote(3)2022 rbx(ring,rseg) = hub_x(inot) + cur_r * rote(1) 2023 rby(ring,rseg) = hub_y(inot) + cur_r * rote(2) 2024 rbz(ring,rseg) = hub_z(inot) + cur_r * rote(3) 2021 2025 2022 2026 !-- !-----------------------------------------------------------! … … 2178 2182 !-- Exchange between PEs (information required on each PE): 2179 2183 #if defined( __parallel ) 2180 CALL MPI_ALLREDUCE( u_int_1_l, u_int, n turbines * MAXVAL(nrings) * &2184 CALL MPI_ALLREDUCE( u_int_1_l, u_int, n_turbines * MAXVAL(nrings) * & 2181 2185 MAXVAL(nsegs), MPI_REAL, MPI_SUM, comm2d, ierr ) 2182 CALL MPI_ALLREDUCE( v_int_1_l, v_int, n turbines * MAXVAL(nrings) * &2186 CALL MPI_ALLREDUCE( v_int_1_l, v_int, n_turbines * MAXVAL(nrings) * & 2183 2187 MAXVAL(nsegs), MPI_REAL, MPI_SUM, comm2d, ierr ) 2184 CALL MPI_ALLREDUCE( w_int_1_l, w_int, n turbines * MAXVAL(nrings) * &2188 CALL MPI_ALLREDUCE( w_int_1_l, w_int, n_turbines * MAXVAL(nrings) * & 2185 2189 MAXVAL(nsegs), MPI_REAL, MPI_SUM, comm2d, ierr ) 2186 2190 #else … … 2194 2198 !-- Loop over number of turbines: 2195 2199 2196 DO inot = 1, n turbines2200 DO inot = 1, n_turbines 2197 2201 pit_loop: DO 2198 2202 … … 2200 2204 torque_total(inot) = 0.0_wp 2201 2205 thrust_rotor(inot) = 0.0_wp 2202 pitch_a dd(inot) = pitch_add(inot) + 0.25_wp2203 ! IF ( myid == 0 ) PRINT*, 'Pitch', inot, pitch_a dd(inot)2206 pitch_angle(inot) = pitch_angle(inot) + 0.25_wp 2207 ! IF ( myid == 0 ) PRINT*, 'Pitch', inot, pitch_angle(inot) 2204 2208 ELSE 2205 cos_yaw = COS( phi_yaw(inot))2206 sin_yaw = SIN( phi_yaw(inot))2209 cos_yaw = COS(yaw_angle(inot)) 2210 sin_yaw = SIN(yaw_angle(inot)) 2207 2211 IF ( pitch_control ) THEN 2208 pitch_a dd(inot) = MAX(pitch_add_old(inot) - pitch_rate * &2212 pitch_angle(inot) = MAX(pitch_angle_old(inot) - pitch_rate * & 2209 2213 dt_3d , 0.0_wp ) 2210 2214 ENDIF … … 2261 2265 ! 2262 2266 !-- Coordinates of the single segments (center points): 2263 rbx(ring,rseg) = rcx(inot) + cur_r * rote(1)2264 2265 rby(ring,rseg) = rcy(inot) + cur_r * rote(2)2266 2267 rbz(ring,rseg) = rcz(inot) + cur_r * rote(3)2267 rbx(ring,rseg) = hub_x(inot) + cur_r * rote(1) 2268 2269 rby(ring,rseg) = hub_y(inot) + cur_r * rote(2) 2270 2271 rbz(ring,rseg) = hub_z(inot) + cur_r * rote(3) 2268 2272 2269 2273 ! … … 2296 2300 2297 2301 phi_rel(rseg) = ATAN2( u_rot , & 2298 ( omega_rot(inot) * cur_r -&2302 ( rotor_speed(inot) * cur_r - & 2299 2303 vtheta(rseg) ) ) 2300 2304 … … 2337 2341 ! 2338 2342 !-- Correct with collective pitch angle: 2339 alpha_attack(rseg) = alpha_attack(rseg) - pitch_a dd(inot)2343 alpha_attack(rseg) = alpha_attack(rseg) - pitch_angle(inot) 2340 2344 2341 2345 ! … … 2343 2347 !-- to the rotor: 2344 2348 vrel(rseg) = SQRT( u_rot**2 + & 2345 ( omega_rot(inot) * cur_r -&2349 ( rotor_speed(inot) * cur_r - & 2346 2350 vtheta(rseg) )**2 ) 2347 2351 … … 2393 2397 ( (2.0_wp*pi) / 360.0_wp ) 2394 2398 2395 IF ( t l_cor) THEN2399 IF ( tip_loss_correction ) THEN 2396 2400 2397 2401 !-- Tip loss correction following Schito … … 2401 2405 !-- 2402 2406 tl_factor = ( 2.0 / pi ) * & 2403 ACOS( EXP( -1.0 * ( 3.0 * ( r r(inot) - cur_r ) / &2407 ACOS( EXP( -1.0 * ( 3.0 * ( rotor_radius(inot) - cur_r ) / & 2404 2408 ( 2.0 * cur_r * abs( sin( phi_rel(rseg) ) ) ) ) ) ) 2405 2409 … … 2476 2480 !-- the force balance of the accelerating torque 2477 2481 !-- and the torque of the rotating rotor and generator 2478 om_rate = ( torque_total(inot) * air_dens * gear_eff -&2479 gear_ratio * torque_gen_old(inot) ) / &2480 ( inertia_rot +&2481 gear_ratio * gear_ratio * inertia_gen) * dt_3d2482 om_rate = ( torque_total(inot) * air_density * gear_efficiency - & 2483 gear_ratio * torque_gen_old(inot) ) / & 2484 ( rotor_inertia + & 2485 gear_ratio * gear_ratio * generator_inertia ) * dt_3d 2482 2486 2483 2487 ! 2484 2488 !-- The generator speed is given by the product of gear gear_ratio 2485 2489 !-- and rotor speed 2486 omega_gen(inot) = gear_ratio * ( omega_rot(inot) + om_rate )2490 generator_speed(inot) = gear_ratio * ( rotor_speed(inot) + om_rate ) 2487 2491 2488 2492 ENDIF … … 2495 2499 !-- maximum pitch rate, then the pitch loop is repeated with a pitch 2496 2500 !-- gain 2497 IF ( ( omega_gen(inot) > rated_genspeed ) .AND.&2498 ( pitch_a dd(inot) < 25.0_wp ) .AND. &2499 ( pitch_a dd(inot) < pitch_add_old(inot) + &2501 IF ( ( generator_speed(inot) > generator_speed_rated ) .AND. & 2502 ( pitch_angle(inot) < 25.0_wp ) .AND. & 2503 ( pitch_angle(inot) < pitch_angle_old(inot) + & 2500 2504 pitch_rate * dt_3d ) ) THEN 2501 2505 pitch_sw = .TRUE. … … 2507 2511 ! 2508 2512 !-- The current pitch is saved for the next time step 2509 pitch_a dd_old(inot) = pitch_add(inot)2513 pitch_angle_old(inot) = pitch_angle(inot) 2510 2514 pitch_sw = .FALSE. 2511 2515 ENDIF … … 2554 2558 !-- 2555 2559 !-- Calculation of the boundaries: 2556 i_smear(inot) = CEILING( ( r r(inot) * ABS( roty(inot,1) ) + &2560 i_smear(inot) = CEILING( ( rotor_radius(inot) * ABS( roty(inot,1) ) + & 2557 2561 eps_min ) / dx ) 2558 j_smear(inot) = CEILING( ( r r(inot) * ABS( roty(inot,2) ) + &2562 j_smear(inot) = CEILING( ( rotor_radius(inot) * ABS( roty(inot,2) ) + & 2559 2563 eps_min ) / dy ) 2560 2564 … … 2654 2658 IF ( start_up ) THEN 2655 2659 WDLON = MAX( 1 , NINT( 30.0_wp / dt_3d ) ) ! 30s running mean array 2656 ALLOCATE( wd30(1:n turbines,1:WDLON) )2660 ALLOCATE( wd30(1:n_turbines,1:WDLON) ) 2657 2661 wd30 = 999.0_wp ! Set to dummy value 2658 2662 ALLOCATE( wd30_l(1:WDLON) ) 2659 2663 2660 2664 WDSHO = MAX( 1 , NINT( 2.0_wp / dt_3d ) ) ! 2s running mean array 2661 ALLOCATE( wd2(1:n turbines,1:WDSHO) )2665 ALLOCATE( wd2(1:n_turbines,1:WDSHO) ) 2662 2666 wd2 = 999.0_wp ! Set to dummy value 2663 2667 ALLOCATE( wd2_l(1:WDSHO) ) … … 2669 2673 !-- 2670 2674 !-- Loop over number of turbines: 2671 DO inot = 1, n turbines2675 DO inot = 1, n_turbines 2672 2676 ! 2673 2677 !-- Find processor at i_hub, j_hub … … 2687 2691 CALL wtm_yawcontrol( inot ) 2688 2692 2689 phi_yaw_l(inot) = phi_yaw(inot)2693 yaw_angle_l(inot) = yaw_angle(inot) 2690 2694 2691 2695 ENDIF … … 2697 2701 !-- Transfer of information to the other cpus 2698 2702 #if defined( __parallel ) 2699 CALL MPI_ALLREDUCE( u_inflow_l, u_inflow, n turbines, MPI_REAL,&2703 CALL MPI_ALLREDUCE( u_inflow_l, u_inflow, n_turbines, MPI_REAL, & 2700 2704 MPI_SUM, comm2d, ierr ) 2701 CALL MPI_ALLREDUCE( wdir_l, wdir, n turbines, MPI_REAL, MPI_SUM,&2705 CALL MPI_ALLREDUCE( wdir_l, wdir, n_turbines, MPI_REAL, MPI_SUM, & 2702 2706 comm2d, ierr ) 2703 CALL MPI_ALLREDUCE( phi_yaw_l, phi_yaw, nturbines, MPI_REAL,&2707 CALL MPI_ALLREDUCE( yaw_angle_l, yaw_angle, n_turbines, MPI_REAL, & 2704 2708 MPI_SUM, comm2d, ierr ) 2705 2709 #else 2706 2710 u_inflow = u_inflow_l 2707 2711 wdir = wdir_l 2708 phi_yaw = phi_yaw_l2712 yaw_angle = yaw_angle_l 2709 2713 2710 2714 2711 2715 #endif 2712 DO inot = 1, n turbines2716 DO inot = 1, n_turbines 2713 2717 ! 2714 2718 !-- Update rotor orientation … … 2722 2726 ! 2723 2727 !-- Transfer of information to the other cpus 2724 ! CALL MPI_ALLREDUCE( omega_gen, omega_gen_old, nturbines, &2728 ! CALL MPI_ALLREDUCE( generator_speed, generator_speed_old, n_turbines, & 2725 2729 ! MPI_REAL,MPI_SUM, comm2d, ierr ) 2726 2730 #if defined( __parallel ) 2727 CALL MPI_ALLREDUCE( torque_gen, torque_gen_old, n turbines,&2731 CALL MPI_ALLREDUCE( torque_gen, torque_gen_old, n_turbines, & 2728 2732 MPI_REAL, MPI_SUM, comm2d, ierr ) 2729 CALL MPI_ALLREDUCE( omega_rot_l, omega_rot, nturbines,&2733 CALL MPI_ALLREDUCE( rotor_speed_l, rotor_speed, n_turbines, & 2730 2734 MPI_REAL, MPI_SUM, comm2d, ierr ) 2731 CALL MPI_ALLREDUCE( omega_gen_f, omega_gen_f_old, nturbines,&2735 CALL MPI_ALLREDUCE( generator_speed_f, generator_speed_f_old, n_turbines, & 2732 2736 MPI_REAL, MPI_SUM, comm2d, ierr ) 2733 2737 #else 2734 2738 torque_gen_old = torque_gen 2735 omega_rot = omega_rot_l2736 omega_gen_f_old = omega_gen_f2739 rotor_speed = rotor_speed_l 2740 generator_speed_f_old = generator_speed_f 2737 2741 #endif 2738 2742 … … 2742 2746 2743 2747 2744 DO inot = 1, n turbines2748 DO inot = 1, n_turbines 2745 2749 2746 2750 … … 2748 2752 IF ( myid == 0 ) THEN 2749 2753 IF ( openfile_turb_mod(400+inot)%opened ) THEN 2750 WRITE ( 400+inot, 106 ) time_since_reference_point, omega_rot(inot),&2751 omega_gen(inot), torque_gen_old(inot),&2752 torque_total(inot), pitch_a dd(inot),&2753 torque_gen_old(inot)* omega_gen(inot)*gen_eff,&2754 torque_total(inot)* omega_rot(inot)*air_dens,&2755 thrust_rotor(inot), &2756 wdir(inot)*180.0_wp/pi, &2757 ( phi_yaw(inot))*180.0_wp/pi2754 WRITE ( 400+inot, 106 ) time_since_reference_point, rotor_speed(inot), & 2755 generator_speed(inot), torque_gen_old(inot), & 2756 torque_total(inot), pitch_angle(inot), & 2757 torque_gen_old(inot)*generator_speed(inot)*generator_efficiency, & 2758 torque_total(inot)*rotor_speed(inot)*air_density, & 2759 thrust_rotor(inot), & 2760 wdir(inot)*180.0_wp/pi, & 2761 (yaw_angle(inot))*180.0_wp/pi 2758 2762 2759 2763 ELSE … … 2764 2768 turbine_id ), FORM='FORMATTED' ) 2765 2769 WRITE ( 400+inot, 105 ) inot 2766 WRITE ( 400+inot, 106 ) time_since_reference_point, omega_rot(inot),&2767 omega_gen(inot), torque_gen_old(inot),&2768 torque_total(inot), pitch_a dd(inot),&2769 torque_gen_old(inot)* omega_gen(inot)*gen_eff,&2770 torque_total(inot)* omega_rot(inot)*air_dens,&2771 thrust_rotor(inot), &2772 wdir(inot)*180.0_wp/pi, &2773 ( phi_yaw(inot))*180.0_wp/pi2770 WRITE ( 400+inot, 106 ) time_since_reference_point, rotor_speed(inot), & 2771 generator_speed(inot), torque_gen_old(inot), & 2772 torque_total(inot), pitch_angle(inot), & 2773 torque_gen_old(inot)*generator_speed(inot)*generator_efficiency, & 2774 torque_total(inot)*rotor_speed(inot)*air_density, & 2775 thrust_rotor(inot), & 2776 wdir(inot)*180.0_wp/pi, & 2777 (yaw_angle(inot))*180.0_wp/pi 2774 2778 ENDIF 2775 2779 ENDIF … … 2843 2847 IF ( .NOT. ANY( wd30_l == 999.) ) THEN 2844 2848 2845 missal = SUM( wd30_l ) / SIZE( wd30_l ) - phi_yaw(inot)2846 ! 2847 !-- Check if turbine is missaligned by more than max_miss2848 IF ( ABS( missal ) > max_miss) THEN2849 missal = SUM( wd30_l ) / SIZE( wd30_l ) - yaw_angle(inot) 2850 ! 2851 !-- Check if turbine is missaligned by more than yaw_misalignment_max 2852 IF ( ABS( missal ) > yaw_misalignment_max ) THEN 2849 2853 ! 2850 2854 !-- Check in which direction to yaw … … 2852 2856 ! 2853 2857 !-- Start yawing of turbine 2854 phi_yaw(inot) = phi_yaw(inot) + yawdir(inot) * yaw_speed * dt_3d2858 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2855 2859 doyaw(inot) = .TRUE. 2856 2860 wd30_l = 999. ! fill with dummies again … … 2863 2867 !-- If turbine is already yawing: 2864 2868 !-- Initialize 2 s running mean and yaw until the missalignment is smaller 2865 !-- than min_miss2869 !-- than yaw_misalignment_min 2866 2870 2867 2871 ELSE … … 2876 2880 ! 2877 2881 !-- Calculate missalignment of turbine 2878 missal = SUM( wd2_l - phi_yaw(inot) ) / SIZE( wd2_l )2882 missal = SUM( wd2_l - yaw_angle(inot) ) / SIZE( wd2_l ) 2879 2883 ! 2880 2884 !-- Check if missalignment is still larger than 0.5 degree and if the 2881 2885 !-- yaw direction is still right 2882 IF ( ( ABS( missal ) > min_miss ) .AND.&2886 IF ( ( ABS( missal ) > yaw_misalignment_min ) .AND. & 2883 2887 ( yawdir(inot) == SIGN( 1.0_wp, missal ) ) ) THEN 2884 2888 ! 2885 2889 !-- Continue yawing 2886 phi_yaw(inot) = phi_yaw(inot) + yawdir(inot) * yaw_speed * dt_3d2890 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2887 2891 ELSE 2888 2892 ! … … 2894 2898 ! 2895 2899 !-- Continue yawing 2896 phi_yaw(inot) = phi_yaw(inot) + yawdir(inot) * yaw_speed * dt_3d2900 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2897 2901 ENDIF 2898 2902 … … 2919 2923 ! 2920 2924 !-- Calculate slope constant for region 15 2921 slope15 = ( slope2 * min_reg2 * min_reg2 ) / ( min_reg2 - min_reg15 ) 2925 slope15 = ( region_2_slope * region_2_min * region_2_min ) / & 2926 ( region_2_min - region_15_min ) 2922 2927 ! 2923 2928 !-- Calculate upper limit of slipage region 2924 vs_sysp = rated_genspeed / 1.1_wp2929 vs_sysp = generator_speed_rated / 1.1_wp 2925 2930 ! 2926 2931 !-- Calculate slope of slipage region 2927 slope25 = ( rated_power / rated_genspeed ) /&2928 ( rated_genspeed - vs_sysp )2932 region_25_slope = ( generator_power_rated / generator_speed_rated ) / & 2933 ( generator_speed_rated - vs_sysp ) 2929 2934 ! 2930 2935 !-- Calculate lower limit of slipage region 2931 min_reg25 = ( slope25 - SQRT( slope25 * ( slope25 - 4.0_wp *&2932 slope2 * vs_sysp ) ) ) /&2933 ( 2.0_wp * slope2)2936 region_25_min = ( region_25_slope - SQRT( region_25_slope * & 2937 ( region_25_slope - 4.0_wp * region_2_slope * vs_sysp ) ) ) / & 2938 ( 2.0_wp * region_2_slope ) 2934 2939 ! 2935 2940 !-- Frequency for the simple low pass filter 2936 Fcorner = 0.25_wp2941 Fcorner = 0.25_wp 2937 2942 ! 2938 2943 !-- At the first timestep the torque is set to its maximum to prevent 2939 2944 !-- an overspeeding of the rotor 2940 2945 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 2941 torque_gen_old(:) = max_torque_gen2946 torque_gen_old(:) = generator_torque_max 2942 2947 ENDIF 2943 2948 … … 2968 2973 !-- for the control of the generator torque 2969 2974 lp_coeff = EXP( -2.0_wp * 3.14_wp * dt_3d * Fcorner ) 2970 omega_gen_f(inot) = ( 1.0_wp - lp_coeff ) * omega_gen(inot) + lp_coeff *&2971 omega_gen_f_old(inot)2972 2973 IF ( omega_gen_f(inot) <= min_reg15) THEN2975 generator_speed_f(inot) = ( 1.0_wp - lp_coeff ) * generator_speed(inot) + lp_coeff *& 2976 generator_speed_f_old(inot) 2977 2978 IF ( generator_speed_f(inot) <= region_15_min ) THEN 2974 2979 ! 2975 2980 !-- Region 1: Generator torque is set to zero to accelerate the rotor: 2976 2981 torque_gen(inot) = 0 2977 2982 2978 ELSEIF ( omega_gen_f(inot) <= min_reg2) THEN2983 ELSEIF ( generator_speed_f(inot) <= region_2_min ) THEN 2979 2984 ! 2980 2985 !-- Region 1.5: Generator torque is increasing linearly with rotor speed: 2981 torque_gen(inot) = slope15 * ( omega_gen_f(inot) - min_reg15)2986 torque_gen(inot) = slope15 * ( generator_speed_f(inot) - region_15_min ) 2982 2987 2983 ELSEIF ( omega_gen_f(inot) <= min_reg25) THEN2988 ELSEIF ( generator_speed_f(inot) <= region_25_min ) THEN 2984 2989 ! 2985 2990 !-- Region 2: Generator torque is increased by the square of the generator 2986 2991 !-- speed to keep the TSR optimal: 2987 torque_gen(inot) = slope2 * omega_gen_f(inot) * omega_gen_f(inot)2992 torque_gen(inot) = region_2_slope * generator_speed_f(inot) * generator_speed_f(inot) 2988 2993 2989 ELSEIF ( omega_gen_f(inot) < rated_genspeed ) THEN2994 ELSEIF ( generator_speed_f(inot) < generator_speed_rated ) THEN 2990 2995 ! 2991 2996 !-- Region 2.5: Slipage region between 2 and 3: 2992 torque_gen(inot) = slope25 * ( omega_gen_f(inot) - vs_sysp )2997 torque_gen(inot) = region_25_slope * ( generator_speed_f(inot) - vs_sysp ) 2993 2998 2994 2999 ELSE … … 2996 3001 !-- Region 3: Generator torque is antiproportional to the rotor speed to 2997 3002 !-- keep the power constant: 2998 torque_gen(inot) = rated_power / omega_gen_f(inot)3003 torque_gen(inot) = generator_power_rated / generator_speed_f(inot) 2999 3004 3000 3005 ENDIF … … 3002 3007 !-- Calculate torque rate and confine with a max 3003 3008 trq_rate = ( torque_gen(inot) - torque_gen_old(inot) ) / dt_3d 3004 trq_rate = MIN( MAX( trq_rate, -1.0_wp * max_trq_rate ), max_trq_rate ) 3009 trq_rate = MIN( MAX( trq_rate, -1.0_wp * generator_torque_rate_max ), & 3010 generator_torque_rate_max ) 3005 3011 ! 3006 3012 !-- Calculate new gen torque and confine with max torque 3007 3013 torque_gen(inot) = torque_gen_old(inot) + trq_rate * dt_3d 3008 torque_gen(inot) = MIN( torque_gen(inot), max_torque_gen)3014 torque_gen(inot) = MIN( torque_gen(inot), generator_torque_max ) 3009 3015 ! 3010 3016 !-- Overwrite values for next timestep 3011 omega_rot_l(inot) = omega_gen(inot) / gear_ratio3017 rotor_speed_l(inot) = generator_speed(inot) / gear_ratio 3012 3018 3013 3019 … … 3194 3200 !-- At the first call of this routine write the spatial coordinates. 3195 3201 IF ( .NOT. initial_write_coordinates ) THEN 3196 ALLOCATE ( output_values_1d_target(1:n turbines) )3197 output_values_1d_target = rcx(1:nturbines)3202 ALLOCATE ( output_values_1d_target(1:n_turbines) ) 3203 output_values_1d_target = hub_x(1:n_turbines) 3198 3204 output_values_1d_pointer => output_values_1d_target 3199 3205 return_value = dom_write_var( nc_filename, & … … 3201 3207 values_realwp_1d = output_values_1d_pointer, & 3202 3208 bounds_start = (/1/), & 3203 bounds_end = (/n turbines/) )3204 3205 output_values_1d_target = rcy(1:nturbines)3209 bounds_end = (/n_turbines/) ) 3210 3211 output_values_1d_target = hub_y(1:n_turbines) 3206 3212 output_values_1d_pointer => output_values_1d_target 3207 3213 return_value = dom_write_var( nc_filename, & … … 3209 3215 values_realwp_1d = output_values_1d_pointer, & 3210 3216 bounds_start = (/1/), & 3211 bounds_end = (/n turbines/) )3212 3213 output_values_1d_target = rcz(1:nturbines)3217 bounds_end = (/n_turbines/) ) 3218 3219 output_values_1d_target = hub_z(1:n_turbines) 3214 3220 output_values_1d_pointer => output_values_1d_target 3215 3221 return_value = dom_write_var( nc_filename, & … … 3217 3223 values_realwp_1d = output_values_1d_pointer, & 3218 3224 bounds_start = (/1/), & 3219 bounds_end = (/n turbines/) )3225 bounds_end = (/n_turbines/) ) 3220 3226 3221 3227 initial_write_coordinates = .TRUE. … … 3225 3231 t_ind = t_ind + 1 3226 3232 3227 ALLOCATE ( output_values_1d_target(1:n turbines) )3228 output_values_1d_target = omega_rot(:)3233 ALLOCATE ( output_values_1d_target(1:n_turbines) ) 3234 output_values_1d_target = rotor_speed(:) 3229 3235 output_values_1d_pointer => output_values_1d_target 3230 3236 … … 3233 3239 values_realwp_1d = output_values_1d_pointer, & 3234 3240 bounds_start = (/1, t_ind/), & 3235 bounds_end = (/n turbines, t_ind /) )3236 3237 output_values_1d_target = omega_gen(:)3241 bounds_end = (/n_turbines, t_ind /) ) 3242 3243 output_values_1d_target = generator_speed(:) 3238 3244 output_values_1d_pointer => output_values_1d_target 3239 3245 return_value = dom_write_var( nc_filename, & … … 3241 3247 values_realwp_1d = output_values_1d_pointer, & 3242 3248 bounds_start = (/1, t_ind/), & 3243 bounds_end = (/n turbines, t_ind /) )3249 bounds_end = (/n_turbines, t_ind /) ) 3244 3250 3245 3251 output_values_1d_target = torque_gen_old(:) … … 3250 3256 values_realwp_1d = output_values_1d_pointer, & 3251 3257 bounds_start = (/1, t_ind/), & 3252 bounds_end = (/n turbines, t_ind /) )3258 bounds_end = (/n_turbines, t_ind /) ) 3253 3259 3254 3260 output_values_1d_target = torque_total(:) … … 3259 3265 values_realwp_1d = output_values_1d_pointer, & 3260 3266 bounds_start = (/1, t_ind/), & 3261 bounds_end = (/n turbines, t_ind /) )3262 3263 output_values_1d_target = pitch_a dd(:)3267 bounds_end = (/n_turbines, t_ind /) ) 3268 3269 output_values_1d_target = pitch_angle(:) 3264 3270 output_values_1d_pointer => output_values_1d_target 3265 3271 … … 3268 3274 values_realwp_1d = output_values_1d_pointer, & 3269 3275 bounds_start = (/1, t_ind/), & 3270 bounds_end = (/n turbines, t_ind /) )3271 3272 output_values_1d_target = torque_gen_old(:)* omega_gen(:)*gen_eff3276 bounds_end = (/n_turbines, t_ind /) ) 3277 3278 output_values_1d_target = torque_gen_old(:)*generator_speed(:)*generator_efficiency 3273 3279 output_values_1d_pointer => output_values_1d_target 3274 3280 … … 3277 3283 values_realwp_1d = output_values_1d_pointer, & 3278 3284 bounds_start = (/1, t_ind/), & 3279 bounds_end = (/n turbines, t_ind /) )3280 3281 DO inot = 1, n turbines3282 output_values_1d_target(inot) = torque_total(inot)* omega_rot(inot)*air_dens3285 bounds_end = (/n_turbines, t_ind /) ) 3286 3287 DO inot = 1, n_turbines 3288 output_values_1d_target(inot) = torque_total(inot)*rotor_speed(inot)*air_density 3283 3289 ENDDO 3284 3290 output_values_1d_pointer => output_values_1d_target … … 3288 3294 values_realwp_1d = output_values_1d_pointer, & 3289 3295 bounds_start = (/1, t_ind/), & 3290 bounds_end = (/n turbines, t_ind /) )3296 bounds_end = (/n_turbines, t_ind /) ) 3291 3297 3292 3298 output_values_1d_target = thrust_rotor(:) … … 3297 3303 values_realwp_1d = output_values_1d_pointer, & 3298 3304 bounds_start = (/1, t_ind/), & 3299 bounds_end = (/n turbines, t_ind /) )3305 bounds_end = (/n_turbines, t_ind /) ) 3300 3306 3301 3307 output_values_1d_target = wdir(:)*180.0_wp/pi … … 3306 3312 values_realwp_1d = output_values_1d_pointer, & 3307 3313 bounds_start = (/1, t_ind/), & 3308 bounds_end = (/n turbines, t_ind /) )3309 3310 output_values_1d_target = ( phi_yaw(:))*180.0_wp/pi3314 bounds_end = (/n_turbines, t_ind /) ) 3315 3316 output_values_1d_target = (yaw_angle(:))*180.0_wp/pi 3311 3317 output_values_1d_pointer => output_values_1d_target 3312 3318 … … 3315 3321 values_realwp_1d = output_values_1d_pointer, & 3316 3322 bounds_start = (/1, t_ind/), & 3317 bounds_end = (/n turbines, t_ind /) )3323 bounds_end = (/n_turbines, t_ind /) ) 3318 3324 3319 3325 output_values_0d_target = time_since_reference_point … … 3323 3329 'time', & 3324 3330 values_realwp_0d = output_values_0d_pointer, & 3325 bounds_start = (/t_ind/), 3331 bounds_start = (/t_ind/), & 3326 3332 bounds_end = (/t_ind/) ) 3327 3333 -
palm/trunk/TESTS/cases/wind_turbine_model/INPUT/wind_turbine_model_p3d
r4069 r4447 38 38 time_turbine_on = 1.0, 39 39 40 n turbines = 4,41 rcx = 300.0, 600.0, 300.0, 600.0,42 rcy = 300.0, 300.0, 600.0, 600.0,43 rcz = 480.0, 480.0, 480.0, 480.0,44 r r= 63.0, 63.0, 63.0, 63.0,45 rnac= 1.5, 1.5, 1.5, 1.5,46 dtow= 4.0, 4.0, 4.0, 4.0,40 n_turbines = 4, 41 hub_x = 300.0, 600.0, 300.0, 600.0, 42 hub_y = 300.0, 300.0, 600.0, 600.0, 43 hub_z = 480.0, 480.0, 480.0, 480.0, 44 rotor_radius = 63.0, 63.0, 63.0, 63.0, 45 nacelle_radius = 1.5, 1.5, 1.5, 1.5, 46 tower_diameter = 4.0, 4.0, 4.0, 4.0, 47 47 48 omega_rot= 0.9, 0.9, 0.9, 0.9,49 phi_yaw= -20.0, 10.0, 0.0, 30.0,50 tilt = 0.0,51 rated_power= 5296610.0,48 rotor_speed = 0.9, 0.9, 0.9, 0.9, 49 yaw_angle = -20.0, 10.0, 0.0, 30.0, 50 tilt_angle = 0.0, 51 generator_power_rated = 5296610.0, 52 52 speed_control = .T., 53 53 pitch_control = .T., 54 pitch_a dd= 0.0,54 pitch_angle = 0.0, 55 55 yaw_control = .T., 56 56 yaw_speed = 0.01, 57 57 58 t urb_cd_tower= 1.2, 1.2, 1.2, 1.2,58 tower_cd = 1.2, 1.2, 1.2, 1.2, 59 59 60 segment_length = 1.0,61 segment_width = 0.5,60 segment_length_tangential = 1.0, 61 segment_width_radial = 0.5, 62 62 63 n airfoils = 8, /63 n_airfoils = 8, /
Note: See TracChangeset
for help on using the changeset viewer.
