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

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