Modifications of cpfem_marc2005r3 according the ones in cpfem_marc2007r1

This commit is contained in:
Luc Hantcherli 2008-03-14 21:56:46 +00:00
parent 1267ce78b6
commit 6d721dc16c
1 changed files with 154 additions and 136 deletions

View File

@ -1,15 +1,15 @@
!********************************************************************
! Material subroutine for MSC.Marc Version 0.1
!
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf
!
! last modified: 08.11.2007
!********************************************************************
! Usage:
! - choose material as hypela2
! - set statevariable 2 to index of material
! - set statevariable 3 to index of texture
! - set statevariable 2 to index of material
! - set statevariable 3 to index of texture
! - choose output of user variables if desired
! - make sure the file "mattex.mpie" exists in the working
! directory
@ -27,108 +27,119 @@
! - creeps: timinc
!********************************************************************
!
include "prec.f90"
include "debug.f90"
include "prec.f90"
include "FEsolving.f90"
include "debug.f90"
include "math.f90"
include "IO.f90"
include "mesh.f90"
include "crystal.f90"
include "constitutive.f90"
include "CPFEM.f90"
SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,&
frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,&
include "IO.f90"
include "mesh.f90"
include "crystal.f90"
include "constitutive.f90"
include "CPFEM.f90"
!
SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,&
frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,&
nnode,jtype,lclass,ifr,ifu)
!********************************************************************
! This is the Marc material routine
!********************************************************************
!
! ************* user subroutine for defining material behavior **************
!
!
! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and
! Rotation tensors at previous and current states, the analysis can be
! computationally expensive. Please use the user subroutine -> hypela
! if these kinematic quantities are not needed in the constitutive model
!
!
! IMPORTANT NOTES :
!
! (1) F,R,U are only available for continuum and membrane elements (not for
! shells and beams).
!
! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(=
! total Lagrange with large disp) in the parameter section of input deck.
! For updated Lagrangian formulation use the -> 'Plasticity,3' card(=
! update+finite+large disp+constant d) in the parameter section of
! input deck.
!
!
! d stress strain law to be formed
! g change in stress due to temperature effects
! e total elastic strain
! de increment of strain
! s stress - should be updated by user
! t state variables (comes in at t=n, must be updated
! to have state variables at t=n+1)
! dt increment of state variables
! ngens size of stress - strain law
! n element number
! nn integration point number
! kcus(1) layer number
! kcus(2) internal layer number
! matus(1) user material identification number
! matus(2) internal material identification number
! ndi number of direct components
! nshear number of shear components
! disp incremental displacements
! dispt displacements at t=n (at assembly, lovl=4) and
! displacements at t=n+1 (at stress recovery, lovl=6)
! coord coordinates
! ncrd number of coordinates
! ndeg number of degrees of freedom
! itel dimension of F and R, either 2 or 3
! nnode number of nodes per element
! jtype element type
! lclass element class
! ifr set to 1 if R has been calculated
! ifu set to 1 if strech has been calculated
!
! at t=n :
!
! ffn deformation gradient
! frotn rotation tensor
! strechn square of principal stretch ratios, lambda(i)
! eigvn(i,j) i principal direction components for j eigenvalues
!
! at t=n+1 :
!
! ffn1 deformation gradient
! frotn1 rotation tensor
! strechn1 square of principal stretch ratios, lambda(i)
! eigvn1(i,j) i principal direction components for j eigenvalues
!
! The following operation obtains U (stretch tensor) at t=n+1 :
!
! call scla(un1,0.d0,itel,itel,1)
! do 3 k=1,3
! do 2 i=1,3
! do 1 j=1,3
! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k)
!1 continue
!2 continue
!3 continue
!
! ************* user subroutine for defining material behavior **************
!
use prec, only: pReal,pInt
!
! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and
! Rotation tensors at previous and current states, the analysis can be
! computationally expensive. Please use the user subroutine -> hypela
! if these kinematic quantities are not needed in the constitutive model
!
!
! IMPORTANT NOTES :
!
! (1) F,R,U are only available for continuum and membrane elements (not for
! shells and beams).
!
! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(=
! total Lagrange with large disp) in the parameter section of input deck.
! For updated Lagrangian formulation use the -> 'Plasticity,3' card(=
! update+finite+large disp+constant d) in the parameter section of
! input deck.
!
!
! d stress strain law to be formed
! g change in stress due to temperature effects
! e total elastic strain
! de increment of strain
! s stress - should be updated by user
! t state variables (comes in at t=n, must be updated
! to have state variables at t=n+1)
! dt increment of state variables
! ngens size of stress - strain law
! n element number
! nn integration point number
! kcus(1) layer number
! kcus(2) internal layer number
! matus(1) user material identification number
! matus(2) internal material identification number
! ndi number of direct components
! nshear number of shear components
! disp incremental displacements
! dispt displacements at t=n (at assembly, lovl=4) and
! displacements at t=n+1 (at stress recovery, lovl=6)
! coord coordinates
! ncrd number of coordinates
! ndeg number of degrees of freedom
! itel dimension of F and R, either 2 or 3
! nnode number of nodes per element
! jtype element type
! lclass element class
! ifr set to 1 if R has been calculated
! ifu set to 1 if strech has been calculated
!
! at t=n :
!
! ffn deformation gradient
! frotn rotation tensor
! strechn square of principal stretch ratios, lambda(i)
! eigvn(i,j) i principal direction components for j eigenvalues
!
! at t=n+1 :
!
! ffn1 deformation gradient
! frotn1 rotation tensor
! strechn1 square of principal stretch ratios, lambda(i)
! eigvn1(i,j) i principal direction components for j eigenvalues
!
! The following operation obtains U (stretch tensor) at t=n+1 :
!
! call scla(un1,0.d0,itel,itel,1)
! do 3 k=1,3
! do 2 i=1,3
! do 1 j=1,3
! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k)
!1 continue
!2 continue
!3 continue
!
use prec, only: pReal,pInt, ijaco
use FEsolving
use CPFEM, only: CPFEM_general
use math, only: invnrmMandel
implicit real(pReal) (a-h,o-z)
!
integer(pInt) computationMode
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2)
! Marc common blocks are in fixed format so they have to be pasted in here
! Beware of changes in newer Marc versions -- these are from 2005r3
! concom is needed for inc, subinc, ncycle, lovl
! include 'concom'
! concom is needed for inc, subinc, ncycle, lovl
! include 'concom'
common/concom/ &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
@ -144,37 +155,44 @@
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, iharmt, inc_incdat, iautspc,ibrake
! creeps is needed for timinc (time increment)
! include 'creeps'
common/creeps/ &
cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
!
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2)
!
! subroutine cpfem_general(mpie_ffn, mpie_ffn1, temperature, mpie_inc, mpie_subinc, mpie_cn,
! mpie_stress_recovery, mpie_tinc, mpie_en, mpie_in, mpie_s, mpie_d, mpie_ngens)
!********************************************************************
! This routine calculates the material behaviour
!********************************************************************
! mpie_ffn deformation gradient for t=t0
! mpie_ffn1 deformation gradient for t=t1
! temperature temperature
! mpie_inc increment number
! mpie_subinc subincrement number
! mpie_cn number of cycle
! mpie_stress_recovery indicates wether we are in stiffness assemly(lovl==4) or stress recovery(lovl==6)
! mpie_tinc time increment
! mpie_en element number
! mpie_in intergration point number
! mpie_s stress vector in Marc notation, i.e. 11 22 33 12, 23, 13
! mpie_d jacoby in Marc notation
! mpie_ngens size of stress strain law
!********************************************************************
call CPFEM_general(ffn, ffn1, t(1), inc, incsub, ncycle, stress_recovery, timinc, n(1), nn, s, d, ngens)
return
ishrink,ioffsflg,isetoff, iharmt, inc_incdat, iautspc,ibrake
! creeps is needed for timinc (time increment)
! include 'creeps'
common/marc_creeps/ &
cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
if (inc == 0) then
cycleCounter = 0
else
if (theInc /= inc .or. theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1
endif
if (theInc /= inc) outdatedByNewInc = .true.
if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle
if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect
if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute
if (computationMode == 2 .and. outdatedByNewInc) then
outdatedByNewInc = .false.
computationMode = 1 ! compute and age former results
endif
theInc = inc
theCycle = ncycle
theLovl = lovl
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! Marc: 11, 22, 33, 12, 23, 13
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)
s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens)
return
END SUBROUTINE
!
@ -199,21 +217,21 @@
! ndi (3) number of direct stress components
! nshear (3) number of shear stress components
!
!********************************************************************
use prec, only: pReal,pInt
use CPFEM, only: CPFEM_results, CPFEM_Nresults
use constitutive, only: constitutive_maxNresults
use mesh, only: mesh_FEasCP
!********************************************************************
use prec, only: pReal,pInt
use CPFEM, only: CPFEM_results, CPFEM_Nresults
use constitutive, only: constitutive_maxNresults
use mesh, only: mesh_FEasCP
implicit none
!
real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*)
real(pReal) v, t(*)
integer(pInt) m, nn, layer, ndi, nshear, jpltcd
!
! assign result variable
v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
nn, mesh_FEasCP('elem', m))
integer(pInt) m, nn, layer, ndi, nshear, jpltcd
!
! assign result variable
v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
nn, mesh_FEasCP('elem', m))
return
END SUBROUTINE
!