1 : C DEVFB4 SOURCE CHAT 05/01/12 22:45:30 5004
2 : SUBROUTINE DEVFB4(ITYP,FTOTB,XPTB,IPALB,IPLIB,XPALB,XVALB,NLIAB,
3 : & NPLB,IND,IND1,INDM1,NPA,NPAM1,IND2,PDT,PDTS2,
4 : & FEXPSM,NPC1,I,IERRD,iannul)
5 : IMPLICIT INTEGER(I-N)
6 : IMPLICIT REAL*8(A-H,O-Z)
7 : *--------------------------------------------------------------------*
8 : * *
9 : * Op{rateur DYNE : algorithme de Fu - de Vogelaere *
10 : * ________________________________________________ *
11 : * *
12 : * Calcul des forces de choc sur base B pour les liaisons *
13 : * de type PROFIL_PROFIL_... *
14 : * *
15 : * Param}tres: *
16 : * *
17 : * e ITYP type de la liaison. *
18 : * es FTOTB Forces ext{rieures totalis{es sur la base B. *
19 : * e XPTB Tableau des d{placements des points *
20 : * e IPALB Renseigne sur la liaison. *
21 : * e IPLIB Tableau contenant les num{ros "DYNE" de la liaison. *
22 : * e XPALB Tableau contenant les param}tres de la liaison. *
23 : * es XVALB Tableau contenant les variables internes de liaisons. *
24 : * e NLIAB Nombre de liaisons sur la base B. *
25 : * e NPLB Nombre total de points intervenant dans les liaisons. *
26 : * e IND Indice du pas. *
27 : * e I num{ro de la liaison. *
28 : * es IERRD num{ro d'erreur. *
29 : * *
30 : * *
31 : * Auteur, date de cr{ation: *
32 : * *
33 : * Lionel VIVAN : le 1 f{vrier 1990 *
34 : * *
35 : *--------------------------------------------------------------------*
36 : *
37 : INTEGER IPALB(NLIAB,*),IPLIB(NLIAB,*)
38 : REAL*8 XPALB(NLIAB,*),XPTB(NPLB,4,*),FTOTB(NPLB,*)
39 : REAL*8 XVALB(NLIAB,4,*),XFOR(3),FEXPSM(NPLB,NPC1,2,*)
40 : LOGICAL LPETPP
41 : PARAMETER ( ZERO = 0.D0 , PRECIS = 1.D-15 )
42 : *
43 : * --- choc {l{mentaire PROFIL_PROFIL_INTERIEUR
44 : * --- choc {l{mentaire PROFIL_PROFIL_EXTERIEUR
45 : *
46 : IF (ITYP.EQ.31 .OR. ITYP.EQ.32) THEN
47 : NPOI = IPLIB(I,1)
48 : IDIM = IPALB(I,3)
49 : NOMBN1 = IPALB(I,4)
50 : NOMBN2 = IPALB(I,5)
51 : XRAID = XPALB(I,1)
52 : XSECT = XPALB(I,2)
53 : XEXPO = XPALB(I,3)
54 : ID1 = 3
55 : ID2 = ID1 + IDIM
56 : ID3 = ID1 + 2*IDIM
57 : ID4 = ID1 + 3*IDIM
58 : ID5 = ID1 + 4*IDIM
59 : ID6 = ID1 + 5*IDIM
60 : ID7 = ID6 + IDIM*NOMBN1
61 : ID8 = ID7 + IDIM*NOMBN2
62 : IP1 = 5
63 : IP2 = IP1 + NOMBN1
64 : IP3 = IP2
65 : IP4 = IP1 + NOMBN1 + NOMBN1*NOMBN2
66 : SOMFN = ZERO
67 : *
68 : * Recherche du nombre de point qui ont travers{ le maillage fixe,
69 : * donc s'il y a eu choc.
70 : IPOIN = 0
71 : ITIN = ID7
72 : DO 10 IN = 1,NOMBN2
73 : IP5 = IP4 + NOMBN1*(IN - 1)
74 : * D{placement du point IN du profil mobile dans le plan des profils
75 : XN1 = ZERO
76 : YN1 = ZERO
77 : DO 12 ID = 1,IDIM
78 : XXX = XPALB(I,ITIN+ID) + XPTB(NPOI,IND,ID)
79 : & - XPALB(I,ID2+ID)
80 : XN1 = XN1 + ( XXX * XPALB(I,ID3+ID) )
81 : YN1 = YN1 + ( XXX * XPALB(I,ID4+ID) )
82 : 12 CONTINUE
83 : * end do
84 : * On boucle sur le nombre d'{l{ment du profil fixe pour savoir si
85 : * le point IN @ travers{ le profil fixe
86 : NUME1 = ID6
87 : NUME2 = ID6 + IDIM
88 : ICOEF = ID8
89 : DO 14 IE = 1,NOMBN1
90 : IPALB(I,IP5+IE) = 0
91 : IPOSI = IPALB(I,IP3+IE)
92 : ICAS = IPALB(I,IP1+IE)
93 : XAIJ = XPALB(I,ICOEF+1)
94 : XBIJ = XPALB(I,ICOEF+2)
95 : IF (ICAS.EQ.1) THEN
96 : IF (XN1.GT.XAIJ) THEN
97 : IPOS2 = 1
98 : ELSE IF (XN1.LT.XAIJ) THEN
99 : IPOS2 = -1
100 : ELSE
101 : IPOS2 = 0
102 : ENDIF
103 : ELSE IF (ICAS.EQ.2) THEN
104 : IF (YN1.GT.XBIJ) THEN
105 : IPOS2 = 1
106 : ELSE IF (YN1.LT.XBIJ) THEN
107 : IPOS2 = -1
108 : ELSE
109 : IPOS2 = 0
110 : ENDIF
111 : ELSE
112 : RES = ( XAIJ * XN1 ) + XBIJ
113 : IF (YN1.GT.RES) THEN
114 : IPOS2 = 1
115 : ELSE IF (YN1.LT.RES) THEN
116 : IPOS2 = -1
117 : ELSE
118 : IPOS2 = 0
119 : ENDIF
120 : ENDIF
121 : IF (IPOSI.NE.IPOS2) THEN
122 : * Le point IN a travers{ le maillage fixe,
123 : * quel {l{ment a {t{ travers{ par ce point ?
124 : IF (IE.EQ.NOMBN1) NUME2 = ID6
125 : XE1 = ZERO
126 : YE1 = ZERO
127 : XE2 = ZERO
128 : YE2 = ZERO
129 : DO 16 ID = 1,IDIM
130 : XXX = XPALB(I,NUME1+ID) - XPALB(I,ID2+ID)
131 : YYY = XPALB(I,NUME2+ID) - XPALB(I,ID2+ID)
132 : XE1 = XE1 + ( XXX * XPALB(I,ID3+ID) )
133 : YE1 = YE1 + ( XXX * XPALB(I,ID4+ID) )
134 : XE2 = XE2 + ( YYY * XPALB(I,ID3+ID) )
135 : YE2 = YE2 + ( YYY * XPALB(I,ID4+ID) )
136 : 16 CONTINUE
137 : * end do
138 : IF (XE2.GT.XE1) THEN
139 : IF (XE1.LE.XN1 .AND. XN1.LT.XE2) THEN
140 : IPOIN = IPOIN + 1
141 : IPALB(I,IP5+IE) = 1
142 : GOTO 18
143 : ENDIF
144 : ELSE
145 : IF (XE2.LT.XN1 .AND. XN1.LE.XE1) THEN
146 : IPOIN = IPOIN + 1
147 : IPALB(I,IP5+IE) = 1
148 : GOTO 18
149 : ENDIF
150 : ENDIF
151 : IF (YE2.GT.YE1) THEN
152 : IF (YE1.LE.YN1 .AND. YN1.LT.YE2) THEN
153 : IPOIN = IPOIN + 1
154 : IPALB(I,IP5+IE) = 1
155 : GOTO 18
156 : ENDIF
157 : ELSE
158 : IF (YE2.LT.YN1 .AND. YN1.LE.YE1) THEN
159 : IPOIN = IPOIN + 1
160 : IPALB(I,IP5+IE) = 1
161 : GOTO 18
162 : ENDIF
163 : ENDIF
164 : ENDIF
165 : NUME1 = NUME2
166 : NUME2 = NUME2 + IDIM
167 : ICOEF = ICOEF + 2
168 : 14 CONTINUE
169 : * end do
170 : 18 CONTINUE
171 : IP3 = IP3 + NOMBN1
172 : ITIN = ITIN + IDIM
173 : 10 CONTINUE
174 : * end do
175 : *
176 : * Si le nombre de points ayant travers{s le profil fixe est {gale
177 : * au nombre de points du maillage mobile, alors on imprime une
178 : * erreur et on arrete le calcul.
179 : IF (IPOIN.EQ.NOMBN2) THEN
180 : IERRD = 1
181 : RETURN
182 : ENDIF
183 : *
184 : * s'il y a eu choc,
185 : * il y a IPOIN qui ont travers{ le maillage fixe
186 : *
187 : IF (IPOIN.NE.0) THEN
188 : *
189 : * Il y a eu choc :
190 : * - recherche du premier point qui a travers{ le maillage fixe, d'o:
191 : * d{termination du point P, intersection maillage fixe-mobile
192 : * - recherche du dernier point qui a travers{ le maillage fixe, d'o:
193 : * d{termination du point PP, intersection maillage fixe-mobile
194 : * - @ partir des points P et PP, d{termination de la normale de choc
195 : * - calcul de l'aire de choc, @ partir du point P, des points du
196 : * maillage fixe, du point PP et des points du maillage mobile qui
197 : * n'ont pas travers{ le maillage fixe
198 : * - calcul de la force de choc.
199 : *
200 : ICOMP = 0
201 : IN2 = 0
202 : 20 CONTINUE
203 : * On boucle sur les points du profil mobile
204 : IN2 = IN2 + 1
205 : IP5 = IP4 + NOMBN1*(IN2 - 1)
206 : IF (IN2.GT.NOMBN2) GOTO 22
207 : LPETPP = .FALSE.
208 : DO 24 IE1 = 1,NOMBN1
209 : IETAT = IPALB(I,IP5+IE1)
210 : * On regarde si le point IN2 a travers{ un {l{ment du profil fixe
211 : IF (IETAT.EQ.1) THEN
212 : ICOMP = ICOMP + 1
213 : IF (IN2.EQ.1) THEN
214 : * Si on traite le premier point, on cherche les points en partant du
215 : * dernier
216 : DO 26 INN2 = NOMBN2,1,-1
217 : IP6 = IP4 + NOMBN1*(INN2 - 1)
218 : DO 28 IEE1 = 1,NOMBN1
219 : IETA2 = IPALB(I,IP6+IEE1)
220 : IF (IETA2.EQ.1) THEN
221 : ICOMP = ICOMP + 1
222 : IEEE1 = IEE1
223 : GOTO 26
224 : ENDIF
225 : 28 CONTINUE
226 : * end do
227 : * Le point INN2 n'a pas travers{ le maillage fixe
228 : * Le point (INN2-1) et la point INN2 forment le segment de droite
229 : * pour d{finir le point d'intersection P avec l'{l{ment IEEE1
230 : IF (INN2.EQ.NOMBN2) THEN
231 : IPOIN1 = NOMBN2
232 : IPOIN2 = 1
233 : IELEM1 = IE1
234 : ELSE
235 : IPOIN1 = INN2
236 : IPOIN2 = INN2 + 1
237 : IELEM1 = IEEE1
238 : ENDIF
239 : GOTO 30
240 : 26 CONTINUE
241 : * end do
242 : 30 CONTINUE
243 : ELSE
244 : IPOIN1 = IN2 - 1
245 : IPOIN2 = IN2
246 : IELEM1 = IE1
247 : ENDIF
248 : * D{termination du point P
249 : * Le point IPOIN1 n'a pas travers{ le maillage fixe
250 : * Le point IPOIN2 a travers{ le maillage fixe par l'{l{ment IELEM1
250 : CALL DYNE31(IPOIN1,IPOIN2,IELEM1,XPALB,IPALB,XPTB,
251 : & NLIAB,NPLB,I,NPOI,IND,ID1,IP1,XP,YP)
252 : *
253 : * Recherche du point PP, on regarde les points suivant
254 : IF (IN2.EQ.NOMBN2) THEN
255 : IPOIN3 = IN2
256 : IPOIN4 = 1
257 : IELEM2 = IE1
258 : GOTO 32
259 : ENDIF
260 : IN22 = IN2 + 1
261 : IEEE1 = IE1
262 : DO 34 INN2 = IN22,NOMBN2
263 : IP7 = IP4 + NOMBN1*(INN2 - 1)
264 : DO 36 IEE1 = 1,NOMBN1
265 : IETA2 = IPALB(I,IP7+IEE1)
266 : IF (IETA2.EQ.1) THEN
267 : ICOMP = ICOMP + 1
268 : IEEE1 = IEE1
269 : IF (INN2.EQ.NOMBN2) THEN
270 : IPOIN3 = INN2
271 : IPOIN4 = 1
272 : IELEM2 = IEE1
273 : ENDIF
274 : GOTO 34
275 : ENDIF
276 : 36 CONTINUE
277 : * end do
278 : * Le point INN2 n'a pas travers{ le maillage fixe
279 : * Le point (INN2-1) et la point INN2 forment le segment de droite
280 : * pour d{finir le point d'intersection P avec l'{l{ment IEEE1
281 : IPOIN3 = INN2 - 1
282 : IPOIN4 = INN2
283 : IELEM2 = IEEE1
284 : GOTO 32
285 : 34 CONTINUE
286 : * end do
287 : 32 CONTINUE
288 : * D{termination du point PP
289 : * Le point IPOIN3 a travers{ le maillage fixe par l'{l{ment IELEM2
290 : * Le point IPOIN4 n'a pas travers{ le maillage fixe
290 : CALL DYNE31(IPOIN3,IPOIN4,IELEM2,XPALB,IPALB,XPTB,
291 : & NLIAB,NPLB,I,NPOI,IND,ID1,IP1,XPP,YPP)
292 : LPETPP = .TRUE.
293 : GOTO 40
294 : ENDIF
295 : 24 CONTINUE
296 : * end do
297 : 40 CONTINUE
298 : *
299 : IF ( LPETPP ) THEN
300 : * On a d{termin{ les points P et PP, calcul de la normale
301 : XNORM1 = ( XP + XPP ) * 0.5D0
302 : YNORM1 = ( YP + YPP ) * 0.5D0
303 : XXX = XP - XPP
304 : YYY = YP - YPP
305 : IPT3 = ID7 + IDIM*(IPOIN3-1)
306 : XN4 = ZERO
307 : YN4 = ZERO
308 : DO 42 ID = 1,IDIM
309 : XXX = XPALB(I,IPT3+ID) + XPTB(NPOI,IND,ID)
310 : & - XPALB(I,ID2+ID)
311 : XN4 = XN4 + ( XXX * XPALB(I,ID3+ID) )
312 : YN4 = YN4 + ( XXX * XPALB(I,ID4+ID) )
313 : 42 CONTINUE
314 : * END DO
315 : IF (ABS(YYY).LT.PRECIS) THEN
316 : YNN = YN4 - YNORM1
317 : XNORM = ZERO
318 : YNORM = YNN / ABS(YNN)
319 : ELSE IF (ABS(XXX).LT.PRECIS) THEN
320 : XNN = XN4 - XNORM1
321 : XNORM = XNN / ABS(XNN)
322 : YNORM = ZERO
323 : ELSE
324 : * le premier vecteur est P-PP
325 : * la normale de choc sera le produit vectoriel P_PP, N
326 : XPPP = XPP - XP
327 : YPPP = YPP - YP
328 : PS = SQRT( XPPP*XPPP + YPPP*YPPP )
329 : XNORM = YPPP * XPALB(I,ID5+3) / PS
330 : YNORM = -XPPP * XPALB(I,ID5+3) / PS
331 : ENDIF
332 : * Calcul de l'aire
332 : CALL DYNE32(IPOIN1,IPOIN4,IELEM1,IELEM2,XP,YP,XPP,YPP,
333 : & XPALB,IPALB,XPTB,NLIAB,NPLB,I,NPOI,ID1,IP1,IND,SOMAIR)
334 : AIRE = XSECT - SOMAIR
335 : * Calcul de la force de choc
336 : * Liaison conditionnelle
337 : IF (iannul.EQ.0) THEN
338 : FN = -XRAID * ( AIRE ** XEXPO )
339 : ELSE
340 : FN = 0
341 : ENDIF
342 : SOMFN = SOMFN + FN
343 : * On revient dans le repere global
344 : XFOR1 = FN * XNORM
345 : XFOR2 = FN * YNORM
346 : DO 38 ID = 1,IDIM
347 : XXX = XFOR1*XPALB(I,ID3+ID) + XFOR2*XPALB(I,ID4+ID)
348 : FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XXX
349 : 38 CONTINUE
350 : * end do
351 : IN2 = INN2
352 : ENDIF
353 : IF (ICOMP.EQ.IPOIN) GOTO 22
354 : GOTO 20
355 : 22 CONTINUE
356 : ENDIF
357 : DO 100 ID = 1,IDIM
358 : XVALB(I,IND,3+ID) = XPTB(NPOI,IND,ID)
359 : 100 CONTINUE
360 : * end do
361 : XVALB(I,IND,1) = SOMFN
362 : *
363 : * --- choc ...........
364 : *
365 : * ELSE IF (ITYP.EQ. ) THEN
366 : * .......
367 : * .......
368 : *
369 : ENDIF
370 : *
371 : END
372 :
373 :
© Cast3M 2003 - All rights reserved.
Disclaimer