diff --git a/src/MNH/eol_adr.f90 b/src/MNH/eol_adr.f90 new file mode 100644 index 0000000000000000000000000000000000000000..62f27088dc2c646855de191ba2e0f22da81e3fd3 --- /dev/null +++ b/src/MNH/eol_adr.f90 @@ -0,0 +1,520 @@ +!MNH_LIC Copyright 2017-2021 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!----------------------------------------------------------------- +! ####################### + MODULE MODI_EOL_ADR +! ####################### +! +INTERFACE +! +SUBROUTINE EOL_ADR(KTCOUNT, PTSTEP, & + PDXX, PDYY, PDZZ, & + PRHO_M, & + PUT_M, PVT_M, PWT_M, & + PFX_RG, PFY_RG, PFZ_RG ) + +! +INTEGER, INTENT(IN) :: KTCOUNT ! iteration count +REAL, INTENT(IN) :: PTSTEP ! timestep except +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ ! mesh size +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_M ! dry Density +REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT_M,PVT_M,PWT_M ! Wind speed at mass point +! +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFX_RG ! Aerodynamic force .. +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFY_RG ! .. cartesian mesh .. +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFZ_RG ! .. global frame) +! +! +END SUBROUTINE EOL_ADR +! +END INTERFACE +! +END MODULE MODI_EOL_ADR +! +! ################################################################### + SUBROUTINE EOL_ADR(KTCOUNT, PTSTEP, & + PDXX, PDYY, PDZZ, & + PRHO_M, & + PUT_M, PVT_M, PWT_M, & + PFX_RG, PFY_RG, PFZ_RG ) +! ################################################################### +! +!!**** *MODI_EOL_ALM* - +!! +!! PURPOSE +!! ------- +!! It is possible to include wind turbines parameterization in Meso-NH, +!! and several methods are available. One of the models is the Actuator +!! Line Method (ALM). It allows to compute aerodynamic forces according +!! the wind speed and the caracteristics of the wind turbine. +!! +!!** METHOD +!! ------ +!! The ALM consists in modeling each blade by one line divided into blade +!! element points (Sørensen and Shen, 2002). These points are applying +!! aerodynamic forces into the flow. +!! Each point carries a two-dimensional (2D) airfoil, and its characteris- +!! tics, as lift and drag coefficients. Knowing these coefficients, and +!! the angle of attack, the lift and drag forces can be evaluated. +!! +!! REFERENCE +!! --------- +!! PA. Joulin PhD Thesis. 2020. +!! +!! +!! AUTHOR +!! ------ +!! PA. Joulin *CNRM & IFPEN* +!! +!! MODIFICATIONS +!! ------------- +!! Original 24/01/17 +!! Modification 14/10/20 (PA. Joulin) Updated for a main version +!! +!!--------------------------------------------------------------- +! +! +!* 0. DECLARATIONS +! ------------ +! +! To work with wind turbines +USE MODD_EOL_ADR +USE MODD_EOL_KINE_ADR +! +USE MODD_EOL_SHARED_IO, ONLY: CINTERP, LTECOUTPTS, LTIPLOSSG +USE MODD_EOL_SHARED_IO, ONLY: XTHRUT, XTORQT, XPOWT +USE MODD_EOL_SHARED_IO, ONLY: XELT_RAD +USE MODD_EOL_SHARED_IO, ONLY: XAOA_GLB, XFLIFT_GLB, XFDRAG_GLB, XFAERO_RG_GLB +! +USE MODI_EOL_MATHS +USE MODI_EOL_READER, ONLY: GET_AIRFOIL_ID_ADR +USE MODI_EOL_PRINTER, ONLY: PRINT_TSPLIT +USE MODI_EOL_DEBUGGER +! Math +USE MODD_CST, ONLY: XPI +! To know the grid +USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZS,XZZ +USE MODE_ll, ONLY: GET_INDICE_ll +USE MODD_PARAMETERS, ONLY: JPVEXT +! MPI stuffs +USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD +USE MODD_PRECISION, ONLY: MNHREAL_MPI +USE MODD_MPIF, ONLY: MPI_SUM +USE MODE_SUM_ll, ONLY: MIN_ll +USE MODD_VAR_ll, ONLY: IP +! +! +IMPLICIT NONE +! +!* 0.1 Declarations of dummy arguments : +! +INTEGER, INTENT(IN) :: KTCOUNT ! iteration count +REAL, INTENT(IN) :: PTSTEP ! timestep except +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ ! mesh size +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_M ! dry Density +REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT_M,PVT_M,PWT_M ! Wind speed at mass point +! +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFX_RG ! Aerodynamic force .. +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFY_RG ! .. cartesian mesh .. +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFZ_RG ! .. global frame) +! +!* 0.2 Declarations of local variables : +! +! Indicies Compteurs +INTEGER :: IIB,IJB,IKB ! Begin of a CPU domain +INTEGER :: IIE,IJE,IKE ! End of a CPU domain +INTEGER :: IKU ! Vertical size of the domain +INTEGER :: JI, JJ, JK ! Loop index +INTEGER :: JROT, JAZI, JRAD ! Rotor, azimutal, and radial element indicies +! + +! -- Wind -- +REAL :: ZRHO_I ! Interpolated density [kg/m3] +REAL :: ZUT_I ! Interpolated wind speed U (RG) [m/s] +REAL :: ZVT_I ! Interpolated wind speed V (RG) [m/s] +REAL :: ZWT_I ! Interpolated wind speed W (RG) [m/s] +REAL, DIMENSION(3) :: ZWIND_VEL_RG ! Wind velocity in RG frame [m/s] +REAL, DIMENSION(3) :: ZWIND_VEL_RA ! Wind velocity in RE frame [m/s] +REAL, DIMENSION(3) :: ZWINDREL_VEL_RA ! Relative wind velocity in RE frame [m/s] +REAL :: ZWINDREL_VEL ! Norm of the relative wind velocity [m/s] +REAL, DIMENSION(SIZE(PUT_M,1),SIZE(PUT_M,2),SIZE(PUT_M,3)) :: ZZH ! True heigth to interpolate 8NB +! +! -- Wind turbine -- +INTEGER :: INB_WT, INB_B ! Total numbers of turbines and blades +INTEGER :: INB_AELT, INB_RELT ! Total numbers of azimutal and radial elt +REAL :: ZRAD ! Blade radius [m] +INTEGER :: IAID ! Airfoil index [-] +! +! -- Aero -- +REAL :: ZAOA ! Attack angle of an element [rad] +REAL :: ZCDRAG ! Drag coefficient of an element [] +REAL :: ZCLIFT ! Lift coefficient of an element [] +REAL :: ZFDRAG ! Drag force of an element, parallel to Urel [N] +REAL :: ZFLIFT ! Lift force of an element, perpendicular to Urel [N] +REAL, DIMENSION(3) :: ZFAERO_RA ! Aerodynamic force (lift+drag) in RA [N] +REAL, DIMENSION(3) :: ZFAERO_RG ! Aerodynamic force (lift+drag) in RG [N] +! Tip loss +REAL :: ZFTIPL ! tip loss function +REAL :: ZPHI ! angle twist+pitch+aa +! +! Thrust, Torque and Power +REAL, DIMENSION(3) :: ZFAERO_RH ! Aerodynamic force (lift+drag) in RH [N] (thrust/torque) +REAL, DIMENSION(3) :: ZDIST_HBELT_RH ! Distance between blade element and hub, in RH [m] +REAL, DIMENSION(3) :: ZDIST_HBELT_RG ! Distances between blade element and hub, in RG [m] +REAL, DIMENSION(3) :: Z3D_TORQT ! Full torque force (3D) of the wind turbine [N] +! +! -- Numerical -- +INTEGER :: IINFO ! code info return +! +! +!* 0.3 Implicit arguments +! +! A. From MODD_EOL_ALM +!TYPE(FARM) :: TFARM +!TYPE(TURBINE) :: TTURBINE +!TYPE(BLADE) :: TBLADE +!TYPE(AIRFOIL), DIMENSION(:), ALLOCATABLE :: TAIRFOIL +! +!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XELT_RAD ! Blade elements radius [m] +!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAOA_GLB ! Angle of attack of an element [rad] +!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFLIFT_GLB ! Lift force, parallel to Urel [N] +!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFDRAG_GLB ! Drag force, perpendicular to Urel [N] +!REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RE_GLB ! Aerodyn. force (lift+drag) in RE [N] +!REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RG_GLB ! Aerodyn. force (lift+drag) in RG [N] +! +!INTEGER :: NNB_BLAELT ! Number of blade elements +!LOGICAL :: LTIMESPLIT ! Flag to apply Time splitting method +!LOGICAL :: LTIPLOSSG ! Flag to apply Glauert's tip loss correction +!LOGICAL :: LTECOUTPTS ! Flag to get Tecplot file output of element points +! +! B. From MODD_EOL_SHARED_IO: +! for namelist NAM_EOL_ALM +!CHARACTER(LEN=100) :: CFARM_CSVDATA ! Farm file to read +!CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! Turbine file to read +!CHARACTER(LEN=100) :: CBLADE_CSVDATA ! Blade file to read +!CHARACTER(LEN=100) :: CAIRFOIL_CSVDATA ! Airfoil file to read +!CHARACTER(LEN=3) :: CINTERP ! Interpolation method for wind speed +! for output +!REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N] +!REAL, DIMENSION(:), ALLOCATABLE :: XTORQT ! Torque [Nm] +!REAL, DIMENSION(:), ALLOCATABLE :: XPOWT ! Power [W] +! +! +!------------------------------------------------------------------------------- +! +! +!* 1. INITIALIZATIONS +! --------------- +! +!* 1.1 Subdomain (CPU) indices +! +CALL GET_INDICE_ll(IIB,IJB,IIE,IJE) ! Get begin and end domain index (CPU) +IKU = SIZE(PUT_M,3) ! Top of the domain end index +IKB=1+JPVEXT ! Vertical begin index +IKE=IKU-JPVEXT ! Vertical end index +! +!* 1.2 Some usefull integers +! +INB_WT = TFARM%NNB_TURBINES +INB_B = TTURBINE%NNB_BLADES +INB_AELT = TBLADE%NNB_AZIELT +INB_RELT = TBLADE%NNB_RADELT +! +!* 1.3 Vertical coordinate in case of interpolation +! +IF (CINTERP=='8NB') THEN + DO JK=1,IKU-1 + ZZH(:,:,JK) = (0.5*(XZZ(:,:,JK)+XZZ(:,:,JK+1))-XZS(:,:)) + END DO + ZZH(:,:,IKU) = 2*ZZH(:,:,IKU-1) - ZZH(:,:,IKU-2) +END IF +! +!* 1.4 Set to zeros at each MNH time steps +! +XAOA_ELT(:,:,:) = 0. +XFLIFT_ELT(:,:,:) = 0. +XFDRAG_ELT(:,:,:) = 0. +XFAERO_RA_ELT(:,:,:,:) = 0. +XFAERO_RG_ELT(:,:,:,:) = 0. +! Global variables (seen by all CPU) +XAOA_GLB(:,:,:) = 0. +XFLIFT_GLB(:,:,:) = 0. +XFDRAG_GLB(:,:,:) = 0. +XFAERO_RA_GLB(:,:,:,:) = 0. +XFAERO_RG_GLB(:,:,:,:) = 0. +! +XFAERO_RA_BLA(:,:,:) = 0. +XAOA_BLA(:,:) = 0. +! +XTHRUT(:) = 0. +XTORQT(:) = 0. +! +!------------------------------------------------------------------------------ +! +!* 2. KINEMATICS COMPUTATIONS +! ----------------------- +! + CALL EOL_KINE_ADR(KTCOUNT, PTSTEP) +! +! +!------------------------------------------------------------------------------- +! +!* 4. COMPUTES AERODYNAMIC FORCES THAT ACTS ON THE BLADES DUE TO THE WIND +! -------------------------------------------------------------- +! +!* 4.1 Finding the position of wind turbines +! +! Loop over domain + DO JK=IKB,IKE + DO JJ=IJB,IJE + DO JI=IIB,IIE + ! Loop over wind turbines + DO JROT=1, INB_WT + DO JAZI=1, INB_AELT + DO JRAD=1, INB_RELT + + ! Position test + IF (XPOS_ELT_RG(JROT,JAZI,JRAD,1) >= XXHAT(JI) .AND. & + XPOS_ELT_RG(JROT,JAZI,JRAD,1) < XXHAT(JI) + PDXX(JI,JJ,JK)) THEN +! + IF (XPOS_ELT_RG(JROT,JAZI,JRAD,2) >= XYHAT(JJ) .AND. & + XPOS_ELT_RG(JROT,JAZI,JRAD,2) < XYHAT(JJ) + PDYY(JI,JJ,JK)) THEN +! + IF (XPOS_ELT_RG(JROT,JAZI,JRAD,3) >= XZZ(JI,JJ,JK) .AND. & + XPOS_ELT_RG(JROT,JAZI,JRAD,3) < XZZ(JI,JJ,JK) + PDZZ(JI,JJ,JK)) THEN +! +!* 4.2 Extracting the wind +! + SELECT CASE(CINTERP) + CASE('CLS') + ZUT_I = PUT_M(JI,JJ,JK) + ZVT_I = PVT_M(JI,JJ,JK) + ZWT_I = PWT_M(JI,JJ,JK) + ZRHO_I = PRHO_M(JI,JJ,JK) + CASE('8NB') + ZUT_I = INTERP_LIN8NB(XPOS_ELT_RG(JROT,JAZI,JRAD,:),& + JI,JJ,JK,PUT_M,ZZH) + ZVT_I = INTERP_LIN8NB(XPOS_ELT_RG(JROT,JAZI,JRAD,:),& + JI,JJ,JK,PVT_M,ZZH) + ZWT_I = INTERP_LIN8NB(XPOS_ELT_RG(JROT,JAZI,JRAD,:),& + JI,JJ,JK,PWT_M,ZZH) + ZRHO_I = INTERP_LIN8NB(XPOS_ELT_RG(JROT,JAZI,JRAD,:),& + JI,JJ,JK,PRHO_M,ZZH) + END SELECT + ZWIND_VEL_RG(1) = ZUT_I + ZWIND_VEL_RG(2) = ZVT_I + ZWIND_VEL_RG(3) = ZWT_I +! +! PRINT*, 'JROT= ', JROT, 'JAZI= ', JAZI, 'JRAD= ', JRAD +! PRINT*, 'WIND_VEL_RG(1)= ', ZWIND_VEL_RG(1),'WIND_VEL_RG(2)= ', ZWIND_VEL_RG(2), 'WIND_VEL_RG(3)= ', ZWIND_VEL_RG(3) +!* 4.3 Calculating the wind in RA frame +! + ZWIND_VEL_RA(:) = MATMUL(XMAT_RA_RG(JROT,JAZI,JRAD,:,:), ZWIND_VEL_RG(:)) +! PRINT*, 'WIND_VEL_RA(1)= ', ZWIND_VEL_RA(1),'WIND_VEL_RA(2)= ', ZWIND_VEL_RA(2), 'WIND_VEL_RA(3)= ', ZWIND_VEL_RA(3) +! +!* 4.4 Calculating the relative wind speed in RE frame + norm +! + ZWINDREL_VEL_RA(:) = ZWIND_VEL_RA(:) - XTVEL_ELT_RA(JROT,JAZI,JRAD,:) + ZWINDREL_VEL = NORM(ZWINDREL_VEL_RA) +! +! PRINT*, 'WINDREL_VEL_RA(1)= ', ZWINDREL_VEL_RA(1),'WINDREL_VEL_RA(2)= ', ZWINDREL_VEL_RA(2) +! PRINT*, 'WINDREL_VEL_RA(3)= ', ZWINDREL_VEL_RA(3) +! PRINT*, 'WINDREL_VEL ', ZWINDREL_VEL +!* 4.5 Calculating the angle of attack +! + ZAOA = ATAN2(-ZWINDREL_VEL_RA(3), ZWINDREL_VEL_RA(2)) +! +! PRINT*, 'AOA= ', ZAOA +!* 4.6 Getting aerodynamic coefficients from tabulated data +! + ZRAD = XELT_RAD(JROT,JAZI,JRAD) ! Radius of the element +! PRINT*, 'RAD= ', ZRAD + IAID = GET_AIRFOIL_ID_ADR(TTURBINE,TBLADE,TAIRFOIL,ZRAD) ! ID of the airfoil + ZCLIFT = INTERP_SPLCUB(ZAOA*180/XPI, & + TAIRFOIL(IAID)%XAA,& + TAIRFOIL(IAID)%XCL) + ZCDRAG = INTERP_SPLCUB(ZAOA*180/XPI, & + TAIRFOIL(IAID)%XAA,& + TAIRFOIL(IAID)%XCD) + +! PRINT*, 'CLIFT= ', ZCLIFT +! PRINT*, 'CDRAG= ', ZCDRAG + +! +!* 4.7 Tip loss correction (Glauert) +! + IF (LTIPLOSSG) THEN + ZPHI = + ZAOA & + - TFARM%XBLA_PITCH(JROT) & + - XTWIST_ELT(JROT,JRAD) + IF (ZPHI > 0.0) THEN + ZFTIPL = (2.0/XPI)*ACOS(MIN( & + 1.0, EXP(-(TTURBINE%NNB_BLADES/2.0) & + *(TTURBINE%XR_MAX-ZRAD)/(ZRAD*SIN(ZPHI))))) + ! PRINT*, 'JROT= ', JROT, 'JAZI= ', JAZI, 'JROT= ', JRAD + ! PRINT*, 'RAD= ', ZRAD, 'PHI= ', ZPHI, 'FTIPL=', ZFTIPL + ! PRINT*, 'EXP= ', EXP(-(TTURBINE%NNB_BLADES/2.0)*(TTURBINE%XR_MAX-ZRAD)/(ZRAD*SIN(ZPHI))) + ! PRINT*, 'ACOS= ', ACOS(MIN(1.0, EXP(-(TTURBINE%NNB_BLADES/2.0)*(TTURBINE%XR_MAX-ZRAD)/(ZRAD*SIN(ZPHI))))) + + ELSE + ZFTIPL = 1.0 + END IF + ZCLIFT = ZFTIPL*ZCLIFT + ZCDRAG = ZFTIPL*ZCDRAG + END IF +! +!* 4.8 Computing aerodynamic forces in relative frame +! that act on blades (wind->blade) + ZFLIFT = 0.5*ZRHO_I*XSURF_APP_ELT(JROT,JRAD)*ZCLIFT*ZWINDREL_VEL**2 + ZFDRAG = 0.5*ZRHO_I*XSURF_APP_ELT(JROT,JRAD)*ZCDRAG*ZWINDREL_VEL**2 +! PRINT*, 'LIFT= ', ZFLIFT, 'DRAG= ', ZFDRAG +! +!* 4.9 Evaluating the aerodynamiques forces in RE frame +! that act on blades (wind->blade) + ZFAERO_RA(1) = .0 + ZFAERO_RA(2) = COS(ZAOA)*ZFDRAG - SIN(ZAOA)*ZFLIFT + ZFAERO_RA(3) = - SIN(ZAOA)*ZFDRAG - COS(ZAOA)*ZFLIFT +! +! PRINT*, 'FAERO_RA(1)= ', ZFAERO_RA(1), 'FAERO_RA(2)= ', ZFAERO_RA(2) +!* 4.10 Evaluating the aerodynamiques forces in RG frame +! that act on blades (wind->blade) + ZFAERO_RG(:) = MATMUL(XMAT_RG_RA(JROT,JAZI,JRAD,:,:), ZFAERO_RA(:)) +! +!* 4.11 Adding it to the cell of Meso-NH + PFX_RG(JI,JJ,JK) = PFX_RG(JI,JJ,JK) + ZFAERO_RG(1) + PFY_RG(JI,JJ,JK) = PFY_RG(JI,JJ,JK) + ZFAERO_RG(2) + PFZ_RG(JI,JJ,JK) = PFZ_RG(JI,JJ,JK) + ZFAERO_RG(3) +! PRINT*, 'FX= ', PFX_RG(JI,JJ,JK), 'FY= ', PFY_RG(JI,JJ,JK), 'FZ= ', PFZ_RG(JI,JJ,JK) +! + + XAOA_ELT(JROT,JAZI,JRAD) = ZAOA + XFLIFT_ELT(JROT,JAZI,JRAD) = ZFLIFT + XFDRAG_ELT(JROT,JAZI,JRAD) = ZFDRAG + XFAERO_RA_ELT(JROT,JAZI,JRAD,:) = ZFAERO_RA(:) + XFAERO_RG_ELT(JROT,JAZI,JRAD,:) = ZFAERO_RG(:) + + + !PRINT*, '---------------------------------------------------- ' + !IF (IP==1) THEN + + !PRINT*, 'TCOUNT= ', KTCOUNT, 'STEP= ', PTSTEP + !PRINT*, 'JROT= ', JROT, 'JAZI= ', JAZI, 'JRAD= ', JRAD + !PRINT*, 'WIND_VEL_RG(1)= ', ZWIND_VEL_RG(1),'WIND_VEL_RG(2)= ', ZWIND_VEL_RG(2), 'WIND_VEL_RG(3)= ', ZWIND_VEL_RG(3) + !PRINT*, 'WIND_VEL_RA(1)= ', ZWIND_VEL_RA(1),'WIND_VEL_RA(2)= ', ZWIND_VEL_RA(2), 'WIND_VEL_RA(3)= ', ZWIND_VEL_RA(3) + !PRINT*, 'WINDREL_VEL_RA(1)= ', ZWINDREL_VEL_RA(1) + !PRINT*, 'WINDREL_VEL_RA(2)= ', ZWINDREL_VEL_RA(2) + !PRINT*, 'WINDREL_VEL_RA(3)= ', ZWINDREL_VEL_RA(3) + !PRINT*, 'WINDREL_VEL ', ZWINDREL_VEL + !PRINT*, 'AOA= ', ZAOA + !PRINT*, 'RAD= ', ZRAD + !PRINT*, 'CLIFT= ', ZCLIFT, 'CDRAG= ', ZCDRAG + !PRINT*, 'LIFT= ', ZFLIFT, 'DRAG= ', ZFDRAG + !PRINT*, 'FAERO_RA(1)= ', ZFAERO_RA(1), 'FAERO_RA(2)= ', ZFAERO_RA(2) + !PRINT*, 'FX= ', PFX_RG(JI,JJ,JK), 'FY= ', PFY_RG(JI,JJ,JK), 'FZ= ', PFZ_RG(JI,JJ,JK) + ! END IF + + !PRINT*, 'FAERO(1)= ', ZFAERO(JRAD,1) , 'FAERO(1)= ', ZFAERO(JRAD,2), 'FAERO(1)= ', ZFAERO(JRAD,3) + ! End of position tests + END IF + END IF + END IF + ! End of wind turbine loops + END DO + END DO + END DO + ! End of domain loops + END DO + END DO + END DO +! +!* 4.13 Top and bottom conditions +PFX_RG(:,:,IKB-1) = PFX_RG(:,:,IKB) +PFX_RG(:,:,IKE+1) = PFX_RG(:,:,IKE) +! +PFY_RG(:,:,IKB-1) = PFY_RG(:,:,IKB) +PFY_RG(:,:,IKE+1) = PFY_RG(:,:,IKE) +! +PFZ_RG(:,:,IKB-1) = PFZ_RG(:,:,IKB) +PFZ_RG(:,:,IKE+1) = PFZ_RG(:,:,IKE) +! +! +!------------------------------------------------------------------------------- +! +!* 5. SHARING THE DATAS OVER THE CPUS +! ------------------------------- +CALL MPI_ALLREDUCE(XAOA_ELT, XAOA_GLB, SIZE(XAOA_GLB), & + MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO) +CALL MPI_ALLREDUCE(XFLIFT_ELT, XFLIFT_GLB, SIZE(XFLIFT_GLB), & + MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO) +CALL MPI_ALLREDUCE(XFDRAG_ELT, XFDRAG_GLB, SIZE(XFDRAG_GLB), & + MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO) +CALL MPI_ALLREDUCE(XFAERO_RA_ELT, XFAERO_RA_GLB, SIZE(XFAERO_RA_GLB),& + MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO) +CALL MPI_ALLREDUCE(XFAERO_RG_ELT, XFAERO_RG_GLB, SIZE(XFAERO_RG_GLB),& + MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO) +! +! +!PRINT*, 'AOA_GLB= ', XAOA_GLB +!PRINT*, 'FLIFT_GLB= ', XFLIFT_GLB +!PRINT*, 'FDRAG_GLB= ', XFDRAG_GLB +!PRINT*, 'FAERO_RA_GLB= ', XFAERO_RA_GLB +!PRINT*, 'FAERO_RG_GLB= ', XFAERO_RG_GLB +!------------------------------------------------------------------------------- +! +!* 6. COMPUTING THRUST, TORQUE AND POWER +! --------------------------------- +! +IF(IP == 1) THEN + DO JROT=1,TFARM%NNB_TURBINES + + DO JRAD=1, TBLADE%NNB_RADELT + DO JAZI=1, TBLADE%NNB_AZIELT +! +!* 6.1 Preliminaries +! Mean AOA on one blade + XAOA_BLA(JROT,JRAD) = XAOA_BLA(JROT,JRAD) + XAOA_GLB(JROT,JAZI,JRAD)/INB_AELT + +! Aerodynamic load (wind->blade) in RA + XFAERO_RA_BLA(JROT,JRAD,:) = XFAERO_RA_BLA(JROT,JRAD,:) + XFAERO_RA_GLB(JROT,JAZI,JRAD,:)/INB_B + + +! PRINT*, '---------------------------------------------------- ' +! PRINT*, 'JROT= ', JROT, 'JAZI= ', JAZI, 'JRAD= ', JRAD +! PRINT*, 'FAERO_RA(1)= ', XFAERO_RA_BLA(JROT,JRAD,1),'FAERO_RA(2)= ', XFAERO_RA_BLA(JROT,JRAD,2) +! PRINT*, 'FAERO_RA(3)= ', XFAERO_RA_BLA(JROT,JRAD,3) +! PRINT*, 'AOA= ', XAOA_BLA(JROT,JRAD) + END DO + END DO + + + DO JAZI=1, TBLADE%NNB_AZIELT + DO JRAD=1, TBLADE%NNB_RADELT +! Aerodynamic load (wind->blade) in RH + ZFAERO_RH(:) = MATMUL(XMAT_RH_RG(JROT,:,:), & + XFAERO_RG_GLB(JROT,JAZI,JRAD,:)) +! Distance between element and hub in RG + ZDIST_HBELT_RG(:) = XPOS_ELT_RG(JROT,JAZI,JRAD,:) - XPOS_HUB_RG(JROT,:) +! Distance between element and hub in RH + ZDIST_HBELT_RH(:) = MATMUL(XMAT_RH_RG(JROT,:,:),ZDIST_HBELT_RG(:)) +! +!* 6.2 Thrust (wind->rotor): in RH + XTHRUT(JROT) = XTHRUT(JROT) + ZFAERO_RH(3) ! Only Z component +!* 6.3 Torque (wind->rotor) in RH + Z3D_TORQT = CROSS(ZDIST_HBELT_RH(:),ZFAERO_RH(:)) + XTORQT(JROT) = XTORQT(JROT) + Z3D_TORQT(3) ! Only Z component + END DO + END DO +! +!* 6.4 Power (wind->rotor) + XPOWT(JROT) = XTORQT(JROT) * TFARM%XOMEGA(JROT) + END DO +END IF +! +! +END SUBROUTINE EOL_ADR diff --git a/src/MNH/eol_kine_adr.f90 b/src/MNH/eol_kine_adr.f90 new file mode 100644 index 0000000000000000000000000000000000000000..7d2de65c7dee7f5c75147cef63b56c82aecff285 --- /dev/null +++ b/src/MNH/eol_kine_adr.f90 @@ -0,0 +1,344 @@ +!MNH_LIC Copyright 2017-2021 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!----------------------------------------------------------------- +! +!####################### +MODULE MODI_EOL_KINE_ADR +! +INTERFACE !---------------------------------------- +! +SUBROUTINE EOL_KINE_ADR(KTCOUNT,PTSTEP) +! + INTEGER, INTENT(IN) :: KTCOUNT ! iteration count + REAL, INTENT(IN) :: PTSTEP ! timestep +! +END SUBROUTINE EOL_KINE_ADR +! +END INTERFACE !------------------------------------ +! +END MODULE MODI_EOL_KINE_ADR +!####################### +! +! +! +! +!################################################################### +SUBROUTINE EOL_KINE_ADR(KTCOUNT,PTSTEP) +! +!!**** *EOL_KINEMATICS * - +!! +!! PURPOSE +!! ------- +!! Compute positions, oritentations and velocities of all the +!! elements of the wind turbine +!! +!!** METHOD +!! ------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! PA. Joulin *CNRM & IFPEN* +!! +!! MODIFICATIONS +!! ------------- +!! Original 04/2017 +!! Modification 10/11/20 (PA. Joulin) Updated for a main version +! P. Wautelet 19/07/2021: replace double precision by real to allow MNH_REAL=4 compilation +!! +!!--------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +!* 0.1 Modules +USE MODD_EOL_KINE_ADR +USE MODD_EOL_ADR +USE MODI_EOL_MATHS +USE MODI_EOL_PRINTER +USE MODD_TIME_n, ONLY : TDTCUR +USE MODD_CST, ONLY : XPI +USE MODD_EOL_SHARED_IO, ONLY: LTECOUTPTS,LCSVOUTFRM +! +IMPLICIT NONE +! +!* 0.2 Declarations of dummy arguments : +! +INTEGER, INTENT(IN) :: KTCOUNT ! iteration count +REAL, INTENT(IN) :: PTSTEP ! timestep +! +!* 0.3 Local variables +REAL, DIMENSION(3,3) :: ZORI_TMAT_X, ZORI_TMAT_Y, ZORI_TMAT_Z +REAL, DIMENSION(3,3) :: ZORI_MAT_X, ZORI_MAT_Y, ZORI_MAT_Z +REAL, DIMENSION(3) :: ZADD_TO_POS +! +REAL, DIMENSION(3) :: ZDIST_TOWO_TELT_RG ! Distance between tower elmt and tower base +REAL, DIMENSION(3) :: ZDIST_TOWO_NELT_RG ! Distance between nacelle and base of tower +REAL, DIMENSION(3) :: ZDIST_NAC_HUB_RG ! Distance between hub and base of nacelle +REAL, DIMENSION(3) :: ZDIST_HUB_CYL_RG ! Distance between cylindrical base and base of hub +REAL, DIMENSION(3) :: ZDIST_CYL_ELT_RG ! Distance between cylindrical base and elements +! +REAL :: ZTIME ! TIME +INTEGER :: JROT, JTELT, JNELT, JAZI ,JRAD ! Loop control +INTEGER :: INB_WT,INB_B,INB_AELT, INB_RELT ! Total numbers +INTEGER :: INB_TELT, INB_NELT ! Total numbers +INTEGER :: ITECOUT ! Unit number for Tecplot file +INTEGER :: ICSVOUT ! Unit number for TOW.csv file +! +! +!--------------------------------------------------------------- +! +!* 1. PRELIMINARIES +! ------------- +! +!* 1.1 Some useful integers +INB_WT = TFARM%NNB_TURBINES +INB_B = TTURBINE%NNB_BLADES +INB_TELT = 2 +INB_NELT = 2 +INB_AELT = TBLADE%NNB_AZIELT +INB_RELT = TBLADE%NNB_RADELT + +! +!* 1.2 Sub-time computation +ZTIME = TDTCUR%xtime +! +!* 1.3 Tecplotfile : opening + headers +IF (LTECOUTPTS) THEN + CALL OPEN_TECOUT_ADR(ITECOUT, KTCOUNT) +END IF +IF (LCSVOUTFRM) THEN + CALL OPEN_3DCSV_ADR(ICSVOUT, KTCOUNT) +END IF +! +! +!--------------------------------------------------------------- +! +!* 2. COMPUTATIONS +! ------------ +! +DO JROT=1, INB_WT +! +! ---- TOWER ---- +! +!* T.0 Update origin positions in RG (Base position + floating) + XPOS_TOWO_RG(JROT,:) = XPOS_REF(:) + XPOSINI_TOWO_RG(JROT,:) & + + XTVEL_TOWO_RG(JROT,:)*ZTIME +! +!* T.1 Update orientation + CALL GET_ORI_MAT_X(XANGINI_TOW_RG(JROT,1) + XRVEL_RT_RG(JROT,1)*ZTIME, ZORI_MAT_X) + CALL GET_ORI_MAT_Y(XANGINI_TOW_RG(JROT,2) + XRVEL_RT_RG(JROT,2)*ZTIME, ZORI_MAT_Y) + CALL GET_ORI_MAT_Z(XANGINI_TOW_RG(JROT,3) + XRVEL_RT_RG(JROT,3)*ZTIME, ZORI_MAT_Z) +! Compute orientation matrix + XMAT_RG_RT(JROT,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z)) +! +!* T.2 Update positions in RG + DO JTELT=1, INB_TELT + XPOS_TELT_RG(JROT,JTELT,:) = XPOS_TOWO_RG(JROT,:) & + + MATMUL(XMAT_RG_RT(JROT,:,:),XPOS_TELT_RT(JROT,JTELT,:)) + END DO +! +!* T.3 Update structural velocities + DO JTELT=1, INB_TELT + ! Rotation of tower already in RG + ! Translation of elements + ZDIST_TOWO_TELT_RG(:) = XPOS_TELT_RG(JROT,JTELT,:) - XPOS_TOWO_RG(JROT,:) + XTVEL_TELT_RG(JROT,JTELT,:) = XTVEL_TOWO_RG(JROT,:) & + + CROSS(XRVEL_RT_RG(JROT,:),ZDIST_TOWO_TELT_RG(:)) + ENDDO +! +!* T.4 Print in tecplot file + IF (LTECOUTPTS) THEN + DO JTELT=1, INB_TELT + CALL PRINT_TECOUT(ITECOUT, XPOS_TELT_RG(JROT,JTELT,:)) + END DO + END IF +! +!* T.5 Print in 3D csv file + IF (LCSVOUTFRM) THEN + CALL PRINT_3DFRM(ICSVOUT,'TOW',XMAT_RG_RT(JROT,:,:),XPOS_TOWO_RG(JROT,:)) + END IF +! +! +! ---- NACELLE ---- +! +!* N.0 Update origin positions in RG + XPOS_NACO_RG(JROT,:) = XPOS_TELT_RG(JROT,INB_TELT,:) & + + MATMUL(XMAT_RG_RT(JROT,:,:),XPOSINI_NACO_RT(JROT,:)) +! +!* N.1 Update orientation + CALL GET_ORI_MAT_X(XANGINI_NAC_RT(JROT,1) + XRVEL_RN_RT(JROT,1)*ZTIME, ZORI_MAT_X) + CALL GET_ORI_MAT_Y(XANGINI_NAC_RT(JROT,2) + XRVEL_RN_RT(JROT,2)*ZTIME, ZORI_MAT_Y) + CALL GET_ORI_MAT_Z(XANGINI_NAC_RT(JROT,3) + XRVEL_RN_RT(JROT,3)*ZTIME, ZORI_MAT_Z) +! +! Orientation matrix + XMAT_RT_RN(JROT,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z)) + XMAT_RG_RN(JROT,:,:) = MATMUL(XMAT_RG_RT(JROT,:,:), XMAT_RT_RN(JROT,:,:)) +! +!* N.2 Update positions in RG + DO JNELT=1, INB_NELT + XPOS_NELT_RG(JROT,JNELT,:) = XPOS_NACO_RG(JROT,:) & + + MATMUL(XMAT_RG_RN(JROT,:,:),XPOS_NELT_RN(JROT,JNELT,:)) + END DO +! +!* N.3 Update structural velocities + ! Rotation of nacelle in RG + XRVEL_RN_RG(JROT,:) = MATMUL(XMAT_RG_RT(JROT,:,:),XRVEL_RN_RT(JROT,:)) & + + XRVEL_RT_RG(JROT,:) + DO JNELT=1, INB_NELT + ! Translation of elements in RG + ZDIST_TOWO_NELT_RG(:) = XPOS_NELT_RG(JROT,JNELT,:) - XPOS_TOWO_RG(JROT,:) + XTVEL_NELT_RG(JROT,JNELT,:) = XTVEL_TOWO_RG(JROT,:) & + + CROSS(XRVEL_RN_RG(JROT,:),ZDIST_TOWO_NELT_RG(:)) + END DO +! +!* N.4 Print in tecplot file + IF (LTECOUTPTS) THEN + DO JNELT=1, INB_NELT + CALL PRINT_TECOUT(ITECOUT, XPOS_NELT_RG(JROT,JNELT,:)) + END DO + END IF +! +!* N.5 Print in 3D csv file + IF (LCSVOUTFRM) THEN + CALL PRINT_3DFRM(ICSVOUT,'NAC',XMAT_RG_RN(JROT,:,:),XPOS_NACO_RG(JROT,:)) + END IF +! +! ---- HUB ---- +! +!* H.1 Update positions + XPOS_HUB_RG(JROT,:) = XPOS_NELT_RG(JROT,INB_NELT,:) & + + MATMUL(XMAT_RG_RN(JROT,:,:),XPOSINI_HUB_RN(JROT,:)) +! +!* H.2 Update orientation + CALL GET_ORI_MAT_X(XANGINI_HUB_RN(JROT,1) + XRVEL_RH_RN(JROT,1)*ZTIME, ZORI_MAT_X) + CALL GET_ORI_MAT_Y(XANGINI_HUB_RN(JROT,2) + XRVEL_RH_RN(JROT,2)*ZTIME, ZORI_MAT_Y) + CALL GET_ORI_MAT_Z(XANGINI_HUB_RN(JROT,3) + XRVEL_RH_RN(JROT,3)*ZTIME, ZORI_MAT_Z) +! Orientation matrix + XMAT_RN_RH(JROT,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z)) + XMAT_RG_RH(JROT,:,:) = MATMUL(XMAT_RG_RN(JROT,:,:), XMAT_RN_RH(JROT,:,:)) + XMAT_RH_RG(JROT,:,:) = TRANSPOSE(XMAT_RG_RH(JROT,:,:)) +! +!* H.3 Update structural velocities +! Rotation of hub in RG + XRVEL_RH_RG(JROT,:) = MATMUL(XMAT_RG_RH(JROT,:,:),XRVEL_RH_RN(JROT,:)) & + + XRVEL_RN_RG(JROT,:) +! Translation of hub in RG + ZDIST_NAC_HUB_RG(:) = XPOS_HUB_RG(JROT,:) - XPOS_NELT_RG(JROT,INB_NELT,:) + XTVEL_HUB_RG(JROT,:) = XTVEL_NELT_RG(JROT,INB_NELT,:) + CROSS(XRVEL_RH_RG(JROT,:),ZDIST_NAC_HUB_RG(:)) +! +!* H.4 Print in tecplot file + IF (LTECOUTPTS) THEN + CALL PRINT_TECOUT(ITECOUT, XPOS_HUB_RG(JROT,:)) + END IF +! +!* H.5 Print in 3D csv file + IF (LCSVOUTFRM) THEN + CALL PRINT_3DFRM(ICSVOUT,'HUB',XMAT_RG_RH(JROT,:,:),XPOS_HUB_RG(JROT,:)) + END IF +! +! ---- CYLINDRICAL ---- +!* C.1 Update positions + XPOS_CYL_RG(JROT,:) = XPOS_HUB_RG(JROT,:) & + + MATMUL(XMAT_RG_RH(JROT,:,:),XPOSINI_CYL_RH(JROT,:)) +! +!* C.2 Update orientation + CALL GET_ORI_MAT_X(XANGINI_CYL_RH(JROT,1) + XRVEL_RC_RH(JROT,1)*ZTIME, ZORI_MAT_X) + CALL GET_ORI_MAT_Y(XANGINI_CYL_RH(JROT,2) + XRVEL_RC_RH(JROT,2)*ZTIME, ZORI_MAT_Y) + CALL GET_ORI_MAT_Z(XANGINI_CYL_RH(JROT,3) + XRVEL_RC_RH(JROT,3)*ZTIME, ZORI_MAT_Z) +! Orientation matrix + XMAT_RH_RC(JROT,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z)) + XMAT_RG_RC(JROT,:,:) = MATMUL(XMAT_RG_RH(JROT,:,:), XMAT_RH_RC(JROT,:,:)) +! +!* C.3 Update structural velocities +! Rotation of cylindrical in RG + XRVEL_RC_RG(JROT,:) = MATMUL(XMAT_RG_RC(JROT,:,:),XRVEL_RC_RH(JROT,:)) & + + XRVEL_RH_RG(JROT,:) +! Translation of hub in RG + ZDIST_HUB_CYL_RG(:) = XPOS_CYL_RG(JROT,:) - XPOS_HUB_RG(JROT,:) + XTVEL_CYL_RG(JROT,:) = XTVEL_HUB_RG(JROT,:) + CROSS(XRVEL_RC_RG(JROT,:),ZDIST_HUB_CYL_RG(:)) +! +!* C.4 Print in tecplot file + IF (LTECOUTPTS) THEN + CALL PRINT_TECOUT(ITECOUT, XPOS_CYL_RG(JROT,:)) + END IF +! +!* C.5 Print in 3D csv file + IF (LCSVOUTFRM) THEN + CALL PRINT_3DFRM(ICSVOUT,'CYL',XMAT_RG_RC(JROT,:,:),XPOS_CYL_RG(JROT,:)) + END IF +! +! ---- ANNULAR ELEMENTS ---- +!* A.0 Positioning sections (cuts) in RG + DO JAZI=1, INB_AELT + DO JRAD=1, INB_RELT+1 + XPOS_SEC_RG(JROT,JAZI,JRAD,:) = XPOS_CYL_RG(JROT,:) & + + MATMUL(XMAT_RG_RC(JROT,:,:),GET_VEC_CART(XPOS_SEC_RC(JROT,JAZI,JRAD,:))) + END DO + + DO JRAD=1, INB_RELT +!* A.1 Positioning sections centers (application points) in RG + XPOS_ELT_RG(JROT,JAZI,JRAD,:) = XPOS_CYL_RG(JROT,:) & + + MATMUL(XMAT_RG_RC(JROT,:,:),GET_VEC_CART(XPOS_ELT_RC(JROT,JAZI,JRAD,:))) +!* A.2 Update orientation +! A.2.1 Orientation in elements transition ref. frame + CALL GET_ORI_MAT_X(XANGINI_ELT_RC(JROT,JAZI,JRAD,1) & + + XRVEL_RA_RC(JROT,JAZI,JRAD,1)*ZTIME, ZORI_TMAT_X) + CALL GET_ORI_MAT_Y(XANGINI_ELT_RC(JROT,JAZI,JRAD,2) & + + XRVEL_RA_RC(JROT,JAZI,JRAD,2)*ZTIME, ZORI_TMAT_Y) + CALL GET_ORI_MAT_Z(XANGINI_ELT_RC(JROT,JAZI,JRAD,3) & + + XRVEL_RA_RC(JROT,JAZI,JRAD,3)*ZTIME, ZORI_TMAT_Z) +! Orientation matrix + XMAT_RC_RTA(JROT,JAZI,JRAD,:,:) = MATMUL(ZORI_TMAT_X, MATMUL(ZORI_TMAT_Y, ZORI_TMAT_Z)) + XMAT_RG_RTA(JROT,JAZI,JRAD,:,:) = MATMUL(XMAT_RG_RC(JROT,:,:), XMAT_RC_RTA(JROT,JAZI,JRAD,:,:)) + XMAT_RTA_RG(JROT,JAZI,JRAD,:,:) = TRANSPOSE(XMAT_RG_RTA(JROT,JAZI,JRAD,:,:)) +! +! A.2.2 Orientation in annular elements ref. frame + + CALL GET_ORI_MAT_X(XANGINI_ELT_RA(JROT,JAZI,JRAD,1), ZORI_MAT_X) + CALL GET_ORI_MAT_Y(XANGINI_ELT_RA(JROT,JAZI,JRAD,2), ZORI_MAT_Y) + CALL GET_ORI_MAT_Z(XANGINI_ELT_RA(JROT,JAZI,JRAD,3), ZORI_MAT_Z) + + XMAT_RTA_RA(JROT,JAZI,JRAD,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z)) + XMAT_RG_RA(JROT,JAZI,JRAD,:,:) = MATMUL(XMAT_RG_RTA(JROT,JAZI,JRAD,:,:), XMAT_RTA_RA(JROT,JAZI,JRAD,:,:)) + XMAT_RA_RG(JROT,JAZI,JRAD,:,:) = TRANSPOSE(XMAT_RG_RA(JROT,JAZI,JRAD,:,:)) +! + +!* A.3 Update structural velocities +! Rotation of elements in RG + XRVEL_RA_RG(JROT,JAZI,JRAD,:) = XRVEL_RC_RG(JROT,:) & + + MATMUL(XMAT_RG_RA(JROT,JAZI,JRAD,:,:),& + XRVEL_RA_RC(JROT,JAZI,JRAD,:)) +! Translation of elements in RG + ZDIST_CYL_ELT_RG(:) = XPOS_ELT_RG(JROT,JAZI,JRAD,:) - XPOS_CYL_RG(JROT,:) + XTVEL_ELT_RG(JROT,JAZI,JRAD,:) = XTVEL_CYL_RG(JROT,:) & + + CROSS(XRVEL_RA_RG(JROT,JAZI,JRAD,:),ZDIST_CYL_ELT_RG(:)) + XTVEL_ELT_RA(JROT,JAZI,JRAD,:) = MATMUL(XMAT_RA_RG(JROT,JAZI,JRAD,:,:),& + XTVEL_ELT_RG(JROT,JAZI,JRAD,:)) +! +!* A.4 Print in tecplot file + IF (LTECOUTPTS) THEN + CALL PRINT_TECOUT(ITECOUT, XPOS_ELT_RG(JROT,JAZI,JRAD,:)) + END IF +! + IF (LCSVOUTFRM) THEN + CALL PRINT_3DFRM(ICSVOUT,'ELT',XMAT_RG_RA(JROT,JAZI,JRAD,:,:),XPOS_ELT_RG(JROT,JAZI,JRAD,:)) + END IF + END DO ! Radial loop + END DO ! Azimutal loop +END DO ! Rotor loop +! +! Closing tec file +IF (LTECOUTPTS) THEN + CLOSE(ITECOUT) +END IF +IF (LCSVOUTFRM) THEN + CLOSE(ICSVOUT) +END IF +! +END SUBROUTINE EOL_KINE_ADR diff --git a/src/MNH/ini_eol_adr.f90 b/src/MNH/ini_eol_adr.f90 new file mode 100644 index 0000000000000000000000000000000000000000..1fee520cdac65855e2574a6943e6d981d6d10e8f --- /dev/null +++ b/src/MNH/ini_eol_adr.f90 @@ -0,0 +1,477 @@ +!MNH_LIC Copyright 2018-2021 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!----------------------------------------------------------------- +! ########################## + MODULE MODI_INI_EOL_ADR +! ########################## +! +INTERFACE +! +SUBROUTINE INI_EOL_ADR(PDXX,PDYY) +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY ! mesh size +! +END SUBROUTINE INI_EOL_ADR +! +END INTERFACE +! +END MODULE MODI_INI_EOL_ADR +! +! ############################################################ + SUBROUTINE INI_EOL_ADR(PDXX,PDYY) +! ############################################################ +! +!!**** *INI_EOL_ADR* - routine to initialize the Actuator Line Model +!! to simulate wind turbines +!! +!! PURPOSE +!! ------- +!! Routine to initialized the ADR (wind turbine model) variables +!! +!!** METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! *MODD_EOL_SHARED_IO: +!! *Namelist NAM_EOL_ADNR (INPUT) : +!! CHARACTER(LEN=100) :: CFARM_CSVDATA ! File to read, with farm data +!! CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! File to read, turbine data +!! CHARACTER(LEN=100) :: CBLADE_CSVDATA ! Blade file to read +!! CHARACTER(LEN=100) :: CAIRFOIL_CSVDATA ! Airfoil file to read +!! *for ouputs : +!! REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N] +!! REAL, DIMENSION(:), ALLOCATABLE :: XTORQT ! Torque [Nm] +!! REAL, DIMENSION(:), ALLOCATABLE :: XPOWT ! Power [W] +!! REAL, DIMENSION(:), ALLOCATABLE :: XTHRU_SUM ! Sum of thrust (N) +!! REAL, DIMENSION(:), ALLOCATABLE :: XTORQ_SUM ! Sum of torque (Nm) +!! REAL, DIMENSION(:), ALLOCATABLE :: XPOW_SUM ! Sum of power (W) +!! +!! INTEGER :: NNB_RADELT ! Number of radial elements +!! INTEGER :: NNB_AZIELT ! Number of azimutal elements +!! +!! *MODD_EOL_ADR (OUTPUT): +!! TYPE(FARM) :: TFARM +!! TYPE(TURBINE) :: TTURBINE +!! TYPE(BLADE) :: TBLADE +!! TYPE(AIRFOIL), DIMENSION(:), ALLOCATABLE :: TAIRFOIL +!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XELT_RAD ! Blade elements radius [m] +!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAOA_GLB ! Angle of attack of an element [rad] +!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFLIFT_GLB ! Lift force, parallel to Urel [N] +!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFDRAG_GLB ! Drag force, perpendicular to Urel [N] +!! REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_REA_GLB ! Aerodyn. force (lift+drag) in RE [N] +!! REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RG_GLB ! Aerodyn. force (lift+drag) in RG [N] +!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAOA_SUM ! Sum of angle of attack [rad] +!! REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_REA_SUM ! Sum of aerodyn. force (lift+drag) in RE [N] +! +!! *MODD_EOL_KINEMATICS +!! Positions +!! Orientations +!! Velocities +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! PA. Joulin * Meteo France & IFPEN * +!! +!! MODIFICATIONS +!! ------------- +!! Original 31/05/18 +!! Modification 10/11/20 (PA. Joulin) Updated for a main version +!! +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +!* 0.1 Modules +! +USE MODD_EOL_ADR +USE MODD_EOL_SHARED_IO, ONLY: CFARM_CSVDATA +USE MODD_EOL_SHARED_IO, ONLY: CTURBINE_CSVDATA +USE MODD_EOL_SHARED_IO, ONLY: CBLADE_CSVDATA +USE MODD_EOL_SHARED_IO, ONLY: CAIRFOIL_CSVDATA +USE MODD_EOL_SHARED_IO, ONLY: XTHRUT, XTORQT, XPOWT +USE MODD_EOL_SHARED_IO, ONLY: XTHRU_SUM, XTORQ_SUM, XPOW_SUM +USE MODD_EOL_SHARED_IO, ONLY: XELT_RAD +USE MODD_EOL_SHARED_IO, ONLY: XAOA_GLB, XFLIFT_GLB, XFDRAG_GLB, XFAERO_RG_GLB +USE MODD_EOL_SHARED_IO, ONLY: XAOA_SUM +USE MODI_EOL_READER, ONLY: READ_CSVDATA_FARM_ADR +USE MODI_EOL_READER, ONLY: READ_CSVDATA_TURBINE_ADR +USE MODI_EOL_READER, ONLY: READ_CSVDATA_BLADE_ADR +USE MODI_EOL_READER, ONLY: READ_CSVDATA_AIRFOIL_ADR +USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_FARM_ADR +USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_TURBINE_ADR +USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_BLADE_ADR +USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_AIRFOIL_ADR +USE MODD_EOL_KINE_ADR +USE MODI_EOL_MATHS +! To print in output listing +USE MODD_LUNIT_n, ONLY: TLUOUT +! Constant +USE MODD_CST, ONLY: XPI +! To know the grid +USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZS +USE MODE_ll, ONLY: GET_INDICE_ll +USE MODD_PARAMETERS, ONLY: JPVEXT +! MPI stuffs +USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD +USE MODD_PRECISION, ONLY: MNHREAL_MPI +USE MODD_MPIF, ONLY: MPI_SUM +! +!* 0.2 Variables +! +IMPLICIT NONE +! Interface +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY ! mesh size +! +! Some loop controlers +! .. for wind turbines +INTEGER :: JROT ! Rotor index +INTEGER :: JTELT ! Tower element index +INTEGER :: JNELT ! Nacelle element index +INTEGER :: JRAD ! Radius index +INTEGER :: JAZI ! Azimut index +! .. for domain +INTEGER :: IIB,IJB,IKB ! Begin of a CPU domain +INTEGER :: IIE,IJE ! End of a CPU domain +INTEGER :: JI,JJ ! Domain index +! Some variables to be coder-friendly +INTEGER :: INB_WT, INB_B ! Total numbers of wind turbines, blade +INTEGER :: INB_RELT, INB_AELT ! Total numbers of radial elt and azimut elt +INTEGER :: INB_TELT, INB_NELT ! Total numbers of tower elt and nacelle elt +REAL :: ZRAD ! Radius along the blade +! Tower base folowing the terrain +REAL,DIMENSION(:,:),ALLOCATABLE :: ZPOSINI_TOWO_RG ! Initial tower origin position +INTEGER :: IINFO ! code info return +! +!------------------------------------------------------------------- +! +!* 1. READING AND ALLOCATING DATA +! --------------------------- +! Reading in csv files +! Allocation of TFARM, TTURBINE, TBLADE and TAIRFOILS inside the function +! +!* 1.1 Wind farm data +! +CALL READ_CSVDATA_FARM_ADR(TRIM(CFARM_CSVDATA),TFARM) +! +!* 1.2 Wind turbine data +! +CALL READ_CSVDATA_TURBINE_ADR(TRIM(CTURBINE_CSVDATA),TTURBINE) +! +!* 1.3 Blade data +! +CALL READ_CSVDATA_BLADE_ADR(TRIM(CBLADE_CSVDATA),TTURBINE,TBLADE) +! +!* 1.4 Airfoil data +! +CALL READ_CSVDATA_AIRFOIL_ADR(TRIM(CAIRFOIL_CSVDATA),TBLADE,TAIRFOIL) +! +! +!------------------------------------------------------------------- +! +!* 2. PRINTING DATA +! ------------- +! +!* 2.1 Wind farm data +! +CALL PRINT_DATA_FARM_ADR(TLUOUT%NLU,TFARM) +! +!* 2.2 Wind turbine data +! +CALL PRINT_DATA_TURBINE_ADR(TLUOUT%NLU,TTURBINE) +! +!* 2.3 Blade data +! +CALL PRINT_DATA_BLADE_ADR(TLUOUT%NLU,TBLADE) +! +!* 2.4 Airfoil data +! +CALL PRINT_DATA_AIRFOIL_ADR(TLUOUT%NLU,TAIRFOIL) +! +! +!------------------------------------------------------------------- +! +!* 3. ALLOCATING VARIABLES +! -------------------- +! +!* 3.0 Preliminaries +INB_WT = TFARM%NNB_TURBINES +INB_B = TTURBINE%NNB_BLADES +INB_AELT = TBLADE%NNB_AZIELT +INB_RELT = TBLADE%NNB_RADELT +! Hard coded variables, but they will be usefull in next updates +INB_TELT = 2 +INB_NELT = 2 +! +XDELTA_AZI = 2d0*XPI/INB_AELT ! Rotor disc azimutal devision +XDELTA_RAD = (TTURBINE%XR_MAX - TTURBINE%XR_MIN)/INB_RELT ! Rotor disc radialal devision + + +!* 3.1 MODD_EOL_ADR variables +! at t +ALLOCATE(XELT_RAD (INB_WT,INB_AELT,INB_RELT )) +ALLOCATE(XELT_AZI (INB_WT,INB_AELT,INB_RELT )) +ALLOCATE(XAOA_ELT (INB_WT,INB_AELT,INB_RELT )) +ALLOCATE(XAOA_GLB (INB_WT,INB_AELT,INB_RELT )) +ALLOCATE(XAOA_BLA (INB_WT,INB_RELT )) +ALLOCATE(XFDRAG_ELT (INB_WT,INB_AELT,INB_RELT )) +ALLOCATE(XFDRAG_GLB (INB_WT,INB_AELT,INB_RELT )) +ALLOCATE(XFLIFT_ELT (INB_WT,INB_AELT,INB_RELT )) +ALLOCATE(XFLIFT_GLB (INB_WT,INB_AELT,INB_RELT )) +ALLOCATE(XFAERO_RA_ELT (INB_WT,INB_AELT,INB_RELT,3)) +ALLOCATE(XFAERO_RA_GLB (INB_WT,INB_AELT,INB_RELT,3)) +ALLOCATE(XFAERO_RA_BLA (INB_WT,INB_RELT,3)) +ALLOCATE(XFAERO_RG_ELT (INB_WT,INB_AELT,INB_RELT,3)) +ALLOCATE(XFAERO_RG_GLB (INB_WT,INB_AELT,INB_RELT,3)) +ALLOCATE(XTHRUT (INB_WT )) +ALLOCATE(XTORQT (INB_WT )) +ALLOCATE(XPOWT (INB_WT )) +! for mean values +ALLOCATE(XAOA_SUM_ADR (INB_WT,INB_RELT )) +ALLOCATE(XFAERO_RA_SUM (INB_WT,INB_RELT,3)) +ALLOCATE(XTHRU_SUM (INB_WT)) +ALLOCATE(XTORQ_SUM (INB_WT)) +ALLOCATE(XPOW_SUM (INB_WT)) +! +!* 3.2 MODD_EOL_KINE_ADR variables +! +! - ORIENTATION MATRIX - +ALLOCATE(XMAT_RG_RT(INB_WT,3,3)) +ALLOCATE(XMAT_RG_RN(INB_WT,3,3)) +ALLOCATE(XMAT_RT_RN(INB_WT,3,3)) +ALLOCATE(XMAT_RG_RH(INB_WT,3,3)) +ALLOCATE(XMAT_RH_RG(INB_WT,3,3)) +ALLOCATE(XMAT_RN_RH(INB_WT,3,3)) +ALLOCATE(XMAT_RG_RC(INB_WT,3,3)) +ALLOCATE(XMAT_RH_RC(INB_WT,3,3)) +ALLOCATE(XMAT_RG_RTA(INB_WT,INB_AELT,INB_RELT,3,3)) +ALLOCATE(XMAT_RTA_RG(INB_WT,INB_AELT,INB_RELT,3,3)) +ALLOCATE(XMAT_RC_RTA(INB_WT,INB_AELT,INB_RELT,3,3)) +ALLOCATE(XMAT_RG_RA(INB_WT,INB_AELT,INB_RELT,3,3)) +ALLOCATE(XMAT_RA_RG(INB_WT,INB_AELT,INB_RELT,3,3)) +ALLOCATE(XMAT_RTA_RA(INB_WT,INB_AELT,INB_RELT,3,3)) +! +! - POSITIONS & ORIENTATIONS - +! Tower +ALLOCATE(ZPOSINI_TOWO_RG (INB_WT,3) ) ! Initial tower origin pos. in RG +ALLOCATE(XPOSINI_TOWO_RG (INB_WT,3) ) ! Initial tower origin pos. in RG +ALLOCATE(XPOS_TOWO_RG (INB_WT,3) ) ! Current tower origin pos. in RG +ALLOCATE(XPOS_TELT_RG (INB_WT,INB_TELT,3) ) ! Current tower element pos. in RG +ALLOCATE(XPOS_TELT_RT (INB_WT,INB_TELT,3) ) ! Current tower element pos. in RT +ALLOCATE(XANGINI_TOW_RG (INB_WT,3) ) ! Initial tower ori. in RG +! Nacelle +ALLOCATE(XPOSINI_NACO_RT (INB_WT,3) ) ! Initial nacelle origin pos. in RT +ALLOCATE(XPOS_NACO_RG (INB_WT,3) ) ! Current nacelle origin pos. in RG +ALLOCATE(XPOS_NELT_RG (INB_WT,INB_NELT,3) ) ! Current nacelle element pos. in RG +ALLOCATE(XPOS_NELT_RN (INB_WT,INB_NELT,3) ) ! Current nacelle element pos. in RN +ALLOCATE(XANGINI_NAC_RT (INB_WT,3) ) ! Initial nacelle ori. in RT +! Hub +ALLOCATE(XPOSINI_HUB_RN (INB_WT,3) ) ! Initial hub pos. in RN +ALLOCATE(XPOS_HUB_RG (INB_WT,3) ) ! Current hub pos. in RG +ALLOCATE(XANGINI_HUB_RN (INB_WT,3) ) ! Initial hub ori. in RN +! Cylindrical +ALLOCATE(XPOSINI_CYL_RH (INB_WT,3) ) ! Initial cylindrical center pos. in RH +ALLOCATE(XPOS_CYL_RG (INB_WT,3) ) ! Current cylindrical center pos. in RG +ALLOCATE(XANGINI_CYL_RH (INB_WT,3) ) ! Initial cylindrical center ori. RH +! Annular Element +ALLOCATE(XPOS_ELT_RC (INB_WT,INB_AELT,INB_RELT,3) ) ! Element pos. in RC +ALLOCATE(XPOS_ELT_RG (INB_WT,INB_AELT,INB_RELT,3) ) ! Element pos. in RG +ALLOCATE(XPOS_SEC_RC (INB_WT,INB_AELT,INB_RELT+1,3)) ! Section pos. in RC +ALLOCATE(XPOS_SEC_RG (INB_WT,INB_AELT,INB_RELT+1,3)) ! Section pos. in RG +ALLOCATE(XANGINI_ELT_RC (INB_WT,INB_AELT,INB_RELT,3)) ! Initial element ori. in RC +ALLOCATE(XANGINI_ELT_RA (INB_WT,INB_AELT,INB_RELT,3)) ! Initial element ori. in transition ref. +ALLOCATE(XTWIST_ELT (INB_WT,INB_RELT) ) ! Element twist in RC +ALLOCATE(XCHORD_ELT (INB_WT,INB_RELT) ) ! Element chord lenght +ALLOCATE(XSURF_ELT (INB_WT,INB_RELT) ) ! Element lift surface +ALLOCATE(XSURF_APP_ELT (INB_WT,INB_RELT) ) ! Element lift surface +! +! - STRUCTURAL VELOCITIES - +! Tower +ALLOCATE(XTVEL_TOWO_RG (INB_WT,3) ) ! Tower base trans. vel. in RG +ALLOCATE(XTVEL_TELT_RG (INB_WT,INB_TELT,3) ) ! Tower element trans. vel. in RG +ALLOCATE(XRVEL_RT_RG (INB_WT,3) ) ! RT/RG rot. vel. +! Nacelle +ALLOCATE(XTVEL_NACO_RT (INB_WT,3) ) ! Nacelle base trans. vel. in RT +ALLOCATE(XTVEL_NELT_RG (INB_WT,INB_NELT,3) ) ! Nacelle element trans. vel. in RG +ALLOCATE(XRVEL_RN_RT (INB_WT,3) ) ! RN/RT rot. vel. +ALLOCATE(XRVEL_RN_RG (INB_WT,3) ) ! RN/RG rot. vel. +! Hub +ALLOCATE(XTVEL_HUB_RN (INB_WT,3) ) ! Hub base trans. vel. in RN +ALLOCATE(XTVEL_HUB_RG (INB_WT,3) ) ! Hub base trans. vel. in RG +ALLOCATE(XRVEL_RH_RN (INB_WT,3) ) ! RH/RT rot. vel. +ALLOCATE(XRVEL_RH_RG (INB_WT,3) ) ! RH/RG rot. vel. +! Cylindrical +ALLOCATE(XTVEL_CYL_RH (INB_WT,3) ) ! Cylindrical base trans. vel. in RH +ALLOCATE(XTVEL_CYL_RG (INB_WT,3) ) ! Cylindrical base trans. vel. in RG +ALLOCATE(XRVEL_RC_RH (INB_WT,3) ) ! RC/RH rot. vel. +ALLOCATE(XRVEL_RC_RG (INB_WT,3) ) ! RC/RG rot. vel. +! Annular Elements +ALLOCATE(XTVEL_ELT_RC (INB_WT,INB_AELT,INB_RELT,3) ) ! Element base trans. vel. in RC +ALLOCATE(XTVEL_ELT_RG (INB_WT,INB_AELT,INB_RELT,3) ) ! Element base trans. vel. in RG +ALLOCATE(XTVEL_ELT_RA (INB_WT,INB_AELT,INB_RELT,3) ) ! Element base trans. vel. in RA +ALLOCATE(XRVEL_RA_RC (INB_WT,INB_AELT,INB_RELT,3) ) ! RA/RB rot. vel. +ALLOCATE(XRVEL_RA_RG (INB_WT,INB_AELT,INB_RELT,3) ) ! RA/RG rot. vel. +! +!------------------------------------------------------------------- +! +!* 4. FIRST BUILDING OF WIND TURBINES +! ------------------------------- +! +!* 4.1 Preliminaries + + +CALL GET_INDICE_ll(IIB,IJB,IIE,IJE) ! Get begin and end domain index (CPU) +IKB=1+JPVEXT ! Vertical begin index +! +XPOS_REF(:) = 0d0 ! Global Origin +! +! +!* 4.2 Tower +DO JROT=1, INB_WT +! Velocities (Rotor, (XYZ)) + XTVEL_TOWO_RG(JROT,:) = 0d0 + XRVEL_RT_RG(JROT,:) = 0d0 +! Positions + ! Init + ZPOSINI_TOWO_RG(JROT,1) = 0 + ZPOSINI_TOWO_RG(JROT,2) = 0 + ZPOSINI_TOWO_RG(JROT,3) = 0 + ! Finding the elevation at the WT position + DO JJ=IJB,IJE + DO JI=IIB,IIE + IF (TFARM%XPOS_X(JROT) >= XXHAT(JI) .AND. & + TFARM%XPOS_X(JROT) < XXHAT(JI) + PDXX(JI,JJ,IKB)) THEN + IF (TFARM%XPOS_Y(JROT) >= XYHAT(JJ) .AND. & + TFARM%XPOS_Y(JROT) < XYHAT(JJ) + PDYY(JI,JJ,IKB)) THEN + ! Horizontal position + ZPOSINI_TOWO_RG(JROT,1) = TFARM%XPOS_X(JROT) ! Tower base position + ZPOSINI_TOWO_RG(JROT,2) = TFARM%XPOS_Y(JROT) + ! Finding the elevation at the WT position + ZPOSINI_TOWO_RG(JROT,3) = XZS(JI,JJ) + END IF + END IF + END DO + END DO + ! Sharing information + CALL MPI_ALLREDUCE(ZPOSINI_TOWO_RG, XPOSINI_TOWO_RG, SIZE(XPOSINI_TOWO_RG),& + MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO) + ! Tower elements + DO JTELT=1, INB_TELT + XPOS_TELT_RT(JROT,JTELT,1) = 0d0 + XPOS_TELT_RT(JROT,JTELT,2) = 0d0 + XPOS_TELT_RT(JROT,JTELT,3) = (JTELT-1)*(TTURBINE%XH_HEIGHT)/(INB_TELT-1) + END DO + ! Angles + XANGINI_TOW_RG(JROT,1) = 0d0 + XANGINI_TOW_RG(JROT,2) = 0d0 + XANGINI_TOW_RG(JROT,3) = 0d0 +! +! +!* 4.3 Nacelle +! Velocities + XTVEL_NACO_RT(JROT,:) = 0d0 + XRVEL_RN_RT(JROT,:) = 0d0 +! Positions (Rotor, (XYZ)) ! From last point of tower + ! Origin + XPOSINI_NACO_RT(JROT,:) = 0d0 ! Distance between nacelle base and tower top + ! Elements + DO JNELT=1, INB_NELT + XPOS_NELT_RN(JROT,JNELT,1) = (JNELT-1)*TTURBINE%XH_DEPORT/(INB_NELT-1) + XPOS_NELT_RN(JROT,JNELT,2) = 0d0 + XPOS_NELT_RN(JROT,JNELT,3) = 0d0 + END DO + ! Angles + XANGINI_NAC_RT(JROT,1) = 0d0 + XANGINI_NAC_RT(JROT,2) = 0d0 + XANGINI_NAC_RT(JROT,3) = XPI + TFARM%XNAC_YAW(JROT) +! +! +!* 4.4 Hub +! Velocities + XTVEL_HUB_RN(JROT,:) = 0d0 + XRVEL_RH_RN(JROT,1) = 0d0 + XRVEL_RH_RN(JROT,2) = 0d0 + XRVEL_RH_RN(JROT,3) = TFARM%XOMEGA(JROT) +! Position (Rotor, (XYZ)) ! From nacelle last point + XPOSINI_HUB_RN(JROT,1) = 0d0 + XPOSINI_HUB_RN(JROT,2) = 0d0 + XPOSINI_HUB_RN(JROT,3) = 0d0 + XANGINI_HUB_RN(JROT,1) = 0d0 + XANGINI_HUB_RN(JROT,2) = XPI/2d0 - TTURBINE%XNAC_TILT + XANGINI_HUB_RN(JROT,3) = 0d0 +! +! +!* 4.5 Cylindrical +! Velocities + XRVEL_RC_RH(JROT,:) = 0d0 + XTVEL_CYL_RH(JROT,:) = 0d0 +! Position + XPOSINI_CYL_RH(JROT,:) = 0d0 + XANGINI_CYL_RH(JROT,1) = 0d0 + XANGINI_CYL_RH(JROT,2) = 0d0 + XANGINI_CYL_RH(JROT,3) = XPI +! +! +!* 4.6 Annular Elements + + DO JAZI=1, INB_AELT + + + DO JRAD=1, INB_RELT+1 + +! - Positioning of sections (cuts) + + XPOS_SEC_RC(JROT,JAZI,JRAD,1) = TTURBINE%XR_MIN + (JRAD-1) * XDELTA_RAD + XPOS_SEC_RC(JROT,JAZI,JRAD,2) = (JAZI-1) * XDELTA_AZI + XPOS_SEC_RC(JROT,JAZI,JRAD,3) = 0d0 + END DO + + DO JRAD=1, INB_RELT +! - Positioning of centers (points of application) + XPOS_ELT_RC(JROT,JAZI,JRAD,1) = XPOS_SEC_RC(JROT,JAZI,JRAD,1) + XDELTA_RAD/2d0 + XPOS_ELT_RC(JROT,JAZI,JRAD,2) = XPOS_SEC_RC(JROT,JAZI,JRAD,2) + XDELTA_AZI/2d0 + XPOS_ELT_RC(JROT,JAZI,JRAD,3) = 0d0 + XELT_RAD(JROT,JAZI,JRAD) = XPOS_ELT_RC(JROT,JAZI,JRAD,1) + XELT_AZI(JROT,JAZI,JRAD) = XPOS_ELT_RC(JROT,JAZI,JRAD,2) + + +! - Calculating chord and twist + ZRAD = XELT_RAD(JROT,JAZI,JRAD) + XCHORD_ELT(JROT,JRAD) = INTERP_SPLCUB(ZRAD, TBLADE%XRAD, TBLADE%XCHORD) + XTWIST_ELT(JROT,JRAD) = INTERP_SPLCUB(ZRAD, TBLADE%XRAD, TBLADE%XTWIST) +! - Calculating lifting surface + XSURF_ELT(JROT,JRAD) = XCHORD_ELT(JROT,JRAD) & + * (XPOS_SEC_RC(JROT,JAZI,JRAD+1,1) & + - XPOS_SEC_RC(JROT,JAZI,JRAD,1)) + XSURF_APP_ELT(JROT,JRAD) = INB_B * XSURF_ELT(JROT,JRAD) & + * XDELTA_AZI/(2d0*XPI) + +! - Velocities + XRVEL_RA_RC(JROT,JAZI,JRAD,:) = 0d0 + XTVEL_ELT_RC(JROT,JAZI,JRAD,:) = 0d0 + +! - Orientation +! - Orientation in the transition reference frame + XANGINI_ELT_RC(JROT,JAZI,JRAD,1) = 0d0 + XANGINI_ELT_RC(JROT,JAZI,JRAD,2) = 0d0 + XANGINI_ELT_RC(JROT,JAZI,JRAD,3) = 0d0 + (JAZI-1) * XDELTA_AZI + XDELTA_AZI/2d0 +! - Orientation in the annular elements reference frame + XANGINI_ELT_RA(JROT,JAZI,JRAD,1) = 0d0 + TFARM%XBLA_PITCH(JROT) + XTWIST_ELT(JROT,JRAD) + XANGINI_ELT_RA(JROT,JAZI,JRAD,2) = 0d0 + XANGINI_ELT_RA(JROT,JAZI,JRAD,3) = 0d0 + END DO ! Loop radial + END DO ! Loop azimutal +END DO ! Loop turbine +! +END SUBROUTINE INI_EOL_ADR diff --git a/src/MNH/modd_eol_adr.f90 b/src/MNH/modd_eol_adr.f90 new file mode 100644 index 0000000000000000000000000000000000000000..7b5520dcfae0ec80c8fae41e03993242cc2b1c12 --- /dev/null +++ b/src/MNH/modd_eol_adr.f90 @@ -0,0 +1,122 @@ +!MNH_LIC Copyright 2017-2021 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!----------------------------------------------------------------- +!! +!! ##################### + MODULE MODD_EOL_ADR +!! ##################### +!! +!!*** *MODD_EOL_ALM* +!! +!! PURPOSE +!! ------- +!! It is possible to include wind turbines parameterization in Meso-NH, +!! and several models are available. One of the models is the Actuator +!! Line Method (ALM). MODD_EOL_ALM contains all the declarations for +!! the ALM model. +!! +!! +!!** AUTHOR +!! ------ +!! PA.Joulin *CNRM & IFPEN* +! +!! MODIFICATIONS +!! ------------- +!! Original 04/01/17 +!! Modification 14/10/20 (PA. Joulin) Updated for a main version +!! +!----------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ----------------- +! +IMPLICIT NONE +! +! ------ TYPES ------ +! +TYPE FARM + INTEGER :: NNB_TURBINES ! Number of wind turbines [-] + REAL, DIMENSION(:), ALLOCATABLE :: XPOS_X ! Tower base position, X coord [m] + REAL, DIMENSION(:), ALLOCATABLE :: XPOS_Y ! Tower base position, Y coord [m] + REAL, DIMENSION(:), ALLOCATABLE :: XOMEGA ! Rotor rotation speed [rad/s] + REAL, DIMENSION(:), ALLOCATABLE :: XNAC_YAW ! Nacelle yaw angle [rad] + REAL, DIMENSION(:), ALLOCATABLE :: XBLA_PITCH ! Blade picth angle [rad] +END TYPE FARM +! +TYPE TURBINE + CHARACTER(LEN=10) :: CNAME ! Nom de la turbine [-] + INTEGER :: NNB_BLADES ! Number of blades [-] + REAL :: XH_HEIGHT ! Hub height [m] + REAL :: XH_DEPORT ! Hub deport [m] + REAL :: XNAC_TILT ! Tilt of the nacelle [m] + REAL :: XR_MIN ! Minimum blade radius [m] + REAL :: XR_MAX ! Maximum blade radius [m] +END TYPE TURBINE +! +TYPE BLADE + INTEGER :: NNB_RADELT ! Number of radial elements + INTEGER :: NNB_AZIELT ! Number of azimutal elements + INTEGER :: NNB_BLADAT ! Number of blade data + REAL, DIMENSION(:), ALLOCATABLE :: XRAD ! Data node radius [m] + REAL, DIMENSION(:), ALLOCATABLE :: XCHORD ! Element chord [m] + REAL, DIMENSION(:), ALLOCATABLE :: XTWIST ! Element twist [rad] + CHARACTER(LEN=20), DIMENSION(:), ALLOCATABLE :: CAIRFOIL ! Arifoil name [-] +END TYPE BLADE +! +TYPE AIRFOIL + CHARACTER(LEN=15) :: CNAME ! Airfoil name [-] + REAL, DIMENSION(:), ALLOCATABLE :: XAA ! Attack Angle [rad] + REAL, DIMENSION(:), ALLOCATABLE :: XRE ! Reynolds Number [-] + REAL, DIMENSION(:), ALLOCATABLE :: XCL ! Lift coefficient [-] + REAL, DIMENSION(:), ALLOCATABLE :: XCD ! Drag coefficient [-] + REAL, DIMENSION(:), ALLOCATABLE :: XCM ! Moment coefficient [-] +END TYPE AIRFOIL +! +! ------ VARIABLES ------ +! --- Farm data --- +TYPE(FARM) :: TFARM +TYPE(TURBINE) :: TTURBINE +TYPE(BLADE) :: TBLADE +TYPE(AIRFOIL), DIMENSION(:), ALLOCATABLE :: TAIRFOIL +! +! --- Global variables (Code & CPU) --- +REAL :: XDELTA_AZI ! Delta azimut [rad] +REAL :: XDELTA_RAD ! Delta radius [rad] +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XELT_AZI ! Elements azimut [rad] +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RA_GLB ! Aerodyn. force (lift+drag) in RA [N] +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RA_ELT ! Aerodyn. force (lift+drag) in RA [N] +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFAERO_RA_BLA ! Aerodyn. force (lift+drag) in RA [N] +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAOA_ELT +REAL, DIMENSION(:,:), ALLOCATABLE :: XAOA_BLA +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFDRAG_ELT +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFLIFT_ELT +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RG_ELT ! Aerodyn. force (lift+drag) in RA [N] + +! Mean values +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFAERO_RA_SUM ! Sum of aerodyn. force (lift+drag) in RA [N] +REAL, DIMENSION(:,:), ALLOCATABLE :: XAOA_SUM_ADR ! Sum of aerodyn. force (lift+drag) in RA [N] +! +! Implicit from MODD_EOL_SHARED_IO : +!REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N] +!REAL, DIMENSION(:), ALLOCATABLE :: XTORQT ! Torque [Nm] +!REAL, DIMENSION(:), ALLOCATABLE :: XPOWT ! Power [W] +!REAL, DIMENSION(:), ALLOCATABLE :: XTHRU_SUM ! Sum of thrust (N) +!REAL, DIMENSION(:), ALLOCATABLE :: XTORQ_SUM ! Sum of torque (Nm) +!REAL, DIMENSION(:), ALLOCATABLE :: XPOW_SUM ! Sum of power (W) +! +! +! --- Namelist NAM_EOL_ADR --- +! Implicit from MODD_EOL_SHARED_IO : +!CHARACTER(LEN=100) :: CFARM_CSVDATA ! Farm file to read +!CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! Turbine file to read +!CHARACTER(LEN=100) :: CBLADE_CSVDATA ! Blade file to read +!CHARACTER(LEN=100) :: CAIRFOIL_CSVDATA ! Airfoil file to read +!CHARACTER(LEN=3) :: CINTERP ! Interpolation method for wind speed +! +INTEGER :: NNB_RADELT ! Number of radial elements +INTEGER :: NNB_AZIELT ! Number of azimutal elements +! + +END MODULE MODD_EOL_ADR diff --git a/src/MNH/modd_eol_kine_adr.f90 b/src/MNH/modd_eol_kine_adr.f90 new file mode 100644 index 0000000000000000000000000000000000000000..df90e6d3e659b30803b153be14c7d9b07608ffcc --- /dev/null +++ b/src/MNH/modd_eol_kine_adr.f90 @@ -0,0 +1,112 @@ +!MNH_LIC Copyright 2018-2021 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!----------------------------------------------------------------- +!! +!! ##################### + MODULE MODD_EOL_KINE_ADR +!! ##################### +!! +!!*** *MODD_EOL_KINE_ALM* +!! +!! PURPOSE +!! ------- +! Declaration to take into account wind turbine motion +! +!! +!!** AUTHOR +!! ------ +!! PA.Joulin *CNRM & IFPEN* +! +!! MODIFICATIONS +!! ------------- +!! Original 04/18 +! P. Wautelet 19/07/2021: replace double precision by real to allow MNH_REAL=4 compilation +!! +!----------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ----------------- +! +! +! - Matrix to move from one fram to an other - +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XMAT_RG_RT ! RG = Repere GLOBAL, RT = Repere TOWER +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XMAT_RG_RN, XMAT_RT_RN ! RN = Repere NACELLE +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XMAT_RG_RH, XMAT_RH_RG, XMAT_RN_RH ! RH = Repere HUB +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XMAT_RG_RC, XMAT_RH_RC ! RC = Repere CYLINDRICAL +REAL, DIMENSION(:,:,:,:,:), ALLOCATABLE :: XMAT_RG_RTA, XMAT_RTA_RG, XMAT_RC_RTA ! RTA = Repere TRANSITION ANNULAR +REAL, DIMENSION(:,:,:,:,:), ALLOCATABLE :: XMAT_RG_RA, XMAT_RA_RG, XMAT_RTA_RA ! RE = Repere ANNULAR + +! - POSITIONS & ORIENTATIONS - +REAL, DIMENSION(3) :: XPOS_REF ! Reference position + +! Tower +REAL, DIMENSION(:,:), ALLOCATABLE :: XPOSINI_TOWO_RG ! Initial tower origin position, in global reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XPOS_TOWO_RG ! Current tower origin real position, in global reference frame +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XPOS_TELT_RG ! Current tower element position, in global reference frame +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XPOS_TELT_RT ! Current tower element position, in tower frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XANGINI_TOW_RG ! Initial tower orientation in global ref frame + +! Nacelle +REAL, DIMENSION(:,:), ALLOCATABLE :: XPOSINI_NACO_RT ! Initial nacelle position, in tower reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XPOS_NACO_RG ! Initial nacelle position, in global reference frame +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XPOS_NELT_RG ! Initial nacelle position, in global reference frame +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XPOS_NELT_RN ! Initial nacelle position, in nacelle reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XANGINI_NAC_RT ! Initial nacelle orientation, in tower reference frame + +! Hub +REAL, DIMENSION(:,:), ALLOCATABLE :: XPOSINI_HUB_RN ! Initial hub position, in nacelle reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XPOS_HUB_RG ! Current hub position, in global reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XANGINI_HUB_RN ! Initial hub orientation, in nacelle reference frame + +! Cylindrical +REAL, DIMENSION(:,:), ALLOCATABLE :: XPOSINI_CYL_RH ! Initial blade root position, in hub reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XPOS_CYL_RG ! Current blade root position, in global reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XANGINI_CYL_RH ! Initial blade orientation, in hub reference frame + +! Annular Element +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XPOS_ELT_RC ! Element position, in cylindrical reference frame +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XPOS_ELT_RG ! Element position, in global reference frame +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XPOS_SEC_RC ! Section position, in cylindrical reference frame +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XPOS_SEC_RG ! Section position, in global reference frame +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XANGINI_ELT_RC ! Initial element orientation in cyl reference frame +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XANGINI_ELT_RA ! Initial element orientation in annular elt transition ref frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XTWIST_ELT ! Element twist, interpolated from data +REAL, DIMENSION(:,:), ALLOCATABLE :: XCHORD_ELT ! Element chord lenght, interpolated from data +REAL, DIMENSION(:,:), ALLOCATABLE :: XSURF_ELT ! Element lift surface +REAL, DIMENSION(:,:), ALLOCATABLE :: XSURF_APP_ELT ! Element lift surface +! +! +! - STRUCTURAL VELOCITIES - +! Tower +REAL, DIMENSION(:,:), ALLOCATABLE :: XTVEL_TOWO_RG ! Tower base translation velocity, in global reference frame +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XTVEL_TELT_RG ! Tower element velocity, in global reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RT_RG ! RT/RG rotational velocity + +! Nacelle +REAL, DIMENSION(:,:), ALLOCATABLE :: XTVEL_NACO_RT ! Nacelle base translation velocity, in tower reference frame +REAL, DIMENSION(:,:,:), ALLOCATABLE :: XTVEL_NELT_RG ! Nacelle element translation velocity, in global reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RN_RT ! RN/RT rotational velocity +REAL, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RN_RG ! RN/RG rotational velocity + +! Hub +REAL, DIMENSION(:,:), ALLOCATABLE :: XTVEL_HUB_RN ! Hub base translation velocity, in global reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XTVEL_HUB_RG ! Hub base translation velocity, in global reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RH_RN ! RH/RN rotational velocity +REAL, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RH_RG ! RH/RG rotational velocity + +! Cylindrical +REAL, DIMENSION(:,:), ALLOCATABLE :: XTVEL_CYL_RH ! Cyllindrical base translation velocity, in hub reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XTVEL_CYL_RG ! Cylindrical base translation velocity, in global reference frame +REAL, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RC_RH ! RC/RH rotational velocity +REAL, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RC_RG ! RC/RG rotational velocity + +! Annular Elements +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XTVEL_ELT_RC ! Element base translation velocity, in cylindrical reference frame +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XTVEL_ELT_RG ! Element base translation velocity, in global reference frame +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XTVEL_ELT_RA ! Element base translation velocity, in annular element reference frame +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XRVEL_RA_RC ! RA/RC rotational velocity +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XRVEL_RA_RG ! RA/RG rotational velocity + +END MODULE MODD_EOL_KINE_ADR diff --git a/src/MNH/modn_eol_adr.f90 b/src/MNH/modn_eol_adr.f90 new file mode 100644 index 0000000000000000000000000000000000000000..e473a8fae6be6a3f03e6bf1caa28d002f6289975 --- /dev/null +++ b/src/MNH/modn_eol_adr.f90 @@ -0,0 +1,46 @@ +!MNH_LIC Copyright 2017-2021 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!----------------------------------------------------------------- +!! +!! ##################### + MODULE MODN_EOL_ADR +!! ##################### +!! +!!*** *MODN_EOL_ALM* +!! +!! PURPOSE +!! ------- +!! NAM_EOL activate the parameterization of wind turbines, and several +!! models are available. One of the models is the Actuator Line +!! Method (ALM). +!! The aim of NAM_EOL_ALM is to specify ALM parameters. +!! +!!** AUTHOR +!! ------ +!! PA. Joulin *CNRM & IFPEN* +! +!! MODIFICATIONS +!! ------------- +!! Original 24/01/17 +!! Modification 14/10/20 (PA. Joulin) Updated for a main version +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +USE MODD_EOL_ADR +USE MODD_EOL_SHARED_IO +!! +!----------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ----------------- +IMPLICIT NONE +SAVE +NAMELIST /NAM_EOL_ADR/ & + CFARM_CSVDATA, CTURBINE_CSVDATA, CBLADE_CSVDATA, CAIRFOIL_CSVDATA, & + NNB_AZIELT, NNB_RADELT, & + CINTERP, LTIPLOSSG, & + LTECOUTPTS,LCSVOUTFRM +! +END MODULE MODN_EOL_ADR