Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
!MNH_LIC Copyright 1994-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_ALM
!
INTERFACE !----------------------------------------
!
SUBROUTINE EOL_KINE_ALM(KTCOUNT,KTSUBCOUNT,PTSUBSTEP,PTSTEP)
!
INTEGER, INTENT(IN) :: KTCOUNT ! iteration count
INTEGER, INTENT(IN) :: KTSUBCOUNT ! sub iteration count
REAL, INTENT(IN) :: PTSUBSTEP ! sub timestep
REAL, INTENT(IN) :: PTSTEP ! timestep
!
END SUBROUTINE EOL_KINE_ALM
!
END INTERFACE !------------------------------------
!
END MODULE MODI_EOL_KINE_ALM
!#######################
!
!
!
!
!###################################################################
SUBROUTINE EOL_KINE_ALM(KTCOUNT,KTSUBCOUNT,PTSUBSTEP,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
!!
!!---------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
!* 0.1 Modules
USE MODD_EOL_KINE_ALM
USE MODD_EOL_ALM
USE MODI_EOL_MATHS
USE MODD_TIME_n, ONLY : TDTCUR
USE MODD_CST, ONLY : XPI
!
IMPLICIT NONE
!
!* 0.2 Declarations of dummy arguments :
!
INTEGER, INTENT(IN) :: KTCOUNT ! iteration count
INTEGER, INTENT(IN) :: KTSUBCOUNT ! sub iteration count
REAL, INTENT(IN) :: PTSUBSTEP ! sub timestep
REAL, INTENT(IN) :: PTSTEP ! timestep
!
!* 0.3 Local variables
DOUBLE PRECISION, DIMENSION(3,3) :: ZORI_MAT_X, ZORI_MAT_Y, ZORI_MAT_Z
DOUBLE PRECISION, DIMENSION(3) :: ZADD_TO_POS
!
DOUBLE PRECISION, DIMENSION(3) :: ZDIST_TOWO_TELT_RG ! Distance between tower elmt and tower base
DOUBLE PRECISION, DIMENSION(3) :: ZDIST_TOWO_NELT_RG ! Distance between nacelle and base of tower
DOUBLE PRECISION, DIMENSION(3) :: ZDIST_NAC_HUB_RG ! Distance between hub and base of nacelle
DOUBLE PRECISION, DIMENSION(3) :: ZDIST_HUB_BLA_RG ! Distance between blade and base of hub
DOUBLE PRECISION, DIMENSION(3) :: ZDIST_BLA_ELT_RG ! Distance between blade and elements
!
DOUBLE PRECISION, DIMENSION(3) :: ZPOS_ELTLE_RE ! Leading Edge (LE) position, in RE
DOUBLE PRECISION, DIMENSION(3) :: ZPOS_ELTLE_RG ! Leading Edge (LE) position, in RG
DOUBLE PRECISION, DIMENSION(3) :: ZPOS_ELTTE_RE ! Trailing Edge (TE) position, in RE
DOUBLE PRECISION, DIMENSION(3) :: ZPOS_ELTTE_RG ! Trailing Edge (TE) position, in RG
!
REAL :: ZTIME ! TIME
INTEGER :: JROT, JBLA, JTELT, JNELT, JBELT ! Loop control
INTEGER :: INB_WT, INB_B, INB_BELT ! Total numbers
INTEGER :: INB_TELT, INB_NELT ! Total numbers
!
!
!---------------------------------------------------------------
!
!* 1. PRELIMINARIES
! -------------
!
!* 1.1 Some usefull integers
INB_WT = TFARM%NNB_TURBINES
INB_B = TTURBINE%NNB_BLADES
INB_TELT = 2
INB_NELT = 2
INB_BELT = TBLADE%NNB_BLAELT
!
!* 1.2 Sub-time computation
ZTIME = TDTCUR%xtime+(KTSUBCOUNT)*PTSUBSTEP
!
!* 1.3 Tecplotfile : opening + headers
IF (LTECOUTPTS) THEN
CALL OPEN_TECOUT(45, KTCOUNT, KTSUBCOUNT)
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(45, XPOS_TELT_RG(JROT,JTELT,:))
END DO
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(45, XPOS_NELT_RG(JROT,JNELT,:))
END DO
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(45, XPOS_HUB_RG(JROT,:))
END IF
!
!
! ---- BLADES ----
!
DO JBLA=1, INB_B
!* B.1 Update positions
XPOS_BLA_RG(JROT,JBLA,:) = XPOS_HUB_RG(JROT,:) &
+ MATMUL(XMAT_RG_RH(JROT,:,:),XPOSINI_BLA_RH(JROT,JBLA,:))
!
!* B.2 Update orientation
CALL GET_ORI_MAT_X(XANGINI_BLA_RH(JROT,JBLA,1) + XRVEL_RB_RH(JROT,JBLA,1)*ZTIME, ZORI_MAT_X)
CALL GET_ORI_MAT_Y(XANGINI_BLA_RH(JROT,JBLA,2) + XRVEL_RB_RH(JROT,JBLA,2)*ZTIME, ZORI_MAT_Y)
CALL GET_ORI_MAT_Z(XANGINI_BLA_RH(JROT,JBLA,3) + XRVEL_RB_RH(JROT,JBLA,3)*ZTIME, ZORI_MAT_Z)
! Orientation matrix
XMAT_RH_RB(JROT,JBLA,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z))
XMAT_RG_RB(JROT,JBLA,:,:) = MATMUL(XMAT_RG_RH(JROT,:,:), XMAT_RH_RB(JROT,JBLA,:,:))
!
!* B.3 Update structural velocities
! Rotation of blade in RG
XRVEL_RB_RG(JROT,JBLA,:) = XRVEL_RH_RG(JROT,:) &
+ MATMUL(XMAT_RG_RB(JROT,JBLA,:,:),XRVEL_RB_RH(JROT,JBLA,:))
! Translation of blade in RG
ZDIST_HUB_BLA_RG(:) = XPOS_BLA_RG(JROT,JBLA,:) - XPOS_HUB_RG(JROT,:)
XTVEL_BLA_RG(JROT,JBLA,:) = XTVEL_HUB_RG(JROT,:) &
+ CROSS(XRVEL_RB_RG(JROT,JBLA,:),ZDIST_HUB_BLA_RG(:))
!
!* B.4 Print in tecplot file
IF (LTECOUTPTS) THEN
CALL PRINT_TECOUT(45, XPOS_BLA_RG(JROT,JBLA,:))
END IF
!
!
! ---- ELEMENTS ----
!* E.0 Positioning sections (cuts) in RG
DO JBELT=1, INB_BELT+1
XPOS_SEC_RG(JROT,JBLA,JBELT,:) = XPOS_BLA_RG(JROT,JBLA,:) &
+ MATMUL(XMAT_RG_RB(JROT,JBLA,:,:),XPOS_SEC_RB(JROT,JBLA,JBELT,:))
ENDDO
!
!
DO JBELT=1, INB_BELT
!* E.1 Positioning sections centers (application points) in RG
XPOS_ELT_RG(JROT,JBLA,JBELT,:) = XPOS_BLA_RG(JROT,JBLA,:) &
+ MATMUL(XMAT_RG_RB(JROT,JBLA,:,:),XPOS_ELT_RB(JROT,JBLA,JBELT,:))
!* E.2 Update orientation
CALL GET_ORI_MAT_X(XANGINI_ELT_RB(JROT,JBLA,JBELT,1) &
+ XRVEL_RE_RB(JROT,JBLA,JBELT,1)*ZTIME, ZORI_MAT_X)
CALL GET_ORI_MAT_Y(XANGINI_ELT_RB(JROT,JBLA,JBELT,2) &
+ XRVEL_RE_RB(JROT,JBLA,JBELT,2)*ZTIME, ZORI_MAT_Y)
CALL GET_ORI_MAT_Z(XANGINI_ELT_RB(JROT,JBLA,JBELT,3) &
+ XRVEL_RE_RB(JROT,JBLA,JBELT,3)*ZTIME, ZORI_MAT_Z)
! Orientation matrix
XMAT_RB_RE(JROT,JBLA,JBELT,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y,ZORI_MAT_Z))
XMAT_RG_RE(JROT,JBLA,JBELT,:,:) = MATMUL(XMAT_RG_RB(JROT,JBLA,:,:), XMAT_RB_RE(JROT,JBLA,JBELT,:,:))
XMAT_RE_RG(JROT,JBLA,JBELT,:,:) = TRANSPOSE(XMAT_RG_RE(JROT,JBLA,JBELT,:,:))
!
!* E.3 Update structural velocities
! Rotation of elements in RG
XRVEL_RE_RG(JROT,JBLA,JBELT,:) = XRVEL_RB_RG(JROT,JBLA,:) &
+ MATMUL(XMAT_RG_RE(JROT,JBLA,JBELT,:,:),&
XRVEL_RE_RB(JROT,JBLA,JBELT,:))
! Translation of elements in RG
ZDIST_BLA_ELT_RG(:) = XPOS_ELT_RG(JROT,JBLA,JBELT,:) - XPOS_BLA_RG(JROT,JBLA,:)
XTVEL_ELT_RG(JROT,JBLA,JBELT,:) = XTVEL_BLA_RG(JROT,JBLA,:) &
+ CROSS(XRVEL_RE_RG(JROT,JBLA,JBELT,:),ZDIST_BLA_ELT_RG(:))
XTVEL_ELT_RE(JROT,JBLA,JBELT,:) = MATMUL(XMAT_RE_RG(JROT,JBLA,JBELT,:,:),&
XTVEL_ELT_RG(JROT,JBLA,JBELT,:))
!
!* E.4 Print in tecplot file
IF (LTECOUTPTS) THEN
CALL PRINT_TECOUT(45, XPOS_ELT_RG(JROT,JBLA,JBELT,:))
END IF
!
! ---- Leading Edge and Trailing Edge ----
! It is just to have a tecplot more beautiful
! For the moment, they are useless for computations
IF (LTECOUTPTS) THEN
! LE.1 Update Positions
ZPOS_ELTLE_RE(1) = 0d0
ZPOS_ELTLE_RE(2) = - 1d0/4d0 * XCHORD_ELT(JROT,JBLA,JBELT)
ZPOS_ELTLE_RE(3) = 0d0
ZPOS_ELTTE_RE(1) = 0d0
ZPOS_ELTTE_RE(2) = + 3d0/4d0 * XCHORD_ELT(JROT,JBLA,JBELT)
ZPOS_ELTTE_RE(3) = 0d0
!
ZPOS_ELTLE_RG(:) = XPOS_ELT_RG(JROT,JBLA,JBELT,:) &
+ MATMUL(XMAT_RG_RE(JROT,JBLA,JBELT,:,:), &
ZPOS_ELTLE_RE(:))
ZPOS_ELTTE_RG(:) = XPOS_ELT_RG(JROT,JBLA,JBELT,:) &
+ MATMUL(XMAT_RG_RE(JROT,JBLA,JBELT,:,:), &
ZPOS_ELTTE_RE(:))
!* LE.2 Print in tecplot file
CALL PRINT_TECOUT(45, ZPOS_ELTLE_RG(:))
CALL PRINT_TECOUT(45, ZPOS_ELTTE_RG(:))
END IF
!
!
END DO ! Blade element loop
END DO ! Blade loop
END DO ! Rotor loop
!
! Closing tec file
IF (LTECOUTPTS) THEN
CLOSE(45)
END IF
!
END SUBROUTINE EOL_KINE_ALM