Télécharger inosc2.eso

Retour à la liste

Numérotation des lignes :

  1. C INOSC2 SOURCE PV 15/04/13 21:15:11 8474
  2. CREAD INOSC1 SOURCE A1 GIBI 06/12/89 12:32:19
  3. C INOSC1 SOURCE CHAT 89/06/01 22:06:52
  4. SUBROUTINE INOSC2(IPT,IPG,DFREQ,XSI,RMAX,TMAX,
  5. 1 DIMALI,DUCTIL,
  6. 2 ALFA,BETA,GAMA,TETA,imod,idema,epse)
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8(A-H,O-Z)
  9. C
  10. C =====================================================
  11. C =
  12. C CALCUL DU MAXIMUM DE LA REPONSE A UN SIGNAL D'UN =
  13. C OSCILLATEUR NON LINEAIRE AVEC CONDITIONS INITIALES =
  14. C EN VITESSE ET EN DEPLACEMENT NULLES. =
  15. C =
  16. C CREATION : 15/06/90 =
  17. C PROGRAMMEUR : A.PINTO AND P.PEGON [ISPRA(I)] =
  18. C INSPIREE SUR LA SUBROUTINE 'INOSC1' DE 03/06/87 =
  19. c
  20. c I.P 05/98 ajout du modele bilineaire, iterations sur
  21. c la ductilite, changement du modèle Takeda
  22. C =
  23. c imod = 1 : modele takeda
  24. c 2 elastoplastique ecrouissage isotrope
  25. c 3 elastoplastique ecrouissage cinematique
  26. c 4 elastique nl
  27. C =====================================================
  28. C
  29. -INC CCOPTIO
  30. -INC SMLREEL
  31. -INC CCREEL
  32.  
  33. * quelques tableaux pour les variables internes du modèle Takeda
  34. REAL*8 VAR0(21),VARF(21)
  35. * qques tableaux mour utiliser le meme s_p que DYNE (dychec)
  36. REAL*8 XABSCI(1,4),XORDON(1,4)
  37.  
  38.  
  39.  
  40. MLREE1=IPT
  41. SEGACT MLREE1
  42. MLREE2=IPG
  43. SEGACT MLREE2
  44. NN=MLREE1.PROG(/1)
  45.  
  46. C
  47. C DIFFERENCE ENTRE LES 2 PREMIERS TEMPS SUCCESSIFS : DLTF
  48. C (PAS DE TEMP CONSTANT)
  49. C
  50. DLTF=MLREE1.PROG(2)-MLREE1.PROG(1)
  51. ** precision
  52. epsi = .05d0
  53. prec = 1d-3
  54. iprec = 0
  55. nitmax = 30
  56.  
  57. C
  58. C CALCUL DE CONSTANTES UTILES
  59. C
  60. W=2*XPI*DFREQ
  61. W2=W*W
  62.  
  63. C
  64. C CALCULO DOS COEFICIENTES DE IWAN
  65. C ================================
  66. DW=XSI*W
  67. D2=XSI*XSI
  68. A0=EXP(-DW*DLTF)
  69. A1=W*SQRT(1.0d0-D2)
  70. AD1=A1*DLTF
  71. A2=SIN(AD1)
  72. A3=COS(AD1)
  73. A4=(2.0d0*D2-1.0d0)/W2
  74. A5=XSI/W
  75. A6=2.0d0*A5/W2
  76. A7=1.0d0/W2
  77. A8=(A1*A3-DW*A2)*A0
  78. A9=-(A1*A2+DW*A3)*A0
  79. A10=A8/A1
  80. A30=A0/A1
  81. A12=A30*A2
  82. A13=A0*A3
  83. A14=A10*A4
  84. A15=A12*A4
  85. A16=A6*A13
  86. A17=A9*A6
  87. A11=A0*(DW*A2/A1+A3)
  88. A21=A10*DW+A9
  89. A22=A10
  90. B11=(-A15-A16+A6)/DLTF-A12*A5-A7*A13
  91. B12=(A15+A16-A6)/DLTF+A7
  92. B21=(-A14-A17-A7)/DLTF-A10*A5-A9*A7
  93. B22=(A14+A17+A7)/DLTF
  94.  
  95. ELAS0=W2
  96. ELAS=ELAS0
  97. ELAS1=ELAS0
  98. if (imod.eq.1 .and. teta.le.0d0) teta = .0001d0
  99. ELAS2=ELAS1*TETA
  100.  
  101. * qques initialisations pour etre compatible avec dychec
  102. xabsci(1,1) = 0.d0
  103. xordon(1,1) = 0.d0
  104. xjeu = 0d0
  105. iperm = imod
  106. if (imod.eq.4) iperm = -2
  107. nl = 1
  108. xamo = 0d0
  109. iannul = 0
  110. xvit = 0d0
  111. nip = 4
  112.  
  113. 400 continue
  114.  
  115. iterm = 0
  116. iterm2 = 0
  117. ibisec = 0
  118. 200 continue
  119. iterm = iterm + 1
  120. iterm2 = iterm2 + 1
  121. C
  122. C INICIALIZACAO DE VARIAVEIS
  123. C ==========================
  124. F00=0.D0
  125. F11=0.D0
  126. D0=0.D0
  127. V0=0.D0
  128. FNMDV1=0.D0
  129. FNMDV=0.D0
  130. DES0=0.D0
  131. FORC0=0.D0
  132. RMAX=0.D0
  133. epse= 0.d0
  134. DY=DIMALI/DUCTIL
  135. fy = dy*elas0
  136. ITER=0
  137. * qques parametres de Takeda
  138. dyn = -dy
  139. fyn = -fy
  140. fc = .999d0*fy
  141. fcn = .999d0*fyn
  142. pinp = beta*fy
  143. pinn = beta*fyn
  144. bp = gama/fy
  145. * qques parametres encore pour dychec
  146. if (imod.ne.1) then
  147. xabsci(1,2) = dy
  148. xordon(1,2) = fy
  149. if (elas2.ge.0d0) then
  150. xabsci(1,3) = dy + dy
  151. xordon(1,3) = fy + elas2*dy
  152. xabsci(1,4) = 3d0*dy
  153. xordon(1,4) = fy + elas2*dy*2d0
  154. else
  155. * si raideur negative on considere un plateau a .01fy
  156. xabsci(1,3) = -.9d0*fy/elas2 + dy
  157. xordon(1,3) = .1d0*fy
  158. xabsci(1,4) = 2d0*xabsci(1,3)
  159. xordon(1,4) = .01d0*fy
  160. endif
  161. xdplas = 0d0
  162. xdplac = 0d0
  163. endif
  164. C
  165. C CREATION DE LA BOUCLE
  166. C
  167. DO 100 I=1,NN
  168. C
  169. C INICIO DAS INTEGRACOES
  170. C ======================
  171. * force exterieure a la fin du pas precedent et courant
  172. F00=F11
  173. F11=-MLREE2.PROG(I)
  174. C INICIO DA CONVERGENCIA
  175. C ======================
  176. 999 ITER=ITER+1
  177. * on ajoute les forces nl au 2nd membre
  178. F0=F00-FNMDV
  179. F1=F11-FNMDV1
  180. D1=A11*D0+A12*V0+B11*F0+B12*F1
  181. V1=A21*D0+A22*V0+B21*F0+B22*F1
  182.  
  183. C
  184. C CALCULO DOS TERMOS NAO-LINEARES NO PASSO
  185. C ========================================
  186. DELTAD = D1 - D0
  187. deltad=max(xpetit,deltad)
  188.  
  189. if (imod.ne.1) then
  190.  
  191. call DYCHEC(d1,XDPLAS,XDPLAC,XJEU,IPERM,XABSCI,
  192. & XORDON,nl,xfla,nl,NIP,XVIT,XAMO,iannul)
  193. forca = -xfla
  194.  
  195. else
  196. * modele Takeda
  197. ** quelques initialisations
  198. forca = forc0
  199. if (iter.eq.1) then
  200. if (i.eq.1) then
  201. do 1 ii = 1,21
  202. var0(ii) = 0.d0
  203. 1 continue
  204. else
  205. do 2 m = 1,21
  206. var0(m) = VARF(m)
  207. 2 continue
  208. endif
  209. endif
  210.  
  211. DO 3 k=1,21
  212. VARF(k)=VAR0(k)
  213. 3 CONTINUE
  214. IFR1=INT(VARF(1))
  215.  
  216. IF(IFR1.EQ.0) THEN
  217. IFR1 = 1
  218. VARF(2)=elas0
  219. ENDIF
  220. DP1 = DELTAD * VARF(2)
  221.  
  222. CALL DDNSTH(IFR1,IFC1,DP1,FORCA,VARF(21),
  223. * VARF(2),elas0,elas2,elas2,fc,fcn,fy,fyn,
  224. * dy,dyn,alfa,alfa,pinp,pinn,bp,bp,
  225. * VARF(3),VARF(4),VARF(5),
  226. * VARF(6),VARF(7),VARF(8),VARF(9),
  227. * VARF(10),VARF(11),VARF(12),VARF(13),
  228. * VARF(14),VARF(15),VARF(16),VARF(17),
  229. * VARF(18),VARF(19),VARF(20))
  230.  
  231. VARF(1)=DBLE(IFR1)
  232.  
  233. endif
  234.  
  235. C ANALISE DA CONVERGENCIA EM FORCA
  236. C ================================
  237. * forces non lineaires à la fin de l'iteration courante
  238. FICT=-D1*ELAS0+FORCA
  239. ERRO=ABS((d1-des0)/deltad)
  240. FNMDV1=FICT
  241. DES0=D1
  242.  
  243. C
  244. C ANALISE DE CONVERGENCIA GLOBAL
  245. C ==============================
  246. IF(ERRO.gt.prec .and. ITER.lt.200) goto 999
  247.  
  248. * valeurs a la fin du pas
  249. ITER=0
  250. depse = abs(deltad - (forca - forc0)/elas0)
  251. epse = epse + depse
  252. FNMDV=FNMDV1
  253. D0=D1
  254. V0=V1
  255. forc0 = forca
  256.  
  257. IF(ABS(D1).GT.RMAX)THEN
  258. RMAX=ABS(D1)
  259. TMAX=MLREE1.PROG(I)
  260. ENDIF
  261.  
  262. IF(ITER.GT.200)THEN
  263. print*,'SPON: temps :',MLREE1.PROG(I),'nb d iterations > 200'
  264. ENDIF
  265. C
  266. 100 CONTINUE
  267.  
  268. if (idema.eq.0) goto 300
  269.  
  270. duct2 = rmax/dy
  271. *
  272. ** petite correction de la limite elastique si elle est trop petite
  273. ** (ca arrive pour des frequences elevees)
  274. ** pour une frequence tres grande ca revient a commencer les iterations
  275. ** avec une ductilite de 1./.9
  276.  
  277. c if (iterm2.eq.1) print*,'dfreq',dfreq
  278. c print*,'duct2=',duct2,'iterm=',iterm,'iterm2=',iterm2,
  279. c & 'freq=',dfreq
  280.  
  281. c print*,dimali,duct2
  282.  
  283.  
  284. if (duct2.gt. (5d0*ductil).and.iterm.eq.1
  285. & .and. (abs(teta)).lt.1d0) then
  286. if (iterm2.eq.1 .or.
  287. & (duct2.gt. (6d0*ductil).and.iterm2.le.3)) then
  288. iterm = iterm - 1
  289. dimali = dimali*(.9d0*ductil -
  290. & ((.9d0*ductil - 1)* exp(-.02*(duct2-5d0*ductil))))
  291. goto 200
  292. endif
  293. endif
  294.  
  295. diaux1 = dimali
  296.  
  297. if ((abs(duct2 - ductil)).gt.
  298. & (epsi*ductil).and.iterm2.le.nitmax) then
  299. if (iterm.eq.1) dimali = rmax
  300. if (iterm.ge.2) then
  301. aa = (duct2 - ductil)*(duct1-ductil)
  302. if (aa.gt.0.d0. and. ibisec.eq.0) then
  303. dimali = rmax
  304. else if (ibisec.eq.0) then
  305. * premiere dichotomie
  306. ibisec = 1
  307. dimali = (dimali + dimal1)/2d0
  308. ducaux = duct1
  309. diaux2 = dimal1
  310. else
  311. * dichotomies suivantes
  312. aa2 = (duct2 - ductil)*(duct3-ductil)
  313. if (aa.lt.0.d0) then
  314. dimali = (dimali + dimal1)/2d0
  315. ducaux = duct1
  316. diaux2 = dimal1
  317. else if (aa2.lt.0.d0) then
  318. dimali = (dimali + dimal3)/2d0
  319. ducaux = duct3
  320. diaux2 = dimal3
  321. else
  322. * dans le cas ou la dichotomie ne marche pas
  323. dimali = rmax
  324. endif
  325. endif
  326. endif
  327. duct1 = duct2
  328. dimal1 = diaux1
  329. if (ibisec.eq.1) then
  330. duct3 = ducaux
  331. dimal3 = diaux2
  332. endif
  333. goto200
  334. else if (iterm2.gt.nitmax) then
  335. * on essaie 2 fois encore avec un critere plus petit
  336. if (iprec.le.1) then
  337. prec = prec/10d0
  338. iprec = iprec + 1
  339. goto 400
  340. endif
  341. REAERR(1)=dfreq
  342. REAERR(2)=duct2
  343. CALL ERREUR(-326)
  344. endif
  345. rmax = dy
  346.  
  347. 300 continue
  348. SEGDES MLREE1
  349. SEGDES MLREE2
  350.  
  351.  
  352. RETURN
  353. END
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales