Changeset 4465 for palm/trunk/SOURCE/wind_turbine_model_mod.f90
- Timestamp:
- Mar 20, 2020 11:35:48 AM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/wind_turbine_model_mod.f90 (modified) (38 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wind_turbine_model_mod.f90
r4460 r4465 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Removed old ASCII outputm, some syntax layout adjustments, added output for 29 ! rotor and tower diameters. Added warning message in case of NetCDF 3 (no 30 ! WTM output file will be produced). 31 ! 32 ! 4460 2020-03-12 16:47:30Z oliver.maas 28 33 ! allow for simulating up to 10 000 wind turbines 29 34 ! … … 258 263 REAL(wp) :: yaw_misalignment_min = 0.008726_wp !< minimum yaw missalignment for which the actuator stops [rad] 259 264 260 !261 !-- Set flag for output files TURBINE_PARAMETERS262 TYPE file_status263 LOGICAL :: opened, opened_before264 END TYPE file_status265 266 TYPE(file_status), DIMENSION(500) :: openfile_turb_mod = &267 file_status(.FALSE.,.FALSE.)268 265 269 266 ! … … 542 539 !-- Set flag that indicates that the wind turbine model is switched on 543 540 wind_turbine = .TRUE. 544 541 545 542 GOTO 12 546 543 … … 1391 1388 ! 1392 1389 !-- Create NetCDF output file 1393 nc_filename = 'DATA_1D_TS_WTM_NETCDF' // TRIM( coupling_char ) 1390 #if defined( __netcdf4 ) 1391 nc_filename = 'DATA_1D_TS_WTM_NETCDF' // TRIM( coupling_char ) 1394 1392 return_value = dom_def_file( nc_filename, 'netcdf4-serial' ) 1395 1393 #else 1394 message_string = 'Wind turbine model output requires NetCDF version 4. ' //& 1395 'No output file will be created.' 1396 CALL message( 'wtm_init_output', 'PA0672', & 1397 0, 1, 0, 6, 0 ) 1398 #endif 1396 1399 IF ( myid == 0 ) THEN 1397 1400 ! … … 1408 1411 DEALLOCATE( ndim ) 1409 1412 1410 !1411 !-- time1412 1413 return_value = dom_def_dim( nc_filename, & 1413 1414 dimension_name = 'time', & … … 1416 1417 values_realwp = (/0.0_wp/) ) 1417 1418 1418 !1419 !-- x1420 1419 variable_name = 'x' 1421 return_value = dom_def_var( nc_filename, &1420 return_value = dom_def_var( nc_filename, & 1422 1421 variable_name = variable_name, & 1423 1422 dimension_names = (/'turbine'/), & 1424 1423 output_type = 'real32' ) 1425 ! 1426 !-- y 1424 1427 1425 variable_name = 'y' 1428 return_value = dom_def_var( nc_filename, &1426 return_value = dom_def_var( nc_filename, & 1429 1427 variable_name = variable_name, & 1430 1428 dimension_names = (/'turbine'/), & … … 1432 1430 1433 1431 variable_name = 'z' 1434 return_value = dom_def_var( nc_filename, &1432 return_value = dom_def_var( nc_filename, & 1435 1433 variable_name = variable_name, & 1436 1434 dimension_names = (/'turbine'/), & 1437 1435 output_type = 'real32' ) 1438 1439 1440 return_value = dom_def_att( nc_filename, & 1436 1437 1438 variable_name = 'rotor_diameter' 1439 return_value = dom_def_var( nc_filename, & 1440 variable_name = variable_name, & 1441 dimension_names = (/'turbine'/), & 1442 output_type = 'real32' ) 1443 1444 variable_name = 'tower_diameter' 1445 return_value = dom_def_var( nc_filename, & 1446 variable_name = variable_name, & 1447 dimension_names = (/'turbine'/), & 1448 output_type = 'real32' ) 1449 1450 return_value = dom_def_att( nc_filename, & 1441 1451 variable_name = 'time', & 1442 attribute_name = 'units', &1452 attribute_name = 'units', & 1443 1453 value = 'seconds since ' // origin_date_time ) 1444 1445 return_value = dom_def_att( nc_filename, &1454 1455 return_value = dom_def_att( nc_filename, & 1446 1456 variable_name = 'x', & 1447 attribute_name = 'units', &1457 attribute_name = 'units', & 1448 1458 value = 'm' ) 1449 1459 1450 return_value = dom_def_att( nc_filename, &1460 return_value = dom_def_att( nc_filename, & 1451 1461 variable_name = 'y', & 1452 attribute_name = 'units', &1462 attribute_name = 'units', & 1453 1463 value = 'm' ) 1454 1464 1455 return_value = dom_def_att( nc_filename, &1465 return_value = dom_def_att( nc_filename, & 1456 1466 variable_name = 'z', & 1457 attribute_name = 'units', & 1458 value = 'm' ) 1459 1460 return_value = dom_def_att( nc_filename, & 1467 attribute_name = 'units', & 1468 value = 'm' ) 1469 1470 return_value = dom_def_att( nc_filename, & 1471 variable_name = 'rotor_diameter', & 1472 attribute_name = 'units', & 1473 value = 'm' ) 1474 1475 return_value = dom_def_att( nc_filename, & 1476 variable_name = 'tower_diameter', & 1477 attribute_name = 'units', & 1478 value = 'm' ) 1479 1480 return_value = dom_def_att( nc_filename, & 1461 1481 variable_name = 'x', & 1462 attribute_name = 'long_name', 1482 attribute_name = 'long_name', & 1463 1483 value = 'x location of rotor center' ) 1464 1484 1465 return_value = dom_def_att( nc_filename, &1485 return_value = dom_def_att( nc_filename, & 1466 1486 variable_name = 'y', & 1467 attribute_name = 'long_name', 1487 attribute_name = 'long_name', & 1468 1488 value = 'y location of rotor center' ) 1469 1489 1470 return_value = dom_def_att( nc_filename, &1490 return_value = dom_def_att( nc_filename, & 1471 1491 variable_name = 'z', & 1472 attribute_name = 'long_name', 1492 attribute_name = 'long_name', & 1473 1493 value = 'z location of rotor center' ) 1474 1475 1476 return_value = dom_def_att( nc_filename, & 1494 1495 return_value = dom_def_att( nc_filename, & 1477 1496 variable_name = 'turbine_name', & 1478 attribute_name = 'long_name', & 1479 value = 'turbine name') 1480 1481 return_value = dom_def_att( nc_filename, & 1497 attribute_name = 'long_name', & 1498 value = 'turbine name') 1499 1500 return_value = dom_def_att( nc_filename, & 1501 variable_name = 'rotor_diameter', & 1502 attribute_name = 'long_name', & 1503 value = 'rotor diameter') 1504 1505 return_value = dom_def_att( nc_filename, & 1506 variable_name = 'tower_diameter', & 1507 attribute_name = 'long_name', & 1508 value = 'tower diameter') 1509 1510 return_value = dom_def_att( nc_filename, & 1482 1511 variable_name = 'time', & 1483 attribute_name = 'standard_name', 1512 attribute_name = 'standard_name', & 1484 1513 value = 'time') 1485 1514 1486 return_value = dom_def_att( nc_filename, &1515 return_value = dom_def_att( nc_filename, & 1487 1516 variable_name = 'time', & 1488 1517 attribute_name = 'axis', & 1489 1518 value = 'T') 1490 1519 1491 return_value = dom_def_att( nc_filename, &1520 return_value = dom_def_att( nc_filename, & 1492 1521 variable_name = 'x', & 1493 1522 attribute_name = 'axis', & 1494 1523 value = 'X' ) 1495 1524 1496 return_value = dom_def_att( nc_filename, &1525 return_value = dom_def_att( nc_filename, & 1497 1526 variable_name = 'y', & 1498 1527 attribute_name = 'axis', & 1499 value = 'Y' ) 1500 1501 1502 variable_name = 'rotor_speed' 1503 return_value = dom_def_var( nc_filename, & 1504 variable_name = variable_name, & 1505 dimension_names = (/ 'turbine', 'time ' /), & 1528 value = 'Y' ) 1529 1530 1531 variable_name = 'generator_power' 1532 return_value = dom_def_var( nc_filename, & 1533 variable_name = variable_name, & 1534 dimension_names = (/ 'turbine', 'time ' /), & 1535 output_type = 'real32' ) 1536 1537 return_value = dom_def_att( nc_filename, & 1538 variable_name = variable_name, & 1539 attribute_name = 'units', & 1540 value = 'W' ) 1541 1542 variable_name = 'generator_speed' 1543 return_value = dom_def_var( nc_filename, & 1544 variable_name = variable_name, & 1545 dimension_names = (/ 'turbine', 'time ' /), & 1546 output_type = 'real32' ) 1547 1548 return_value = dom_def_att( nc_filename, & 1549 variable_name = variable_name, & 1550 attribute_name = 'units', & 1551 value = 'rad/s' ) 1552 1553 variable_name = 'generator_torque' 1554 return_value = dom_def_var( nc_filename, & 1555 variable_name = variable_name, & 1556 dimension_names = (/ 'turbine', 'time ' /), & 1506 1557 output_type = 'real32' ) 1507 1558 1508 return_value = dom_def_att( nc_filename, & 1509 variable_name = variable_name, & 1510 attribute_name = 'units', & 1511 value = 'rad/s' ) 1512 1513 variable_name = 'generator_speed' 1514 return_value = dom_def_var( nc_filename, & 1515 variable_name = variable_name, & 1516 dimension_names = (/ 'turbine', 'time ' /), & 1559 return_value = dom_def_att( nc_filename, & 1560 variable_name = variable_name, & 1561 attribute_name = 'units', & 1562 value = 'Nm' ) 1563 1564 variable_name = 'pitch_angle' 1565 return_value = dom_def_var( nc_filename, & 1566 variable_name = variable_name, & 1567 dimension_names = (/ 'turbine', 'time ' /), & 1568 output_type = 'real32' ) 1569 1570 return_value = dom_def_att( nc_filename, & 1571 variable_name = variable_name, & 1572 attribute_name = 'units', & 1573 value = 'degrees' ) 1574 1575 variable_name = 'rotor_power' 1576 return_value = dom_def_var( nc_filename, & 1577 variable_name = variable_name, & 1578 dimension_names = (/ 'turbine', 'time ' /), & 1517 1579 output_type = 'real32' ) 1518 1580 1519 return_value = dom_def_att( nc_filename, & 1520 variable_name = variable_name, & 1521 attribute_name = 'units', & 1522 value = 'rad/s' ) 1523 1524 1525 variable_name = 'generator_torque' 1526 return_value = dom_def_var( nc_filename, & 1527 variable_name = variable_name, & 1528 dimension_names = (/ 'turbine', 'time ' /), & 1581 return_value = dom_def_att( nc_filename, & 1582 variable_name = variable_name, & 1583 attribute_name = 'units', & 1584 value = 'W' ) 1585 1586 variable_name = 'rotor_speed' 1587 return_value = dom_def_var( nc_filename, & 1588 variable_name = variable_name, & 1589 dimension_names = (/ 'turbine', 'time ' /), & 1529 1590 output_type = 'real32' ) 1530 1591 1531 return_value = dom_def_att( nc_filename, &1532 variable_name = variable_name, &1533 attribute_name = 'units', &1534 value = ' Nm' )1535 1536 variable_name = 'rotor_t orque'1537 return_value = dom_def_var( nc_filename, &1538 variable_name = variable_name, &1539 dimension_names = (/ 'turbine', 'time ' /), 1592 return_value = dom_def_att( nc_filename, & 1593 variable_name = variable_name, & 1594 attribute_name = 'units', & 1595 value = 'rad/s' ) 1596 1597 variable_name = 'rotor_thrust' 1598 return_value = dom_def_var( nc_filename, & 1599 variable_name = variable_name, & 1600 dimension_names = (/ 'turbine', 'time ' /), & 1540 1601 output_type = 'real32' ) 1541 1602 1542 return_value = dom_def_att( nc_filename, &1543 variable_name = variable_name, &1544 attribute_name = 'units', &1545 value = 'N m' )1546 1547 variable_name = ' pitch_angle'1548 return_value = dom_def_var( nc_filename, &1549 variable_name = variable_name, &1550 dimension_names = (/ 'turbine', 'time ' /), 1603 return_value = dom_def_att( nc_filename, & 1604 variable_name = variable_name, & 1605 attribute_name = 'units', & 1606 value = 'N' ) 1607 1608 variable_name = 'rotor_torque' 1609 return_value = dom_def_var( nc_filename, & 1610 variable_name = variable_name, & 1611 dimension_names = (/ 'turbine', 'time ' /), & 1551 1612 output_type = 'real32' ) 1552 1613 1553 return_value = dom_def_att( nc_filename, &1554 variable_name = variable_name, &1555 attribute_name = 'units', &1556 value = ' degrees' )1557 1558 variable_name = ' generator_power'1559 return_value = dom_def_var( nc_filename, &1560 variable_name = variable_name, &1561 dimension_names = (/ 'turbine', 'time ' /), 1614 return_value = dom_def_att( nc_filename, & 1615 variable_name = variable_name, & 1616 attribute_name = 'units', & 1617 value = 'Nm' ) 1618 1619 variable_name = 'wind_direction' 1620 return_value = dom_def_var( nc_filename, & 1621 variable_name = variable_name, & 1622 dimension_names = (/ 'turbine', 'time ' /), & 1562 1623 output_type = 'real32' ) 1563 1624 1564 return_value = dom_def_att( nc_filename, &1565 variable_name = variable_name, &1566 attribute_name = 'units', &1567 value = ' W' )1568 1569 variable_name = ' rotor_power'1570 return_value = dom_def_var( nc_filename, &1571 variable_name = variable_name, &1572 dimension_names = (/ 'turbine', 'time ' /), 1625 return_value = dom_def_att( nc_filename, & 1626 variable_name = variable_name, & 1627 attribute_name = 'units', & 1628 value = 'degrees' ) 1629 1630 variable_name = 'yaw_angle' 1631 return_value = dom_def_var( nc_filename, & 1632 variable_name = variable_name, & 1633 dimension_names = (/ 'turbine', 'time ' /), & 1573 1634 output_type = 'real32' ) 1574 1635 1575 return_value = dom_def_att( nc_filename, & 1576 variable_name = variable_name, & 1577 attribute_name = 'units', & 1578 value = 'W' ) 1579 1580 variable_name = 'rotor_thrust' 1581 return_value = dom_def_var( nc_filename, & 1582 variable_name = variable_name, & 1583 dimension_names = (/ 'turbine', 'time ' /), & 1584 output_type = 'real32' ) 1585 1586 return_value = dom_def_att( nc_filename, & 1587 variable_name = variable_name, & 1588 attribute_name = 'units', & 1589 value = 'N' ) 1590 1591 1592 variable_name = 'wind_direction' 1593 return_value = dom_def_var( nc_filename, & 1594 variable_name = variable_name, & 1595 dimension_names = (/ 'turbine', 'time ' /), & 1596 output_type = 'real32' ) 1597 1598 return_value = dom_def_att( nc_filename, & 1599 variable_name = variable_name, & 1600 attribute_name = 'units', & 1636 return_value = dom_def_att( nc_filename, & 1637 variable_name = variable_name, & 1638 attribute_name = 'units', & 1601 1639 value = 'degrees' ) 1602 1603 variable_name = 'yaw_angle' 1604 return_value = dom_def_var( nc_filename, & 1605 variable_name = variable_name, & 1606 dimension_names = (/ 'turbine', 'time ' /), & 1607 output_type = 'real32' ) 1608 1609 return_value = dom_def_att( nc_filename, & 1610 variable_name = variable_name, & 1611 attribute_name = 'units', & 1612 value = 'degrees' ) 1613 1614 ENDIF 1640 1641 ENDIF 1615 1642 END SUBROUTINE 1616 1643 … … 1897 1924 !-- the rotor. If tilt = 0, rot_coord_trans is a unit matrix. 1898 1925 1899 rot_coord_trans(inot,1,1) = rot_eigen_rad(1)**2 * &1926 rot_coord_trans(inot,1,1) = rot_eigen_rad(1)**2 * & 1900 1927 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1901 rot_coord_trans(inot,1,2) = rot_eigen_rad(1) * rot_eigen_rad(2) * &1928 rot_coord_trans(inot,1,2) = rot_eigen_rad(1) * rot_eigen_rad(2) * & 1902 1929 ( 1.0_wp - COS( tilt_angle ) ) - & 1903 1930 rot_eigen_rad(3) * SIN( tilt_angle ) 1904 rot_coord_trans(inot,1,3) = rot_eigen_rad(1) * rot_eigen_rad(3) * &1931 rot_coord_trans(inot,1,3) = rot_eigen_rad(1) * rot_eigen_rad(3) * & 1905 1932 ( 1.0_wp - COS( tilt_angle ) ) + & 1906 1933 rot_eigen_rad(2) * SIN( tilt_angle ) 1907 rot_coord_trans(inot,2,1) = rot_eigen_rad(2) * rot_eigen_rad(1) * &1934 rot_coord_trans(inot,2,1) = rot_eigen_rad(2) * rot_eigen_rad(1) * & 1908 1935 ( 1.0_wp - COS( tilt_angle ) ) + & 1909 1936 rot_eigen_rad(3) * SIN( tilt_angle ) 1910 rot_coord_trans(inot,2,2) = rot_eigen_rad(2)**2 * &1937 rot_coord_trans(inot,2,2) = rot_eigen_rad(2)**2 * & 1911 1938 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1912 rot_coord_trans(inot,2,3) = rot_eigen_rad(2) * rot_eigen_rad(3) * &1939 rot_coord_trans(inot,2,3) = rot_eigen_rad(2) * rot_eigen_rad(3) * & 1913 1940 ( 1.0_wp - COS( tilt_angle ) ) - & 1914 1941 rot_eigen_rad(1) * SIN( tilt_angle ) 1915 rot_coord_trans(inot,3,1) = rot_eigen_rad(3) * rot_eigen_rad(1) * &1942 rot_coord_trans(inot,3,1) = rot_eigen_rad(3) * rot_eigen_rad(1) * & 1916 1943 ( 1.0_wp - COS( tilt_angle ) ) - & 1917 1944 rot_eigen_rad(2) * SIN( tilt_angle ) 1918 rot_coord_trans(inot,3,2) = rot_eigen_rad(3) * rot_eigen_rad(2) * &1945 rot_coord_trans(inot,3,2) = rot_eigen_rad(3) * rot_eigen_rad(2) * & 1919 1946 ( 1.0_wp - COS( tilt_angle ) ) + & 1920 1947 rot_eigen_rad(1) * SIN( tilt_angle ) 1921 rot_coord_trans(inot,3,3) = rot_eigen_rad(3)**2 * &1948 rot_coord_trans(inot,3,3) = rot_eigen_rad(3)**2 * & 1922 1949 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1923 1950 … … 2409 2436 2410 2437 IF ( tip_loss_correction ) THEN 2411 2438 ! 2412 2439 !-- Tip loss correction following Schito 2413 2440 !-- Schito applies the tip loss correction only to the lift force 2414 2441 !-- Therefore, the tip loss correction is only applied to the lift 2415 2442 !-- coefficient and not to the drag coefficient in our case 2416 !-- 2443 !-- 2417 2444 IF ( phi_rel(rseg) == 0.0_wp ) THEN 2418 2445 tl_factor = 1.0_wp 2419 2446 ELSE 2420 tl_factor = ( 2.0 / pi ) * &2421 ACOS( EXP( -1.0 * ( 3.0 * ( rotor_radius(inot) - cur_r ) / 2447 tl_factor = ( 2.0 / pi ) * & 2448 ACOS( EXP( -1.0 * ( 3.0 * ( rotor_radius(inot) - cur_r ) / & 2422 2449 ( 2.0 * cur_r * abs( sin( phi_rel(rseg) ) ) ) ) ) ) 2423 2450 ENDIF 2424 2425 turb_cl(rseg) = tl_factor * turb_cl(rseg) 2426 2427 ENDIF 2451 2452 turb_cl(rseg) = tl_factor * turb_cl(rseg) 2453 2454 ENDIF 2428 2455 ! 2429 2456 !-- !-----------------------------------------------------! … … 2464 2491 !$OMP CRITICAL 2465 2492 thrust_rotor(inot) = thrust_rotor(inot) + & 2466 thrust_seg(rseg) 2467 2493 thrust_seg(rseg) 2494 2468 2495 2469 2496 torque_total(inot) = torque_total(inot) + (torque_seg * cur_r) … … 2486 2513 CALL cpu_log( log_point_s(62), 'wtm_controller', 'start' ) 2487 2514 2488 2515 2489 2516 IF ( speed_control ) THEN 2490 2517 ! 2491 2518 !-- Calculation of the current generator speed for rotor speed control 2492 2493 ! 2519 2520 ! 2494 2521 !-- The acceleration of the rotor speed is calculated from 2495 2522 !-- the force balance of the accelerating torque … … 2504 2531 !-- and rotor speed 2505 2532 generator_speed(inot) = gear_ratio * ( rotor_speed(inot) + om_rate ) 2506 2533 2507 2534 ENDIF 2508 2535 2509 2536 IF ( pitch_control ) THEN 2510 2537 … … 2523 2550 CYCLE pit_loop 2524 2551 ENDIF 2525 2552 2526 2553 ! 2527 2554 !-- The current pitch is saved for the next time step … … 2535 2562 ! 2536 2563 !-- Call the rotor speed controller 2537 2564 2538 2565 IF ( speed_control ) THEN 2539 2566 ! … … 2546 2573 ENDIF 2547 2574 ENDIF 2548 2575 2549 2576 ENDIF 2550 2577 … … 2592 2619 !-- Predetermine flag to mask topography 2593 2620 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 2594 2621 2595 2622 ! 2596 2623 !-- Determine the square of the distance between the … … 2636 2663 ENDDO ! End of loop over rseg 2637 2664 ENDDO ! End of loop over ring 2638 2665 2639 2666 ! 2640 2667 !-- Rotation of force components: … … 2644 2671 torque_z(k,j,i)*rotz(inot,1) & 2645 2672 ) * flag 2646 2673 2647 2674 rot_tend_y(k,j,i) = rot_tend_y(k,j,i) + ( & 2648 2675 thrust(k,j,i)*rotx(inot,2) + & … … 2650 2677 torque_z(k,j,i)*rotz(inot,2) & 2651 2678 ) * flag 2652 2679 2653 2680 rot_tend_z(k,j,i) = rot_tend_z(k,j,i) + ( & 2654 2681 thrust(k,j,i)*rotx(inot,3) + & … … 2663 2690 2664 2691 CALL cpu_log( log_point_s(63), 'wtm_smearing', 'stop' ) 2665 2692 2666 2693 ENDDO !-- end of loop over turbines 2667 2694 2668 2695 2669 2696 IF ( yaw_control ) THEN 2670 2697 ! … … 2676 2703 wd30 = 999.0_wp ! Set to dummy value 2677 2704 ALLOCATE( wd30_l(1:WDLON) ) 2678 2705 2679 2706 WDSHO = MAX( 1 , NINT( 2.0_wp / dt_3d ) ) ! 2s running mean array 2680 2707 ALLOCATE( wd2(1:n_turbines,1:WDSHO) ) … … 2682 2709 ALLOCATE( wd2_l(1:WDSHO) ) 2683 2710 start_up = .FALSE. 2684 ENDIF 2711 ENDIF 2685 2712 2686 2713 ! … … 2707 2734 2708 2735 yaw_angle_l(inot) = yaw_angle(inot) 2709 2736 2710 2737 ENDIF 2711 2738 ENDIF 2712 2739 2713 2740 ENDDO !-- end of loop over turbines 2714 2741 … … 2726 2753 wdir = wdir_l 2727 2754 yaw_angle = yaw_angle_l 2728 2729 2755 2730 2756 #endif 2731 2757 DO inot = 1, n_turbines 2732 ! 2733 !-- Update rotor orientation 2758 ! 2759 !-- Update rotor orientation 2734 2760 CALL wtm_rotate_rotor( inot ) 2735 2761 2736 2762 ENDDO ! End of loop over turbines 2737 2763 2738 2764 ENDIF ! end of yaw control 2739 2765 2740 2766 IF ( speed_control ) THEN 2741 2767 ! … … 2755 2781 generator_speed_f_old = generator_speed_f 2756 2782 #endif 2757 2783 2758 2784 ENDIF 2759 2785 2760 2761 2762 2763 DO inot = 1, n_turbines 2764 2765 2766 2767 IF ( myid == 0 ) THEN 2768 IF ( openfile_turb_mod(400+inot)%opened ) THEN 2769 WRITE ( 400+inot, 106 ) time_since_reference_point, rotor_speed(inot), & 2770 generator_speed(inot), torque_gen_old(inot), & 2771 torque_total(inot), pitch_angle(inot), & 2772 torque_gen_old(inot)*generator_speed(inot)*generator_efficiency, & 2773 torque_total(inot)*rotor_speed(inot)*air_density, & 2774 thrust_rotor(inot), & 2775 wdir(inot)*180.0_wp/pi, & 2776 (yaw_angle(inot))*180.0_wp/pi 2777 2778 ELSE 2779 2780 WRITE ( turbine_id,'(A2,I2.2)') '_T', inot 2781 OPEN ( 400+inot, FILE=( 'WTM_OUTPUT_DATA' // & 2782 TRIM( coupling_char ) // & 2783 turbine_id ), FORM='FORMATTED' ) 2784 WRITE ( 400+inot, 105 ) inot 2785 WRITE ( 400+inot, 106 ) time_since_reference_point, rotor_speed(inot), & 2786 generator_speed(inot), torque_gen_old(inot), & 2787 torque_total(inot), pitch_angle(inot), & 2788 torque_gen_old(inot)*generator_speed(inot)*generator_efficiency, & 2789 torque_total(inot)*rotor_speed(inot)*air_density, & 2790 thrust_rotor(inot), & 2791 wdir(inot)*180.0_wp/pi, & 2792 (yaw_angle(inot))*180.0_wp/pi 2793 ENDIF 2794 ENDIF 2795 2796 !-- Set open flag 2797 openfile_turb_mod(400+inot)%opened = .TRUE. 2798 ENDDO !-- end of loop over turbines 2799 2800 2801 2802 ENDIF 2803 2804 2786 ENDIF 2787 2805 2788 CALL cpu_log( log_point_s(61), 'wtm_forces', 'stop' ) 2806 2807 ! 2808 !-- Formats 2809 105 FORMAT ('Turbine control data for turbine ',I2,1X,':'/ & 2810 &'----------------------------------------'/ & 2811 &' Time RSpeed GSpeed GenTorque AeroTorque ', & 2812 'Pitch Power(Gen) Power(Rot) RotThrust WDirection ', & 2813 'YawOrient') 2814 2815 106 FORMAT (F9.2,2X,F7.3,2X,F7.2,2X,F12.1,3X,F12.1,1X,F6.2,2X,F12.1,2X, & 2816 F12.1,1X,F12.1,4X,F7.2,4X,F7.2) 2817 2818 2819 2789 2790 2820 2791 END SUBROUTINE wtm_forces 2821 2792 2822 2793 2823 2794 !------------------------------------------------------------------------------! 2824 2795 ! Description: … … 2827 2798 !------------------------------------------------------------------------------! 2828 2799 SUBROUTINE wtm_yawcontrol( inot ) 2829 2800 2830 2801 USE kinds 2831 2802 2832 2803 IMPLICIT NONE 2833 2804 2834 2805 INTEGER(iwp) :: inot 2835 2806 INTEGER(iwp) :: i_wd_30 … … 2876 2847 ENDIF 2877 2848 ENDIF 2878 2849 2879 2850 wd30(inot,:) = wd30_l 2880 2851 2881 ! 2852 ! 2882 2853 !-- If turbine is already yawing: 2883 2854 !-- Initialize 2 s running mean and yaw until the missalignment is smaller … … 2890 2861 wd2_l = CSHIFT( wd2_l, SHIFT = -1 ) 2891 2862 wd2_l(1) = wdir(inot) 2892 ! 2863 ! 2893 2864 !-- Check if array is full ( no more dummies ) 2894 2865 IF ( .NOT. ANY( wd2_l == 999.0_wp ) ) THEN … … 2915 2886 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2916 2887 ENDIF 2917 2888 2918 2889 wd2(inot,:) = wd2_l 2919 2890 2920 2891 ENDIF 2921 2892 2922 2893 END SUBROUTINE wtm_yawcontrol 2923 2894 … … 2976 2947 2977 2948 INTEGER(iwp) :: inot 2978 2979 2949 2980 2950 2981 2951 ! … … 2992 2962 2993 2963 IF ( generator_speed_f(inot) <= region_15_min ) THEN 2994 ! 2964 ! 2995 2965 !-- Region 1: Generator torque is set to zero to accelerate the rotor: 2996 2966 torque_gen(inot) = 0 2997 2967 2998 2968 ELSEIF ( generator_speed_f(inot) <= region_2_min ) THEN 2999 ! 2969 ! 3000 2970 !-- Region 1.5: Generator torque is increasing linearly with rotor speed: 3001 2971 torque_gen(inot) = slope15 * ( generator_speed_f(inot) - region_15_min ) 3002 2972 3003 2973 ELSEIF ( generator_speed_f(inot) <= region_25_min ) THEN 3004 2974 ! … … 3006 2976 !-- speed to keep the TSR optimal: 3007 2977 torque_gen(inot) = region_2_slope * generator_speed_f(inot) * generator_speed_f(inot) 3008 2978 3009 2979 ELSEIF ( generator_speed_f(inot) < generator_speed_rated ) THEN 3010 ! 2980 ! 3011 2981 !-- Region 2.5: Slipage region between 2 and 3: 3012 2982 torque_gen(inot) = region_25_slope * ( generator_speed_f(inot) - vs_sysp ) 3013 2983 3014 2984 ELSE 3015 ! 2985 ! 3016 2986 !-- Region 3: Generator torque is antiproportional to the rotor speed to 3017 2987 !-- keep the power constant: 3018 2988 torque_gen(inot) = generator_power_rated / generator_speed_f(inot) 3019 2989 3020 2990 ENDIF 3021 ! 2991 ! 3022 2992 !-- Calculate torque rate and confine with a max 3023 2993 trq_rate = ( torque_gen(inot) - torque_gen_old(inot) ) / dt_3d 3024 2994 trq_rate = MIN( MAX( trq_rate, -1.0_wp * generator_torque_rate_max ), & 3025 2995 generator_torque_rate_max ) 3026 ! 3027 !-- Calculate new gen torque and confine with max torque 2996 ! 2997 !-- Calculate new gen torque and confine with max torque 3028 2998 torque_gen(inot) = torque_gen_old(inot) + trq_rate * dt_3d 3029 torque_gen(inot) = MIN( torque_gen(inot), generator_torque_max ) 3030 ! 3031 !-- Overwrite values for next timestep 2999 torque_gen(inot) = MIN( torque_gen(inot), generator_torque_max ) 3000 ! 3001 !-- Overwrite values for next timestep 3032 3002 rotor_speed_l(inot) = generator_speed(inot) / gear_ratio 3033 3003 3034 3004 3035 3005 END SUBROUTINE wtm_speed_control 3036 3006 … … 3090 3060 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3091 3061 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * & 3092 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3062 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3093 3063 tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * & 3094 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3064 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3095 3065 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) & 3096 3066 - tend_nac_y - tend_tow_y ) & … … 3156 3126 !-- Calculate the thrust generated by the nacelle and the tower 3157 3127 tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) * & 3158 SIGN( u(k,j,i)**2 , u(k,j,i) ) 3128 SIGN( u(k,j,i)**2 , u(k,j,i) ) 3159 3129 tend_tow_x = 0.5_wp * tow_cd_surf(k,j,i) * & 3160 SIGN( u(k,j,i)**2 , u(k,j,i) ) 3130 SIGN( u(k,j,i)**2 , u(k,j,i) ) 3161 3131 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i) & 3162 3132 - tend_nac_x - tend_tow_x ) & … … 3172 3142 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3173 3143 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * & 3174 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3144 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3175 3145 tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * & 3176 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3146 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3177 3147 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) & 3178 3148 - tend_nac_y - tend_tow_y ) & … … 3217 3187 ALLOCATE ( output_values_1d_target(1:n_turbines) ) 3218 3188 output_values_1d_target = hub_x(1:n_turbines) 3219 output_values_1d_pointer => output_values_1d_target 3189 output_values_1d_pointer => output_values_1d_target 3220 3190 return_value = dom_write_var( nc_filename, & 3221 3191 'x', & … … 3225 3195 3226 3196 output_values_1d_target = hub_y(1:n_turbines) 3227 output_values_1d_pointer => output_values_1d_target 3197 output_values_1d_pointer => output_values_1d_target 3228 3198 return_value = dom_write_var( nc_filename, & 3229 3199 'y', & … … 3233 3203 3234 3204 output_values_1d_target = hub_z(1:n_turbines) 3235 output_values_1d_pointer => output_values_1d_target 3205 output_values_1d_pointer => output_values_1d_target 3236 3206 return_value = dom_write_var( nc_filename, & 3237 3207 'z', & 3238 3208 values_realwp_1d = output_values_1d_pointer, & 3239 3209 bounds_start = (/1/), & 3240 bounds_end = (/n_turbines/) ) 3241 3210 bounds_end = (/n_turbines/) ) 3211 3212 output_values_1d_target = rotor_radius(1:n_turbines) * 2.0_wp 3213 output_values_1d_pointer => output_values_1d_target 3214 return_value = dom_write_var( nc_filename, & 3215 'rotor_diameter', & 3216 values_realwp_1d = output_values_1d_pointer, & 3217 bounds_start = (/1/), & 3218 bounds_end = (/n_turbines/) ) 3219 3220 output_values_1d_target = tower_diameter(1:n_turbines) 3221 output_values_1d_pointer => output_values_1d_target 3222 return_value = dom_write_var( nc_filename, & 3223 'tower_diameter', & 3224 values_realwp_1d = output_values_1d_pointer, & 3225 bounds_start = (/1/), & 3226 bounds_end = (/n_turbines/) ) 3227 3242 3228 initial_write_coordinates = .TRUE. 3243 3229 DEALLOCATE ( output_values_1d_target ) 3244 3230 ENDIF 3245 3231 3246 3232 t_ind = t_ind + 1 3247 3233
Note: See TracChangeset
for help on using the changeset viewer.
