1 : C DEVLFA SOURCE CHAT 05/01/12 22:46:27 5004
2 : SUBROUTINE DEVLFA(Q1,Q2,FTOTA,XM,NA1,IPALA,IPLIA,XPALA,XVALA,
3 : & NLIAA,XDT,NPAS,IND,FINERT,IVINIT,FTEST,FTOTA0)
4 : IMPLICIT INTEGER(I-N)
5 : IMPLICIT REAL*8(A-H,O-Z)
6 : *--------------------------------------------------------------------*
7 : * *
8 : * Operateur DYNE : algorithme de Fu - de Vogelaere *
9 : * ________________________________________________ *
10 : * *
11 : * Calcul des forces de choc base A. *
12 : * *
13 : * Param}tres: *
14 : * *
15 : * e Q1(.,.) Vecteur des deplacements generalises. *
16 : * e Q2(.,.) Vecteur des vitesses generalisees. *
17 : * es FTOTA Forces ext{rieures totalis{es sur la base A. *
18 : * e XM Vecteur des masses g{n{ralis{es. *
19 : * e NA1 Nombre total d'inconnues en base A. *
20 : * e IPALA Renseigne sur le type de la liaison. *
21 : * e IPLIA Tableau contenant les num{ros "DYNE" de la liaison. *
22 : * e XPALA Tableau contenant les param}tres de la liaison. *
23 : * es XVALA Tableau contenant les variables internes des liaisons *
24 : * e NLIAA Nombre de liaisons sur la base A. *
25 : * e XDT Valeurs des pas de temps *
26 : * e NPAS Numero du pas de temps *
27 : * e IND Indice du demi-pas de temps *
28 : * = 2 si 1er demi-pas = 1 si 2eme demi-pas *
29 : * es FINERT Forces d'inertie *
30 : * e IVINIT =1 si vitesses initiales, =0 sinon *
31 : * *
32 : * Auteur, date de cr{ation: *
33 : * *
34 : * Lionel VIVAN, le 20 aout 1989. *
35 : * E. de LANGRE, 15 fevrier 91 correction jeu negatif *
36 : * Lionel VIVAN, 12 mars 1991 ajout de la liaison FLUIDE *
37 : * Denis ROBERT, 30 avril 1992 ajout liaison POLYNOMIALE *
38 : * *
39 : * remarque: *
40 : * ========= *
41 : * Si jeu n{gatif (cas particulier de la base A ou il n'y a pas de *
42 : * normale), on renverse les variables avec XNORM. *
43 : * *
44 : *--------------------------------------------------------------------*
45 : *
46 : INTEGER IPALA(NLIAA,*),IPLIA(NLIAA,*)
47 : REAL*8 XPALA(NLIAA,*),Q1(NA1,*),Q2(NA1,*),FTOTA(NA1,*),XM(*)
48 : REAL*8 XVALA(NLIAA,4,*),XDT(*),FINERT(NA1,*)
49 : PARAMETER (XZERO = 0.D0)
50 : REAL*8 FTest(nA1,4)
51 : REAL*8 ftOTa0(NA1,4)
52 : *
53 : * Boucle sur le nombre de liaisons
54 : *
55 : XFIN = 0.D0
56 : PDT = XDT(NPAS)
57 : PDTS2 = 0.5D0 * PDT
58 : DO 10 I = 1,NLIAA
59 : ITYP = IPALA(I,1)
60 : icond= IPALa(I,2)
61 : iannul= 0
62 : * --- cas des liaisons conditionnelles
63 : if (icond .eq. 1 ) then
64 : DO 101 j = 4,20
65 : jliai = ipala(i,j)
66 : jpliai = abs ( jliai)
67 : if ( jliai . EQ. 0 ) then
68 : goto 101
69 : else
70 : jtyp = ipala(jpliai,1)
71 : do 102 ik = 1,nA1
72 : do 103 jk = 1,4
73 : ftest(ik,jk) = 0d0
74 : ftota0 (ik,jk) = ftota(ik,jk)
75 : 103 continue
76 : 102 continue
77 : * CALCUL DES LIAISONS POUR TEST
78 : *
79 : * ****** choc {l{mentaire POINT_PLAN sans amortissement
80 : *
81 : IF (JTYP.EQ.1) THEN
82 : XRAID = XPALA(jpliai,1)
83 : XJEU = XPALA(jpliai,2)
84 : XNORM = 1.
85 : IF (XJEU.LT.0D0) THEN
86 : XNORM = -1.
87 : XJEU = -1*XJEU
88 : ENDIF
89 : INA1 = IPLIA(jpliai,1)
90 : XDEP = XNORM*Q1(INA1,IND)
90 : CALL DYCHEL(XDEP,XRAID,XJEU, XFL,iannul)
91 : XVALA(jpliai,IND,1) = XNORM*XFL
92 : XVALA(jpliai,IND,4) = XNORM*XDEP
93 : FTest(INA1,IND) =Ftest(INA1,IND) + XNORM*XFL
94 : *
95 : * ------ choc {l{mentaire POINT_PLAN avec amortissement
96 : *
97 : ELSE IF (JTYP.EQ.2) THEN
98 : XRAID = XPALA(jpliai,1)
99 : XJEU = XPALA(jpliai,2)
100 : XAMO = XPALA(jpliai,3)
101 : XNORM = 1.
102 : IPERM = 0
103 : IF (XJEU.LT.XZERO) THEN
104 : XNORM = -1.
105 : XJEU = -1*XJEU
106 : ENDIF
107 : INA1 = IPLIA(jpliai,1)
108 : XDEP = XNORM * Q1(INA1,IND)
109 : IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
110 : XVIT = XNORM * Q2(INA1,IND)
111 : ELSE
112 : IND2 = IND + 1
113 : XDEPM1 = XNORM * Q1(INA1,IND2)
114 : XVIT = (XDEP - XDEPM1) / PDTS2
115 : ENDIF
116 : XVALA(jpliai,IND,3) = XNORM*XVIT
116 : CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFL,IPERM,iannul)
117 : *+*
118 : IF (XDEP.GT.0.D0 .AND. XFL.GT.0.D0) XFL = 0.D0
119 : IF (XDEP.LT.0.D0 .AND. XFL.LT.0.D0) XFL = 0.D0
120 : *+*
121 : XVALA(jpliai,IND,1) = XNORM*XFL
122 : XVALA(jpliai,IND,4) = XNORM*XDEP
123 : FTest(INA1,IND) =Ftest(INA1,IND) + XNORM*XFL
124 : *
125 : * ------ choc {l{mentaire POINT_PLAN_FLUIDE
126 : *
127 : ELSE IF (JTYP.EQ.3) THEN
128 : XINER = XPALA(jpliai,1)
129 : XCONV = XPALA(jpliai,2)
130 : XVISC = XPALA(jpliai,3)
131 : XPCEL = XPALA(jpliai,4)
132 : XPCRA = XPALA(jpliai,5)
133 : XJEU = XPALA(jpliai,6)
134 : INA1 = IPLIA(jpliai,1)
135 : XDEP = Q1(INA1,IND)
136 : IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
137 : XVIT = Q2(INA1,IND)
138 : ELSE
139 : IND2 = IND + 1
140 : XDEPM1 = Q1(INA1,IND2)
141 : XVIT = (XDEP - XDEPM1) / PDTS2
142 : ENDIF
143 : IF (XJEU.GT.0D0) THEN
144 : XDH= XJEU - XDEP
145 : XNORM = 1.
146 : ELSE
147 : XDH= XDEP - XJEU
148 : XNORM = -1.
149 : ENDIF
150 : * Calcul de la masse ajout{e
151 : XXIN = -XINER / XDH
152 : FINERT(INA1,IND) = FINERT(INA1,IND) + XXIN
153 : * Calcul de la force de convection
153 : CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFCO,iannul)
154 : * Calcul de la force de viscosit{
154 : CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFVI,iannul)
155 : * Calcul de la force de perte de charge
155 : CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFPE,iannul)
156 : XFL = (XFCO* XNORM ) + XFVI + XFPE
157 : XVALA(jpliai,IND,1) = XDEP
158 : XVALA(jpliai,IND,2) = XVIT
159 : XVALA(jpliai,IND,3) = XXIN
160 : XVALA(jpliai,IND,4) = XFCO*XNORM
161 : XVALA(jpliai,IND,5) = XFVI
162 : XVALA(jpliai,IND,6) = XFPE
163 : FTest(INA1,IND) =Ftest(iNA1,IND) + XFL
164 : *
165 : * ------ force {l{mentaire de COUPLAGE EN VITESSE
166 : *
167 : ELSE IF (JTYP.EQ.4) THEN
168 : INA1 = IPLIA(jpliai,1)
169 : INA2 = IPLIA(jpliai,2)
170 : XDEP = Q1(INA2,IND)
171 : IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
172 : XVIT = Q2(INA2,IND)
173 : ELSE
174 : IND2 = IND + 1
175 : XDEPM1 = Q1(INA2,IND2)
176 : XVIT = (XDEP - XDEPM1) / PDTS2
177 : ENDIF
178 : XCPLGE = XPALA(jpliai,1)
178 : CALL DYCPLV(XVIT,XCPLGE,XFL,iannul)
179 : XVALA(jpliai,IND,3) = XVIT
180 : XVALA(jpliai,IND,1) = XFL
181 : XVALA(jpliai,IND,4) = XDEP
182 : FTest(INA1,IND) =Ftest(INA1,IND) + XFL
183 : *
184 : * ------ force {l{mentaire de COUPLAGE EN DEPLACEMENT
185 : *
186 : ELSE IF (JTYP.EQ.5) THEN
187 : INA1 = IPLIA(jpliai,1)
188 : INA2 = IPLIA(jpliai,2)
189 : XDEP = Q1(INA2,IND)
190 : XCPLGE = XPALA(jpliai,1)
190 : CALL DYCPLD(XDEP,XCPLGE,XFL,iannul)
191 : XVALA(jpliai,IND,1) = XFL
192 : XVALA(jpliai,IND,4) = XDEP
193 : FTest(INA1,IND) =Ftest(INA1,IND) + XFL
194 : *
195 : * ------ force elementaire de type POLYNOMIALE
196 : *
197 : ELSE IF (JTYP.EQ.6) THEN
198 : * WRITE(IOIMP,*)'DEVLFA : IPLIA(',jpliai,'1)=',IPLIA(jpliai,1)
199 : * nombre de modes "origine"
200 : NMOD = IPALA(I,2)
200 : CALL DYPOL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,XDT,
201 : & NPAS,jpliai,NMOD,FTest,IVINIT,iannul)
202 : *
203 : * ------ choc {l{mentaire ...
204 : *
205 : * ELSE IF (JTYP.EQ. ) THEN
206 : * .......
207 : * .......
208 : *
209 : ENDIF
210 : ***** tests pour voir si on annulera
211 : xff = 0d0
212 : do 104 ik = 1,na1
213 : do 105 jk = 1,4
214 : xff = xff + ( ftest(ik,jk) ** 2)
215 : 105 continue
216 : 104 continue
217 : xff = xff ** .5
218 : if ( ((xff .le. 1e-20 ) .and. ( jliai .gt. 0) )
219 : & .OR. ((xff .gt. 1e-20 ) .and. ( jliai .lt. 0) ) )
220 : & then
221 : iannul = 1
222 : endif
223 : endif
224 : 101 continue
225 : * --- fin du cas des liaisons conditionnelles
226 : endif
227 : *
228 : * ------ choc {l{mentaire POINT_PLAN sans amortissement
229 : *
230 : IF (ITYP.EQ.1) THEN
231 : XRAID = XPALA(I,1)
232 : XJEU = XPALA(I,2)
233 : XNORM = 1.
234 : IF (XJEU.LT.0D0) THEN
235 : XNORM = -1.
236 : XJEU = -1*XJEU
237 : ENDIF
238 : INA1 = IPLIA(I,1)
239 : XDEP = XNORM*Q1(INA1,IND)
239 : CALL DYCHEL(XDEP,XRAID,XJEU, XFL,iannul)
240 : XVALA(I,IND,1) = XNORM*XFL
241 : XVALA(I,IND,4) = XNORM*XDEP
242 : FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL
243 : *
244 : * ------ choc {l{mentaire POINT_PLAN avec amortissement
245 : *
246 : ELSE IF (ITYP.EQ.2) THEN
247 : XRAID = XPALA(I,1)
248 : XJEU = XPALA(I,2)
249 : XAMO = XPALA(I,3)
250 : XNORM = 1.
251 : IPERM = 0
252 : IF (XJEU.LT.XZERO) THEN
253 : XNORM = -1.
254 : XJEU = -1*XJEU
255 : ENDIF
256 : INA1 = IPLIA(I,1)
257 : XDEP = XNORM * Q1(INA1,IND)
258 : IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
259 : XVIT = XNORM * Q2(INA1,IND)
260 : ELSE
261 : IND2 = IND + 1
262 : XDEPM1 = XNORM * Q1(INA1,IND2)
263 : XVIT = (XDEP - XDEPM1) / PDTS2
264 : ENDIF
265 : XVALA(I,IND,3) = XNORM*XVIT
265 : CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFL,IPERM,iannul)
266 : *+*
267 : IF (XDEP.GT.0.D0 .AND. XFL.GT.0.D0) XFL = 0.D0
268 : IF (XDEP.LT.0.D0 .AND. XFL.LT.0.D0) XFL = 0.D0
269 : *+*
270 : XVALA(I,IND,1) = XNORM*XFL
271 : XVALA(I,IND,4) = XNORM*XDEP
272 : FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL
273 : *
274 : * ------ choc {l{mentaire POINT_PLAN_FLUIDE
275 : *
276 : ELSE IF (ITYP.EQ.3) THEN
277 : XINER = XPALA(I,1)
278 : XCONV = XPALA(I,2)
279 : XVISC = XPALA(I,3)
280 : XPCEL = XPALA(I,4)
281 : XPCRA = XPALA(I,5)
282 : XJEU = XPALA(I,6)
283 : INA1 = IPLIA(I,1)
284 : XDEP = Q1(INA1,IND)
285 : IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
286 : XVIT = Q2(INA1,IND)
287 : ELSE
288 : IND2 = IND + 1
289 : XDEPM1 = Q1(INA1,IND2)
290 : XVIT = (XDEP - XDEPM1) / PDTS2
291 : ENDIF
292 : IF (XJEU.GT.0D0) THEN
293 : XDH= XJEU - XDEP
294 : XNORM = 1.
295 : ELSE
296 : XDH= XDEP - XJEU
297 : XNORM = -1.
298 : ENDIF
299 : * Calcul de la masse ajout{e
300 : XXIN = -XINER / XDH
301 : FINERT(INA1,IND) = FINERT(INA1,IND) + XXIN
302 : * Calcul de la force de convection
302 : CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFCO,iannul)
303 : * Calcul de la force de viscosit{
303 : CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFVI,iannul)
304 : * Calcul de la force de perte de charge
304 : CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFPE,
305 : & iannul)
306 : XFL = (XFCO* XNORM ) + XFVI + XFPE
307 : XVALA(I,IND,1) = XDEP
308 : XVALA(I,IND,2) = XVIT
309 : XVALA(I,IND,3) = XXIN
310 : XVALA(I,IND,4) = XFCO*XNORM
311 : XVALA(I,IND,5) = XFVI
312 : XVALA(I,IND,6) = XFPE
313 : FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
314 : *
315 : * ------ force {l{mentaire de COUPLAGE EN VITESSE
316 : *
317 : ELSE IF (ITYP.EQ.4) THEN
318 : INA1 = IPLIA(I,1)
319 : INA2 = IPLIA(I,2)
320 : XDEP = Q1(INA2,IND)
321 : IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
322 : XVIT = Q2(INA2,IND)
323 : ELSE
324 : IND2 = IND + 1
325 : XDEPM1 = Q1(INA2,IND2)
326 : XVIT = (XDEP - XDEPM1) / PDTS2
327 : ENDIF
328 : XCPLGE = XPALA(I,1)
328 : CALL DYCPLV(XVIT,XCPLGE,XFL,iannul)
329 : XVALA(I,IND,3) = XVIT
330 : XVALA(I,IND,1) = XFL
331 : XVALA(I,IND,4) = XDEP
332 : FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
333 : *
334 : * ------ force {l{mentaire de COUPLAGE EN DEPLACEMENT
335 : *
336 : ELSE IF (ITYP.EQ.5) THEN
337 : INA1 = IPLIA(I,1)
338 : INA2 = IPLIA(I,2)
339 : XDEP = Q1(INA2,IND)
340 : XCPLGE = XPALA(I,1)
340 : CALL DYCPLD(XDEP,XCPLGE,XFL,iannul)
341 : XVALA(I,IND,1) = XFL
342 : XVALA(I,IND,4) = XDEP
343 : FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
344 : *
345 : * ------ force elementaire de type POLYNOMIALE
346 : *
347 : ELSE IF (ITYP.EQ.6) THEN
348 : * WRITE(IOIMP,*)'DEVLFA : IPLIA(',I,'1)=',IPLIA(I,1)
349 : * nombre de modes "origine"
350 : NMOD = IPALA(I,2)
350 : CALL DYPOL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,XDT,
351 : & NPAS,I,NMOD,FTOTA,IVINIT,iannul)
352 : *
353 : * ------ choc {l{mentaire ...
354 : *
355 : * ELSE IF (ITYP.EQ. ) THEN
356 : * .......
357 : * .......
358 : *
359 : ENDIF
360 :
361 :
362 : *
363 : * la suite n'est plus utile car on passe iannul aux s_p de calcul des
364 : * forces de liaisons.
365 :
366 : * si la liaisn était annulée on l'annule
367 : * if ( ( icond.eq. 1 ) .and. ( iannul.eq.1)) then
368 : * on annulle l'increment de ftotb
369 : * do 112 ik = 1,na1
370 : * do 113 jk = 1,4
371 : * ftota (ik,jk) = ftota0(ik,jk)
372 : * 113 continue
373 : * 112 continue
374 :
375 : * end if
376 : 10 CONTINUE
377 : * END DO
378 : *
379 : END
380 :
381 :
382 :
383 :
384 :
© Cast3M 2003 - All rights reserved.
Disclaimer