Télécharger issgra.eso

Retour à la liste

Numérotation des lignes :

  1. C ISSGRA SOURCE PV 17/12/08 21:17:29 9660
  2. C
  3. SUBROUTINE ISSGRA(WRK52,WRK53,WRK54,WRK27,IB,IGAU,NBPGAU)
  4. C
  5. C______________________________________________________________________
  6. C
  7. C MODELE ISSGRA
  8. C______________________________________________________________________
  9. C
  10. IMPLICIT INTEGER(I-N)
  11. IMPLICIT REAL*8(A-H,O-Z)
  12. c inutile? -INC CCOPTIO
  13. -INC CCREEL
  14. -INC DECHE
  15. *
  16. C GLOBAL MODEL
  17. DIMENSION XMKel(5,5)
  18. C DIMENSION XMCel(5,5)
  19. DIMENSION uo(5), duo(5), Fo(5)
  20. DIMENSION XMK(5,5), Ftr(5), force(5)
  21. *
  22. DIMENSION upl(5), uel(5)
  23. DIMENSION uoi(5), duoi(5), dueli(5), force_preci(5)
  24. C DIMENSION ueli(5)
  25.  
  26. C Modèle de Plasticité
  27. DIMENSION ddlambda(2), yield(2)
  28. DIMENSION Ho(2,2), H(2,2), Hinv(2,2)
  29. DIMENSION duP(5)
  30. *
  31. C Modèle de plasticité 1
  32. DIMENSION q1(5), upl1(5)
  33. C DIMENSION q1_prec(5)
  34. DIMENSION Gradgs(5),Gradfs(5),Gradfq(5)
  35. DIMENSION hh1(5), XKGradgs(5)
  36. *
  37. C Modèle de plasticité vertical
  38. DIMENSION upl2(5), Gradf2s(5), XKGradf2s(5)
  39. *
  40. C Modèle de Décollement
  41. DIMENSION dforce(5), forcepred(5), force_prec(5), duel(5)
  42. *
  43. C Amortissment
  44. DIMENSION upl_preci(5),q1_preci(5)
  45. *
  46. DIMENSION da(20)
  47. C DIMENSION TIME(2)
  48. *
  49. C-----------------------------------------------------------------------
  50. C-----------------------------------------------------------------------
  51. C Affectation des données du modèle
  52. da(1) =XMAT(13)
  53. da(2) =XMAT(14)
  54. da(3) =XMAT(15)
  55. da(4) =XMAT(7)
  56. da(5) =XMAT(8)
  57. da(6) =XMAT(12)
  58. da(7) =XMAT(11)
  59. da(20)=XMAT(10)
  60. da(8) =XMAT(16)
  61. da(9) =XMAT(17)
  62. da(10)=XMAT(18)
  63. da(11)=XMAT(19)
  64. da(12)=XMAT(20)
  65. da(13)=XMAT(21)
  66. da(14)=XMAT(22)
  67. da(15)=XMAT(23)
  68. da(16)=XMAT(24)
  69. da(17)=XMAT(25)
  70. da(18)=XMAT(26)
  71. da(19)=XMAT(27)
  72. *
  73. Diam=da(1)
  74. XLx=da(2)
  75. XLy=da(3)
  76. XKelz=da(4)
  77. XKelh=da(5)
  78. XKelry=da(6)
  79. XKelrx=da(7)
  80. C a=da(8)
  81. C b=da(9)
  82. C c=da(10)
  83. C d=da(11)
  84. C e=da(12)
  85. C f=da(13)
  86. qmax=da(14)
  87. a9=da(15)
  88. C a6=da(16)
  89. eta3=da(17)
  90. dtime=da(18)
  91. a8=da(19)
  92. XKelt=da(20)
  93. *
  94. *
  95. upl1(1)=VAR0(2)
  96. upl1(2)=VAR0(3)
  97. upl1(3)=VAR0(4)
  98. upl1(4)=VAR0(5)
  99. upl1(5)=VAR0(6)
  100. q1(1)=VAR0(7)
  101. q1(2)=VAR0(8)
  102. q1(3)=VAR0(9)
  103. q1(4)=VAR0(10)
  104. q1(5)=VAR0(11)
  105. force_prec(1)=-SIG0(1)
  106. force_prec(2)=-SIG0(2)
  107. force_prec(3)=SIG0(6)
  108. force_prec(4)=SIG0(3)
  109. force_prec(5)=-SIG0(5)
  110. delta_prec=VAR0(12)
  111. Xkapa_prec=VAR0(13)
  112. Xksi_prec=VAR0(14)
  113. distmin_prec=VAR0(15)
  114. upl2(1)=VAR0(17)
  115. upl2(2)=VAR0(18)
  116. upl2(3)=VAR0(19)
  117. upl2(4)=VAR0(20)
  118. upl2(5)=VAR0(21)
  119. q2=VAR0(22)
  120. DO I=1,5
  121. upl(I)=upl1(I)+upl2(I)
  122. ENDDO
  123. *
  124. *deplx joint --> -deplz ISS
  125. uo(1)=-epst0(1)
  126. *deply joint --> -deplx ISS
  127. uo(2)=-epst0(2)
  128. *rotz joint --> -roty ISS
  129. uo(3)=-epst0(6)
  130. *deplz joint --> deply ISS
  131. uo(4)=epst0(3)
  132. *roty joint --> -rotx ISS
  133. uo(5)=-epst0(5)
  134. *
  135. duo(1)=-DEPST(1)
  136. duo(2)=-DEPST(2)
  137. duo(3)=-DEPST(6)
  138. duo(4)=DEPST(3)
  139. duo(5)=-DEPST(5)
  140.  
  141. C Correction données d'entrées Castem
  142. DO I=1,5
  143. uo(I)=uo(I)+duo(I)
  144. * epstf(I) = epst0(I)+DEPST(I)
  145. ENDDO
  146. *
  147. C on adimensionne les déplacements
  148. if (a9.eq.2) then
  149. DO I=1,5
  150. uo(I) = uo(I)*sqrt(XLx**2+XLy**2)/(XLx*XLy)
  151. ENDDO
  152. uo(3) = uo(3)*XLx
  153. uo(5) = uo(5)*XLy
  154. DO I=1,5
  155. duo(I) = duo(I)*sqrt(XLx**2+XLy**2)/(XLx*XLy)
  156. ENDDO
  157. duo(3)= duo(3)*XLx
  158. duo(5)= duo(5)*XLy
  159. else
  160. DO I=1,5
  161. uo(I) = uo(I)/Diam
  162. ENDDO
  163. uo(3) = uo(3)*Diam
  164. uo(5) = uo(5)*Diam
  165. DO I=1,5
  166. duo(I) = duo(I)/Diam
  167. ENDDO
  168. duo(3)= duo(3)*Diam
  169. duo(5)= duo(5)*Diam
  170. endif
  171. *
  172. C Adimensionnement des efforts
  173. if (a9.eq.1) then
  174. C semelle filante
  175. DO I=1,5
  176. force_prec(I)=force_prec(I)/(Diam*qmax)
  177. ENDDO
  178. force_prec(3)=force_prec(3)/Diam
  179. force_prec(5)=force_prec(5)/Diam
  180. else
  181. if (a9.eq.2) then
  182. C semelle rectangulaire
  183. DO I=1,5
  184. force_prec(I)=force_prec(I)/(XLx*XLy*qmax)
  185. ENDDO
  186. force_prec(3)=force_prec(3)/XLx
  187. force_prec(5)=force_prec(5)/XLy
  188. else
  189. C semelle circulaire
  190. DO I=1,5
  191. force_prec(I)=force_prec(I)/(XPI*Diam**2/4*qmax)
  192. ENDDO
  193. force_prec(3)=force_prec(3)/Diam
  194. force_prec(5)=force_prec(5)/Diam
  195. endif
  196. endif
  197. C ***********************************************************************
  198. C Macro-élément d'ISS
  199. C adimensionnement de la matrice de rigidité élastique suivant la forme de la fondation (paramètre a9)
  200. call MKADIM(da,XMKel)
  201. do i=1,5
  202. do j=1,5
  203. XMK(i,j)=XMKel(i,j)
  204. enddo
  205. enddo
  206.  
  207. C Calcul de l'effort de prédiction élastique en considérant le décollement élastique non-linéaire.
  208. duP(1)=0
  209. duP(2)=0
  210. duP(3)=0
  211. duP(4)=0
  212. duP(5)=0
  213.  
  214. C Réduction du pas de calcul
  215. dtimei=dtime/10.
  216. DO I=1,5
  217. duoi(I)=1./10.*duo(I)
  218. force_preci(I)=force_prec(I)
  219. ENDDO
  220. delta_preci=delta_prec
  221. Xkapa_preci=Xkapa_prec
  222. Xksi_preci=Xksi_prec
  223. distmin_preci=distmin_prec
  224. do k=1,10
  225. q2_preci=q2
  226. *
  227. Xkapa_k=Xkapa_preci
  228. Xksi_k=Xksi_preci
  229. DO I=1,5
  230. q1_preci(I)=q1(I)
  231. upl_preci(I)=upl(I)
  232. uoi(I)=uo(I)-(10.-k)*duoi(I)
  233. uel(I)=uoi(I)-upl(I)
  234. duel(I)=duoi(I)-duP(I)
  235. ENDDO
  236. *si a8=1 on désactive le décollement
  237. if (a8.eq.1) then
  238. call MATVE1 (XMK, uel,5,5,Ftr,2)
  239. delta=delta_prec
  240. deltamax=1
  241. else
  242. call MATVE1 (XMK, duel,5,5,dforce,2)
  243. DO I=1,5
  244. force(I)=force_preci(I)+dforce(I)
  245. ENDDO
  246. call ELNLIN (da, uel, force, XMKel,
  247. & delta_preci, Ftr, delta, deltamax)
  248. endif
  249. C Calcul des fonctions de charge au point de charge élastique
  250. call LOADSF (da, Ftr,q1,yield1)
  251. call LOASFV (Ftr,q2,yield2)
  252. DO I=1,5
  253. force(I)=Ftr(I)
  254. ENDDO
  255. C Boucle de plasticité (Return Mapping)
  256. epsilon=1E-9
  257. if ((yield1.le.epsilon) .and. (yield2.le.epsilon)) then
  258. DO I=1,5
  259. force(I)=Ftr(I)
  260. ENDDO
  261. Xkapa=Xkapa_k
  262. Xksi=Xksi_k
  263. distmin=distmin_preci
  264. else
  265. m=0
  266. do while ((yield1.gt.epsilon) .or. (yield2.gt.epsilon))
  267. m=m+1
  268. IF (m.gt.1000) THEN
  269. KERRE = 22
  270. RETURN
  271. ENDIF
  272. IF ((force(1).LT.0).or.(force(1).GE.1)) THEN
  273. KERRE = 23
  274. RETURN
  275. ENDIF
  276. call GRADIF(da,force,q1,Gradfs,Gradfq)
  277. call GRADIG(da,force, q1, force_preci, Xkapa_k,
  278. & Xksi_k, Xkapa_preci, Xksi_preci, Gradgs, Xkapa, Xksi)
  279. call FVEH(da,force,q1, XMKel,
  280. & Gradgs,distmin_preci,hh1,distmin)
  281. call GRADFV(force,q2,Gradf2s,Gradf2q)
  282. call FVEHV(force,q2, XMKel,Gradf2s,hh2)
  283.  
  284. call MATVE1(XMK,Gradgs,5,5,XKGradgs,2)
  285. call MATVE1(XMK,Gradf2s,5,5,XKGradf2s,2)
  286.  
  287. call SCALPR(5,Gradfs,XKGradgs,Ho11)
  288. call SCALPR(5,Gradfs,XKGradf2s,Ho12)
  289. call SCALPR(5,Gradf2s,XKGradgs,Ho21)
  290. call SCALPR(5,Gradf2s,XKGradf2s,Ho22)
  291.  
  292. call SCALPR(5,Gradfq,hh1,H1)
  293. H2=Gradf2q*hh2
  294. H(1,1)=Ho11+H1
  295. H(1,2)=Ho12
  296. H(2,1)=Ho21
  297. H(2,2)=Ho22+H2
  298. yield(1)=yield1
  299. yield(2)=yield2
  300. *initialisation de ces variables necessaire pour les tests sur les mécanismes
  301. ddlambda(1)=1
  302. ddlambda(2)=1
  303. if ((yield(1).gt.epsilon).
  304. & and.(yield(2).gt.epsilon)) then
  305. call HINV22 (H,Hinv)
  306. call MATVE1(Hinv,yield,2,2,ddlambda,2)
  307. ddlambda1=ddlambda(1)
  308. ddlambda2=ddlambda(2)
  309. endif
  310. C Tests sur les mécanismes
  311. if ((yield(1).le.epsilon).or.(ddlambda(1).le.0)) then
  312. ddlambda1=0
  313. ddlambda2=yield2/H(2,2)
  314. if (ddlambda2.le.0) then
  315. ddlambda2=0
  316. endif
  317. endif
  318. if ((yield(2).le.epsilon).or.(ddlambda(2).le.0)) then
  319. ddlambda2=0
  320. ddlambda1=yield1/H(1,1)
  321. if (ddlambda1.le.0) then
  322. ddlambda1=0
  323. endif
  324. endif
  325.  
  326. ddlambda(1)=ddlambda1
  327. ddlambda(2)=ddlambda2
  328.  
  329.  
  330. q2=q2-ddlambda(2)*hh2
  331. DO I=1,5
  332. q1(I)=q1(I)-ddlambda(1)*hh1(I)
  333. upl1(I)=upl1(I)+ddlambda(1)*Gradgs(I)
  334. upl2(I)=upl2(I)+ddlambda(2)*Gradf2s(I)
  335. upl(I) =upl1(I)+upl2(I)
  336. *
  337. uel(I)=uoi(I)-upl(I)
  338. duP(I)=upl(I)-upl_preci(I)
  339. duel(I)=duoi(I)-duP(I)
  340. ENDDO
  341. *si a8=1 on désactive le décollement
  342. if (a8.eq.1) then
  343. call MATVE1 (XMK, uel,5,5, force,2)
  344. delta=delta_prec
  345. else
  346. call MATVE1 (XMK, duel, 5,5, dforce,2)
  347. DO I=1,5
  348. forcepred(I)=force_preci(I)+dforce(I)
  349. ENDDO
  350. *forcepred et delta_prec sont là pour une initialisation à la boucle de Newton
  351. call ELNLIN (da, uel, forcepred, XMKel,
  352. & delta_preci, force, delta, deltamax)
  353. endif
  354. Xkapa_k=Xkapa
  355. Xksi_k=Xksi
  356. call LOADSF (da,force,q1,yield1)
  357. call LOASFV (force,q2,yield2)
  358. enddo
  359. endif
  360. *
  361. C Correction visqueuse
  362. * calcul de la norme de la force dans l'hyperplan (pour voir si on est dans la phase de chargement statique poids propre)
  363. Xnorm_f25=sqrt(force(2)**2+force(3)**2+force(4)**2+force(5)**2)
  364. if (dtime.ne.0) then
  365. if ((eta3.ne.0) .and. (Xnorm_f25.gt.0.001)) then
  366. DO I=1,5
  367. force(I) = (Ftr(I)+dtime/eta3*force(I))/(1+dtime/eta3)
  368. upl(I) = (upl_preci(I)+dtime/eta3*upl(I))/(1+dtime/eta3)
  369. q1(I) = (q1_preci(I)+dtime/eta3*q1(I))/(1+dtime/eta3)
  370. ENDDO
  371. q2 = (q2_preci+dtime/eta3*q2)/(1+dtime/eta3)
  372. endif
  373. endif
  374. *
  375. IF (force(1).GE.1) THEN
  376. KERRE = 25
  377. RETURN
  378. ENDIF
  379. delta_preci=delta
  380. DO I=1,5
  381. force_preci(I)=force(I)
  382. ENDDO
  383. Xkapa_preci=Xkapa
  384. Xksi_preci=Xksi
  385. distmin_preci=distmin
  386. *
  387. enddo
  388. *
  389. call SCALPR(5,upl,upl,EPSE)
  390. EPSE=sqrt(EPSE)
  391. *
  392. VARF(1)=EPSE
  393. VARF(2)=upl1(1)
  394. VARF(3)=upl1(2)
  395. VARF(4)=upl1(3)
  396. VARF(5)=upl1(4)
  397. VARF(6)=upl1(5)
  398. VARF(7)=q1(1)
  399. VARF(8)=q1(2)
  400. VARF(9)=q1(3)
  401. VARF(10)=q1(4)
  402. VARF(11)=q1(5)
  403. VARF(12)=delta
  404. VARF(13)=Xkapa
  405. VARF(14)=Xksi
  406. VARF(15)=distmin
  407. *deltamax ne sert pas, on la garde juste pour la tracer (provisoire)
  408. VARF(16)=deltamax
  409. VARF(17)=upl2(1)
  410. VARF(18)=upl2(2)
  411. VARF(19)=upl2(3)
  412. VARF(20)=upl2(4)
  413. VARF(21)=upl2(5)
  414. VARF(22)=q2
  415.  
  416. C Redimensionnement des variables
  417. if (a9.eq.1) then
  418. C semelle filante
  419. DO I=1,5
  420. force(I)=force(I)*Diam*qmax
  421. ENDDO
  422. force(3)=force(3)*Diam
  423. force(5)=force(5)*Diam
  424.  
  425. else
  426. if (a9.eq.2) then
  427. C semelle rectangulaire
  428. DO I=1,5
  429. force(I)=force(I)*XLx*XLy*qmax
  430. ENDDO
  431. force(3)=force(3)*XLx
  432. force(5)=force(5)*XLy
  433. else
  434. C semelle circulaire
  435. DO I=1,5
  436. force(I)=force(I)*XPI*Diam**2/4*qmax
  437. ENDDO
  438. force(3)=force(3)*Diam
  439. force(5)=force(5)*Diam
  440. endif
  441. endif
  442. *
  443. call MKREDI(da,XMKel)
  444. do i=1,5
  445. do j=1,5
  446. XMK(i,j)=XMKel(i,j)
  447. enddo
  448. enddo
  449. *
  450. *Fz ISS --> -Fx joint
  451. SIGF(1)=-force(1)
  452. *Fx ISS --> -Fy joint
  453. SIGF(2)=-force(2)
  454. *Fy ISS --> Fz joint
  455. SIGF(3)=force(4)
  456. *Mz ISS --> -Mx joint (torsion linéaire)
  457. SIGF(4)=-XKelt*epst0(4)
  458. *Mx ISS --> -My joint
  459. SIGF(5)=-force(5)
  460. *My ISS --> Mz joint
  461. SIGF(6)=force(3)
  462. *
  463. return
  464. end
  465.  
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  

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