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