Télécharger fefp2.eso

Retour à la liste

Numérotation des lignes :

fefp2
  1. C FEFP2 SOURCE CB215821 20/07/29 21:15:34 10668
  2.  
  3. *******************************************************************
  4. SUBROUTINE FEFP2(MATE,INPLAS,MELE,MELEME,MINTE,xMATRI,
  5. . NBELEM,NBPTEL,NBNN,LRE,MFR,
  6. . IVADESP,IVADPI,IVARI,IVAMAT,
  7. . IVASTF,IVARIF,IVADPF,LHOOK,IRIGE7,
  8. . NDEP,NDEF,NSTRS,NVARI,NMATT,PRECIS,NITMAX,NUPDATE,
  9. . KERRE)
  10. *******************************************************************
  11. *
  12. * entrees :
  13. * mate = numero de material elastico
  14. * inplas = numero de material inelastico
  15. * mele = numero de elemento finito
  16. * meleme = puntero a la malla
  17. * minte = puntero valores integracion
  18. * xmatri = puntero de matrices elementares
  19. * nbelem = numero de elementos
  20. * nbptel = nombre de puntos por elemento
  21. * nbnn =
  22. * lre =
  23. * mfr =
  24. * ivadesp= puntero a un segmento mptval de desplazamientos
  25. * ivadpi = puntero a un segmento mptval de deformaciones inelasticas
  26. * ivari = puntero a un segmento mptval de variables internas
  27. * ivamat = puntero a un segmento mptval de material
  28. * lhook = dimensiones de la matriz de hook
  29. * irige7 = simetria de la matriz tangente
  30. * ndep = numero de componentes de desplazamientos
  31. * ndef = numero de componentes de deformaciones inelasticas
  32. * nstrs = numero de componentes de tensiones
  33. * nvari = numero de componentes de variables internas
  34. * nmatt = numero de componentes de propiedades materiales
  35. * precis = precision de las iteraciones internas
  36. * nitmax = numero maximo de iteraciones
  37. * nupdate = total (0) / update (1) lagrangian formulation
  38. *
  39. * sorties :
  40. * ivastf = puntero a un segmento mptval de tensiones
  41. * ivarif = puntero a un segmento mptval de variables internas
  42. * ivadpf = puntero a un segmento mptval de deformaciones inelasticas
  43. * kerre = indicador de error
  44. *
  45. ************************************************************************
  46.  
  47. IMPLICIT INTEGER(I-N)
  48. IMPLICIT REAL*8(A-H,O-Z)
  49.  
  50.  
  51. -INC PPARAM
  52. -INC CCOPTIO
  53. -INC SMCHAML
  54. -INC SMELEME
  55. -INC SMCOORD
  56. -INC SMMODEL
  57. -INC SMINTE
  58. -INC SMRIGID
  59.  
  60. SEGMENT MPTVAL
  61. INTEGER IPOS(NS),NSOF(NS),IVAL(NCOSOU)
  62. CHARACTER*16 TYVAL(NCOSOU)
  63. ENDSEGMENT
  64. cc hay +4 en XMAT para compatibilidad con feap
  65. SEGMENT WRK0
  66. REAL*8 XMAT(NMATT+4)
  67. REAL*8 DDHOOK(LHOOK,LHOOK),DEFPINI(NDEF),DEFPFIN(NDEF)
  68. REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI),BE(NDEF)
  69. REAL*8 REL(LRE,LRE),REL2(LRE,LRE),SHPWRK(6,NBNN)
  70. REAL*8 BGENE(NSTRS,LRE)
  71. REAL*8 XENEW(3,NBNN),XEOLD(3,NBNN),XMOVI(3,NBNN)
  72. REAL*8 FFI(3,3),FF(3,3),FFIv(9),FFv(9)
  73. ENDSEGMENT
  74.  
  75. SEGINI WRK0
  76.  
  77. ******************
  78. nescri =0
  79. nues =199
  80. c numerical differentiation parameters
  81. nnumer =0
  82. deltax =0.D0
  83. c line-search parameters
  84. level =1
  85. kmax =3
  86. iaugla =0
  87. caugla =0.D0
  88. c number of internal variable 'relative dendity'
  89. nvardet =4
  90. c model number
  91. if (inplas.eq.114) then
  92. nmodel =8
  93. else if (inplas.eq.115) then
  94. nmodel =2
  95. else if (inplas.eq.116) then
  96. nmodel =5
  97. c DO NOT USE LINE-SEARCH
  98. kmax =0
  99. else if (inplas.eq.117) then
  100. nmodel =6
  101. else
  102. write(nues,*)' Model not defined',inplas
  103. endif
  104. c*****************************************************************
  105. c INICIO bucle elementos
  106. c*****************************************************************
  107. DO 1000 IB=1,NBELEM
  108. c*****************************************************************
  109.  
  110. CALL ZERO(REL,LRE,LRE)
  111. CALL ZERO(REL2,LRE,LRE)
  112.  
  113. c coord del elemento ib (en t_(n+1) por el form)
  114. CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XENEW)
  115.  
  116. c desplazamientos impuestos y coordenadas en t_n
  117. MPTVAL=IVADESP
  118. DO INO=1,NBNN
  119. DO ID=1,IDIM
  120. MELVAL=IVAL(ID)
  121. IBMN=MIN(IB,VELCHE(/2))
  122. INMN=MIN(INO,VELCHE(/1))
  123. XMOVI(ID,INO)=VELCHE(INMN,IBMN)
  124. XEOLD(ID,INO)=XENEW(ID,INO)-XMOVI(ID,INO)
  125. ENDDO
  126. ENDDO
  127. cc
  128. c*****************************************************************
  129. c INICIO bucle puntos de gauss
  130. c*****************************************************************
  131. DO 1100 IGAU=1,NBPTEL
  132. c*****************************************************************
  133.  
  134. * var0 = variables internas iniciales
  135. MPTVAL=IVARI
  136. DO IC=1,NVARI
  137. MELVAL=IVAL(IC)
  138. IBMN=MIN(IB,VELCHE(/2))
  139. IGMN=MIN(IGAU,VELCHE(/1))
  140. VAR0(IC)=VELCHE(IGMN,IBMN)
  141. enddo
  142.  
  143. * defplini = deformacion inelastica inicial
  144. MPTVAL=IVADPI
  145. DO IC=1,NDEF
  146. MELVAL=IVAL(IC)
  147. IBMN=MIN(IB,VELCHE(/2))
  148. IGMN=MIN(IGAU,VELCHE(/1))
  149. DEFPINI(IC)=VELCHE(IGMN,IBMN)
  150. enddo
  151.  
  152. * xmat = caracteristicas materiales
  153. MPTVAL=IVAMAT
  154. DO IC=1,NMATT-5
  155. MELVAL =IVAL(IC)
  156. IGMN =MIN(IGAU,VELCHE(/1))
  157. IBMN =MIN(IB ,VELCHE(/2))
  158. XMAT(IC+2)=VELCHE(IGMN,IBMN)
  159. ENDDO
  160. XMAT(1)=XMAT(3)
  161. XMAT(2)=XMAT(4)
  162. XMAT(3)=0.D0
  163. XMAT(4)=0.D0
  164. c compatibilidad con feap, ogden model
  165. if (nmodel.eq.8) then
  166. XMAT(10)=XMAT(9)
  167. XMAT(9) =0.D0
  168. XMAT(13)=-1.D0
  169. endif
  170.  
  171. ****************************************
  172. * integrar
  173. ****************************************
  174.  
  175. CALL HPRIME(XENEW,NBNN,IDIM,SHPTOT,IGAU,SHPWRK,DJAC)
  176.  
  177. C Matriz Finv (d x_old / d x_new)
  178. CALL ZERO(FFI,3,3)
  179. DO ID=1,NBNN
  180. DO IE=1,IDIM
  181. DO IF=1,IDIM
  182. FFI(IE,IF)=FFI(IE,IF)-SHPWRK(IF+1,ID)*XMOVI(IE,ID)
  183. ENDDO
  184. ENDDO
  185. ENDDO
  186. FFI(1,1)=1.D0+FFI(1,1)
  187. FFI(2,2)=1.D0+FFI(2,2)
  188. FFI(3,3)=1.D0+FFI(3,3)
  189. IF (IDIM.EQ.2) THEN
  190. FFI(3,3)=1.D0
  191. IF (IFOUR.EQ.0) THEN
  192. c Cas axisim
  193. R1=0.D0
  194. R2=0.D0
  195. DO ID=1,NBNN
  196. R1=R1+SHPWRK(1,ID)*XENEW(1,ID)
  197. R2=R2+SHPWRK(1,ID)*XEOLD(1,ID)
  198. ENDDO
  199. FFI(3,3)=R2/(R1+1.E-20)
  200. ENDIF
  201. ENDIF
  202.  
  203. C Matrix F (d x_new / d x_old)
  204. C* CALL ZERO(FF ,3,3)
  205. DO I=1,3
  206. DO J=1,3
  207. FF(I,J)=FFI(I,J)
  208. ENDDO
  209. ENDDO
  210. CALL INVALM(FF,3,3,KERRE,1.D-12)
  211. IF (KERRE.NE.0) THEN
  212. WRITE(*,*) ' MATRICE SINGULIERE'
  213. RETURN
  214. ENDIF
  215.  
  216. c Almacenamiento vectorial
  217. DO i = 1,3
  218. FFv(i) =FF(i,i)
  219. FFIv(i)=FFI(i,i)
  220. ENDDO
  221. FFv(4) =FF(1,2)
  222. FFv(5) =FF(2,1)
  223. FFIv(4)=FFI(1,2)
  224. FFIv(5)=FFI(2,1)
  225. IF (IDIM.EQ.3) THEN
  226. write(*,*)' 3D FeFp not available'
  227. return
  228. ENDIF
  229.  
  230. C Determinante de Finv
  231. IF (IDIM.EQ.2) THEN
  232. DETFFI=(FFI(1,1)*FFI(2,2)-FFI(1,2)*FFI(2,1))*FFI(3,3)
  233. ELSEIF (IDIM.EQ.3) THEN
  234. DETFFI= FFI(1,1)*(FFI(2,2)*FFI(3,3)-FFI(3,2)*FFI(2,3))
  235. DETFFI=DETFFI-FFI(2,1)*(FFI(1,2)*FFI(3,3)-FFI(3,2)*FFI(1,3))
  236. DETFFI=DETFFI+FFI(3,1)*(FFI(1,2)*FFI(2,3)-FFI(1,3)*FFI(2,2))
  237. ENDIF
  238.  
  239. ********************************
  240. c Initialize plastic internal variables in first time step
  241. aux = abs(DEFPINI(1))+abs(DEFPINI(2))+abs(DEFPINI(3))
  242. if (aux.lt.1.d-14) then
  243. call zero(DEFPINI,ndef,1)
  244. do i = 1,3
  245. DEFPINI(i) = 1.d0
  246. end do
  247. endif
  248.  
  249. c Inicializar la variable interna 'det(FF_tot)' si es Update_Lagrangian
  250. if ((nupdate.eq.1).AND.(VAR0(nvardet).lt.1.D-14))
  251. . VAR0(nvardet) = 1.D0
  252.  
  253. c Push-Forward plastic internal variables
  254. call pushf35(DEFPINI,FFv,BE)
  255.  
  256. c Copy internal variables to t_n+1
  257. do i = 1,NVARI
  258. VARF(i) = VAR0(i)
  259. end do
  260.  
  261. c density dependent plasticity, POWDER & PODERCAP
  262. xdensi=-1.D0
  263. if (inplas.eq.116.or.inplas.eq.117) then
  264. xdensi = VAR0(nvardet) * DETFFI
  265. endif
  266.  
  267. c Compute KIRCHHOFF stress + tangent
  268. call apf_driver_fefp(BE,VARF,SIGF,DDHOOK,
  269. . NDEF,NVARI,NSTRS,LHOOK,
  270. . XMAT,xdensi,PRECIS,NITMAX,KERRE,
  271. . nescri,nues,nmodel,nnumer,deltax,
  272. . level,kmax,iaugla,caugla,ib,igau,izone)
  273.  
  274. c Total Lagrangian
  275. if (nupdate.eq.0) then
  276. c Guardar deformaciones elasticas en referencia
  277. call pushf35(BE,FFIv,DEFPFIN)
  278.  
  279. c Update Lagrangian
  280. else
  281. c Guardar deformaciones elasticas en t_(n+1)
  282. call equv(DEFPFIN,BE,NDEF)
  283. c Actualizar el determinante de FF^(-1)_(total)
  284. VARF(nvardet) = VAR0(nvardet) * DETFFI
  285. DETFFI = VARF(nvardet)
  286. endif
  287.  
  288. c Get CAUCHY stresses and tangent moduli
  289. do i = 1,NSTRS
  290. SIGF(i) = SIGF(i)*DETFFI
  291. do j = 1,NSTRS
  292. DDHOOK(i,j) = DDHOOK(i,j)*DETFFI
  293. enddo
  294. enddo
  295.  
  296. c add pressure contribution, log model (only deviatoric models)
  297. c al ser deviatoric models => det(Fe) = det (F)
  298. c por este motivo se puede calcular con este DETFFI
  299. c los casos no isocoricos calculan a partir de BE, que solo tiene
  300. c la parte elastica
  301. if (nmodel.eq.8) then
  302. if (abs(DETFFI).lt.1.D-14) DETFFI = 1.D-14
  303. uppj = xmat(1)/3.D0/(1.D0-2.D0*xmat(2)) * DETFFI
  304. press = uppj * LOG(ABS(1.D0/DETFFI))
  305. do i = 1,3
  306. SIGF(i) = SIGF(i) + press
  307. DDHOOK(i,i) = DDHOOK(i,i) - 2.d0*press + uppj
  308. enddo
  309. do i = 4,NSTRS
  310. DDHOOK(i,i) = DDHOOK(i,i) - press
  311. enddo
  312. DDHOOK(1,2) = DDHOOK(1,2) + uppj
  313. DDHOOK(1,3) = DDHOOK(1,3) + uppj
  314. DDHOOK(2,3) = DDHOOK(2,3) + uppj
  315. endif
  316.  
  317. IF (IRIGE7.EQ.0)THEN
  318. c compute lower part by symmetry
  319. do i = 2,NSTRS
  320. do j = 1,i-1
  321. DDHOOK(i,j) = DDHOOK(j,i)
  322. enddo
  323. enddo
  324. ENDIF
  325. ****************************************
  326. * guardar resultados: Remember: Cauchy stresses
  327. ****************************************
  328.  
  329. * sigf = tensiones finales
  330. MPTVAL=IVASTF
  331. DO IC=1,NSTRS
  332. MELVAL=MPTVAL.IVAL(IC)
  333. MELVAL.VELCHE(IGAU,IB)=SIGF(IC)
  334. enddo
  335.  
  336. * varf = variables internas finales
  337. MPTVAL=IVARIF
  338. DO IC=1,NVARI
  339. MELVAL=MPTVAL.IVAL(IC)
  340. MELVAL.VELCHE(IGAU,IB)=VARF(IC)
  341. enddo
  342.  
  343. * defpfin = deformations plasticas finales
  344. MPTVAL=IVADPF
  345. DO IC=1,NDEF
  346. MELVAL=MPTVAL.IVAL(IC)
  347. MELVAL.VELCHE(IGAU,IB)=DEFPFIN(IC)
  348. enddo
  349.  
  350. ****************************************
  351. c calcula la matriz B = BGENE y el jacobiano DJAC
  352. ****************************************
  353. XDPGE = 0.D0
  354. YDPGE = 0.D0
  355. DIM3 = 1.D0
  356. CALL BMATST(IGAU,NBPTEL,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  357. . MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NIFOUR,DIM3,
  358. . XENEW,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
  359. IF (abs(DJAC).LT.1.E-17) then
  360. write(*,*)' Jacobiano cero, en elem', ib,' gauss',igau
  361. ELSEIF (DJAC.LT.0.D0) then
  362. write(*,*)' Jacobiano negativo, en elem', ib,' gauss',igau
  363. return
  364. endif
  365. DJAC=ABS(DJAC)*MINTE.POIGAU(IGAU)
  366. c Ktan = BDB
  367. IF (IRIGE7.EQ.2)THEN
  368. CALL BDBSTS(BGENE,DJAC,DDHOOK,LRE,NSTRS,REL)
  369. ELSE
  370. CALL BDBST (BGENE,DJAC,DDHOOK,LRE,NSTRS,REL)
  371. ENDIF
  372. c Ksigma
  373. DO IA=1,NBNN
  374. DO IO=1,6
  375. SHPWRK(IO,IA)=SHPTOT(IO,IA,IGAU)
  376. ENDDO
  377. ENDDO
  378. CALL DEVOLU(XENEW,SHPWRK,MFR,NBNN,IFOUR,
  379. . NIFOUR,IDIM,DIM3,RR,DJAC)
  380. DJAC=ABS(DJAC)*MINTE.POIGAU(IGAU)
  381. IF(IFOUR.EQ.1) THEN
  382. IF(NIFOUR.EQ.0) THEN
  383. CALL THSIG1(SHPWRK,DJAC,SIGF,NBNN,LRE,REL2,RR)
  384. ELSE
  385. CALL THSIG2(SHPWRK,DJAC,SIGF,NBNN,LRE,REL2,NIFOUR,RR)
  386. ENDIF
  387. ELSEIF(IFOUR.EQ.0) THEN
  388. CALL THSIG3(SHPWRK,DJAC,SIGF,NBNN,LRE,REL2,RR)
  389. ELSE
  390. CALL THSIGH(SHPWRK,DJAC,SIGF,NBNN,IDIM,LRE,REL2)
  391. ENDIF
  392. c*****************************************************************
  393. c FIN fin bucle puntos de gauss
  394. 1100 continue
  395. c*****************************************************************
  396. c guarda la matriz de rigidez elemental REL en XMATRI.RE
  397. C* NLIGRP=LRE
  398. C* NLIGRD=LRE
  399. IF (IRIGE7.EQ.2)THEN
  400. DO IA=1,LRE
  401. DO IIB=1,IA
  402. RE(IA,IIB,ib)=REL(IA,IIB)+REL2(IA,IIB)
  403. RE(IIB,IA,ib)=REL(IIB,IA)+REL2(IA,IIB)
  404. ENDDO
  405. ENDDO
  406. ELSE
  407. DO IA=1,LRE
  408. DO IIB=1,IA
  409. RE(IA,IIB,ib)=REL(IA,IIB)+REL2(IA,IIB)
  410. RE(IIB,IA,ib)=RE(IA,IIB,ib)
  411. ENDDO
  412. ENDDO
  413. ENDIF
  414. c*****************************************************************
  415. c FIN bucle elementos
  416. 1000 continue
  417. segdes xmatri
  418. c*****************************************************************
  419. c desactivar segmentos de trabajo
  420. SEGSUP WRK0
  421.  
  422. RETURN
  423. END
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  

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