Télécharger fefp2.eso

Retour à la liste

Numérotation des lignes :

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

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