Ignore:
Timestamp:
Feb 26, 2018 1:40:05 PM (4 years ago)
Author:
Giersch
Message:

Bugfix in interpolation of lift and drag coefficients during initialization, default value of dt_dopr_listing is set to zero to allow header output

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/wind_turbine_model_mod.f90

    r2792 r2836  
    2626! -----------------
    2727! $Id$
     28! Bugfix in interpolation of lift and drag coefficients on fine grid of radius
     29! segments and angles of attack, speed-up of the initialization of the wind
     30! turbine model
     31!
     32! 2792 2018-02-07 06:45:29Z Giersch
    2833! omega_rot_l has to be calculated after determining the indices for the hub in
    2934! case of restart runs
     
    14061411
    14071412!
    1408 !--    For each possible radial position (resolution: 0.1 m --> 630 values) and
    1409 !--    each possible angle of attack (resolution: 0.01 degrees --> 36000 values!)
     1413!--    For each possible radial position (resolution: 0.1 m --> 631 values if rr(1)=63m)
     1414!--    and each possible angle of attack (resolution: 0.1 degrees --> 3601 values!)
    14101415!--    determine the lift and drag coefficient by interpolating between the
    14111416!--    tabulated values of each table (interpolate to current angle of attack)
    14121417!--    and between the tables (interpolate to current radial position):
    14131418
    1414        ALLOCATE( turb_cl_sel1(0:dlenbl) ) 
    1415        ALLOCATE( turb_cl_sel2(0:dlenbl) ) 
    1416        ALLOCATE( turb_cd_sel1(0:dlenbl) )
    1417        ALLOCATE( turb_cd_sel2(0:dlenbl) )
     1419       ALLOCATE( turb_cl_sel1(1:dlenbl) ) 
     1420       ALLOCATE( turb_cl_sel2(1:dlenbl) ) 
     1421       ALLOCATE( turb_cd_sel1(1:dlenbl) )
     1422       ALLOCATE( turb_cd_sel2(1:dlenbl) )
    14181423
    14191424       radres     = INT( rr(1) * 10.0_wp ) + 1_iwp
    14201425       dlenbl_int = INT( 360.0_wp / accu_cl_cd_tab ) + 1_iwp
    14211426
    1422 
    1423        ALLOCATE( turb_cl_tab(0:dlenbl_int,1:radres) )
    1424        ALLOCATE( turb_cd_tab(0:dlenbl_int,1:radres) )
    1425 
     1427       ALLOCATE( turb_cl_tab(1:dlenbl_int,1:radres) )
     1428       ALLOCATE( turb_cd_tab(1:dlenbl_int,1:radres) )
    14261429
    14271430       DO iir = 1, radres ! loop over radius
    14281431
    1429           DO iialpha = 1, dlenbl_int  ! loop over angles
    1430 
    1431              cur_r = ( iir - 1_iwp ) * 0.1_wp             
     1432          cur_r = ( iir - 1_iwp ) * 0.1_wp
     1433!
     1434!--       Find position in table 1
     1435          lct = MINLOC( ABS( trad1 - cur_r ) )
     1436!             lct(1) = lct(1)
     1437
     1438          IF ( ( trad1(lct(1)) - cur_r ) .GT. 0.0 )  THEN
     1439             lct(1) = lct(1) - 1
     1440          ENDIF
     1441
     1442          trow = lct(1)
     1443!
     1444!--       Calculate weights for radius interpolation
     1445          weight_a = ( trad2(trow) - cur_r ) / ( trad2(trow) - trad1(trow) )
     1446          weight_b = ( cur_r - trad1(trow) ) / ( trad2(trow) - trad1(trow) )
     1447          t1 = ttoint1(trow)
     1448          t2 = ttoint2(trow)
     1449
     1450          IF ( t1 .EQ. t2 ) THEN  ! if both are the same, the weights are NaN
     1451             weight_a = 0.5_wp    ! then do interpolate in between same twice
     1452             weight_b = 0.5_wp    ! using 0.5 as weight
     1453          ENDIF
     1454
     1455          IF ( t1 == 0 .AND. t2 == 0 ) THEN
     1456             turb_cd_sel1 = 0.0_wp
     1457             turb_cd_sel2 = 0.0_wp
     1458             turb_cl_sel1 = 0.0_wp
     1459             turb_cl_sel2 = 0.0_wp
     1460
     1461             turb_cd_tab(1,iir) = 0.0_wp  ! For -180° (iialpha=1) the values   
     1462             turb_cl_tab(1,iir) = 0.0_wp  ! for each radius has to be set
     1463                                          ! explicitly             
     1464          ELSE
     1465             turb_cd_sel1 = turb_cd_table(:,t1)
     1466             turb_cd_sel2 = turb_cd_table(:,t2)
     1467             turb_cl_sel1 = turb_cl_table(:,t1)
     1468             turb_cl_sel2 = turb_cl_table(:,t2)
     1469!
     1470!--          For -180° (iialpha=1) the values for each radius has to be set
     1471!--          explicitly
     1472             turb_cd_tab(1,iir) = ( weight_a * turb_cd_table(1,t1) + weight_b  &
     1473                                  * turb_cd_table(1,t2) )   
     1474             turb_cl_tab(1,iir) = ( weight_a * turb_cl_table(1,t1) + weight_b  &
     1475                                  * turb_cl_table(1,t2) )
     1476          ENDIF
     1477
     1478          DO iialpha = 2, dlenbl_int  ! loop over angles
     1479             
    14321480             alpha_attack_i = -180.0_wp + REAL( iialpha-1 ) * accu_cl_cd_tab
    14331481             ialpha = 1
    1434 
     1482             
    14351483             DO WHILE ( ( alpha_attack_i > alpha_attack_tab(ialpha) ) .AND. (ialpha <= dlen ) )
    14361484                ialpha = ialpha + 1
    14371485             ENDDO
    1438 !
    1439 !--          Find position in table
    1440              lct = MINLOC( ABS( trad1 - cur_r ) )
    1441 !                lct(1) = lct(1)
    1442 
    1443              IF ( ( trad1(lct(1)) - cur_r ) .GT. 0.0 )  THEN
    1444                 lct(1) = lct(1) - 1
    1445              ENDIF
    1446 
    1447              trow = lct(1)
    1448 !
    1449 !--          Calculate weights for interpolation
    1450              weight_a = ( trad2(trow) - cur_r ) / ( trad2(trow) - trad1(trow) )
    1451              weight_b = ( cur_r - trad1(trow) ) / ( trad2(trow) - trad1(trow) )
    1452              t1 = ttoint1(trow)
    1453              t2 = ttoint2(trow)
    1454 
    1455              IF ( t1 .EQ. t2 ) THEN  ! if both are the same, the weights are NaN
    1456                 weight_a = 0.5_wp    ! then do interpolate in between same twice
    1457                 weight_b = 0.5_wp    ! using 0.5 as weight
    1458              ENDIF
    1459 
    1460              IF ( t1 == 0 .AND. t2 == 0 ) THEN
    1461                 turb_cd_sel1 = 0.0_wp
    1462                 turb_cd_sel2 = 0.0_wp
    1463                 turb_cl_sel1 = 0.0_wp
    1464                 turb_cl_sel2 = 0.0_wp
    1465              ELSE
    1466                 turb_cd_sel1 = turb_cd_table(:,t1)
    1467                 turb_cd_sel2 = turb_cd_table(:,t2)
    1468                 turb_cl_sel1 = turb_cl_table(:,t1)
    1469                 turb_cl_sel2 = turb_cl_table(:,t2)
    1470              ENDIF
    14711486
    14721487!
     
    14971512                                        ( weight_a * turb_cd_sel1(ialpha) +    &
    14981513                                          weight_b * turb_cd_sel2(ialpha) )
    1499    
     1514            
    15001515          ENDDO   ! end loop over angles of attack
    1501 
     1516         
    15021517       ENDDO   ! end loop over radius
    1503    
     1518
     1519
    15041520    END SUBROUTINE wtm_read_blade_tables
    15051521
     
    20162032!--                Determine index of current angle of attack, needed for
    20172033!--                finding the appropriate interpolated values of the lift and
    2018 !--                drag coefficients (-180.0 degrees = 0, +180.0 degrees = 36000,
    2019 !--                so one index every 0.01 degrees):
     2034!--                drag coefficients (-180.0 degrees = 1, +180.0 degrees = 3601,
     2035!--                so one index every 0.1 degrees):
    20202036                   iialpha = CEILING( ( alpha_attack(rseg) + 180.0_wp )        &
    2021                                       * ( 1.0_wp / accu_cl_cd_tab ) )
     2037                                      * ( 1.0_wp / accu_cl_cd_tab ) ) + 1.0_wp
    20222038!
    20232039!--                Determine index of current radial position, needed for
Note: See TracChangeset for help on using the changeset viewer.