Changeset 2836 for palm/trunk/SOURCE/wind_turbine_model_mod.f90
 Timestamp:
 Feb 26, 2018 1:40:05 PM (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/wind_turbine_model_mod.f90
r2792 r2836 26 26 !  27 27 ! $Id$ 28 ! Bugfix in interpolation of lift and drag coefficients on fine grid of radius 29 ! segments and angles of attack, speedup of the initialization of the wind 30 ! turbine model 31 ! 32 ! 2792 20180207 06:45:29Z Giersch 28 33 ! omega_rot_l has to be calculated after determining the indices for the hub in 29 34 ! case of restart runs … … 1406 1411 1407 1412 ! 1408 ! For each possible radial position (resolution: 0.1 m > 63 0 values) and1409 ! each possible angle of attack (resolution: 0.01 degrees > 36000values!)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!) 1410 1415 ! determine the lift and drag coefficient by interpolating between the 1411 1416 ! tabulated values of each table (interpolate to current angle of attack) 1412 1417 ! and between the tables (interpolate to current radial position): 1413 1418 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) ) 1418 1423 1419 1424 radres = INT( rr(1) * 10.0_wp ) + 1_iwp 1420 1425 dlenbl_int = INT( 360.0_wp / accu_cl_cd_tab ) + 1_iwp 1421 1426 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) ) 1426 1429 1427 1430 DO iir = 1, radres ! loop over radius 1428 1431 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 1432 1480 alpha_attack_i = 180.0_wp + REAL( iialpha1 ) * accu_cl_cd_tab 1433 1481 ialpha = 1 1434 1482 1435 1483 DO WHILE ( ( alpha_attack_i > alpha_attack_tab(ialpha) ) .AND. (ialpha <= dlen ) ) 1436 1484 ialpha = ialpha + 1 1437 1485 ENDDO 1438 !1439 ! Find position in table1440 lct = MINLOC( ABS( trad1  cur_r ) )1441 ! lct(1) = lct(1)1442 1443 IF ( ( trad1(lct(1))  cur_r ) .GT. 0.0 ) THEN1444 lct(1) = lct(1)  11445 ENDIF1446 1447 trow = lct(1)1448 !1449 ! Calculate weights for interpolation1450 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 NaN1456 weight_a = 0.5_wp ! then do interpolate in between same twice1457 weight_b = 0.5_wp ! using 0.5 as weight1458 ENDIF1459 1460 IF ( t1 == 0 .AND. t2 == 0 ) THEN1461 turb_cd_sel1 = 0.0_wp1462 turb_cd_sel2 = 0.0_wp1463 turb_cl_sel1 = 0.0_wp1464 turb_cl_sel2 = 0.0_wp1465 ELSE1466 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 ENDIF1471 1486 1472 1487 ! … … 1497 1512 ( weight_a * turb_cd_sel1(ialpha) + & 1498 1513 weight_b * turb_cd_sel2(ialpha) ) 1499 1514 1500 1515 ENDDO ! end loop over angles of attack 1501 1516 1502 1517 ENDDO ! end loop over radius 1503 1518 1519 1504 1520 END SUBROUTINE wtm_read_blade_tables 1505 1521 … … 2016 2032 ! Determine index of current angle of attack, needed for 2017 2033 ! 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): 2020 2036 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 2022 2038 ! 2023 2039 ! Determine index of current radial position, needed for
Note: See TracChangeset
for help on using the changeset viewer.