OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Updated material routines from KMM, including a new one for plasticity

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@847 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2011-02-25 16:57:42 +00:00
committed by Knut Morten Okstad
parent 4270cbc306
commit 86cff9c8c8
23 changed files with 1277 additions and 311 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Apps/FiniteDefElasticity/CMakeLists.txt
View File

@@ -92,7 +92,7 @@ IF(NOT WIN32)
ENDIF(NOT WIN32)
# Enable all warnings
#SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall")
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall")
FILE(GLOB_RECURSE NonLinEl_SRCS ${PROJECT_SOURCE_DIR}/../SIMLinEl?D.C
${PROJECT_SOURCE_DIR}/*.[Cf])

27
Apps/FiniteDefElasticity/Material/amat3d.f
View File

@@ -1,5 +1,4 @@
subroutine amat3d (taup, aap, q, tau, aa,
& ipswb3, iwr, ierr)
subroutine amat3d (ipsw, iwr, taup, aap, q, tau, aa, ierr)
C
C ---------------------------------------------------------------------
C
@@ -12,7 +11,9 @@ C
C METHOD:
C
C
C ARGUMENTS INPUT:
C ARGUMENTS INPUT:
C ipsw - Print switch
C iwr - Write unit number
C taup(3) - Principal deviatoric Kirchhoff stresses
C aap(6,6) - Deviatoric tangent moduli in principal basis
C q(6,6) - Transformation matrix: [principal]-->[std] basis
@@ -34,10 +35,10 @@ C
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -57,19 +58,19 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer i, j, k, ipswb3, iwr, ierr
integer i, j, k, ipsw, iwr, ierr
C
real*8 taup(3), aap(6,6), q(6,6), tau(6), aa(6,6), vi(6)
C
include 'const.h'
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE AMAT3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
call rprin0(taup, 1, 3, 'taup ', iwr)
call rprin0(aap , 6, 6, 'aap ', iwr)
@@ -142,9 +143,9 @@ C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE AMAT3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
call rprin0(tau, 1, 3, 'tau ', iwr)
call rprin0(aa , 6, 6, 'aa ', iwr)

23
Apps/FiniteDefElasticity/Material/cons2d.f
View File

@@ -1,12 +1,13 @@
subroutine cons2d (mTYP,mVER,detF,F,pmat,U,Sig,Cst,ipsw,iwr,ierr)
subroutine cons2d (ipsw,iter,iwr,lfirst,mTYP,mVER,
& detF,F,Fp,pmat,HV,U,Sig,Cst,ierr)
C
C CONS2D: 2D wrapper on the 3D material of FENRIS.
C
implicit none
C
integer mTYP, mVER, ipsw, iwr, ierr
real*8 detF, U, F(4), pmat(*), Sig(3), Cst(3,3)
real*8 F3D(3,3), S3D(6), C3D(6,6)
integer ipsw, iter, iwr, lfirst, mTYP, mVER, ierr
real*8 detF, F(4), Fp(4), pmat(*), HV, U, Sig(3), Cst(3,3)
real*8 F3D(3,3), F3P(3,3), S3D(6), C3D(6,6)
C
F3D = 0.0D0
F3D(1,1) = F(1)
@@ -14,9 +15,16 @@ C
F3D(1,2) = F(3)
F3D(2,2) = F(4)
F3D(3,3) = 1.0D0
call cons3d (mTYP,mVER,detF,F3D,pmat,U,S3D,C3D,ipsw,iwr,ierr)
F3P = 0.0D0
F3P(1,1) = Fp(1)
F3P(2,1) = Fp(2)
F3P(1,2) = Fp(3)
F3P(2,2) = Fp(4)
F3P(3,3) = 1.0D0
C
call cons3d (ipsw,iter,iwr,lfirst,mTYP,mVER,
& detF,F3D,F3P,pmat,HV,U,S3D,C3D,ierr)
C
Cst(1:2,1:2) = C3D(1:2,1:2)
Cst(1:2,3) = C3D(1:2,4)
Cst(3,1:2) = C3D(4,1:2)
@@ -26,4 +34,3 @@ C
C
return
end

113
Apps/FiniteDefElasticity/Material/cons3d.f
View File

@@ -1,7 +1,7 @@
subroutine cons3d (mTYP, mVER, detF, F, pmat, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
subroutine cons3d (ipsw, iter, iwr, lfirst, mTYP, mVER, nHV, nTM,
& detF, F, pMAT, HV, Engy, Sig, Cst, ierr)
C
C ---------------------------------------------------------------------
C ----------------------------------------------------------------------
C
C F E N R I S ROUTINE: CONS3D NTNU/SINTEF/DNV
C
@@ -21,13 +21,24 @@ C combined isotropic and kinematic hardening,
C 3-D stress
C
C ARGUMENTS INPUT:
C mTYP - Material type number
C mVER - Material version number
C detF - Determinant of the deformation gradient
C F - Deformation gradient
C pmat(2) - Material property array
C ipsw - Print switch
C iter - Iteration number
C iwr - Write unit number
C lFIRST - Control variable
C = 0 : The number of computed equilibrium configurations
C are greater or equal to one
C = 1 : No equilibrium configuration has been computed
C mTYP - Material type number
C mVER - Material version number
C nHV - Number of history parameters/stress points
C nTM - Number og stress/strain components per point
C detF - Determinant of the deformation gradient
C F - Deformation gradient
C pMAT - Material property array
C HV - History parameters at t_n
C
C ARGUMENTS OUTPUT:
C HV - History parameters at t_n+1
C Engy - Strain energy density
C Sig(6) - Cauchy stresses (spatial stresses)
C Cst(6,6) - Spatial constitutive tensor
@@ -45,14 +56,20 @@ C Bmod - Material property (Bulk modulus)
C Emod - Material property (Youngs modulus)
C Pratio - Material property (Poissons ratio)
C Smod - Material property (Shear modulus)
C GLe(6) - Green-Lagrangse strains
C PK2(6) - 2nd Piola-Kirchhoff stresses (material stresses)
C Cmt(6,6) - Material constitutive tensor
C incomp - logical variable
C = .true. if incompressible material ): Pratio = 0
C = .true. if incompressible material ): Pratio = 0.5
C istrt - Start state:
C = 0 : Elastic
C = 1 : Last solution
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -72,26 +89,32 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer mTYP, mVER, ipswb3, iwr, ierr
real*8 detF, Engy, F(9), pmat(2), Sig(6), Cst(6,6)
real*8 Bmod, Emod, Pratio, Smod
logical incomp
integer mTYP, mVER, ipsw, iter, iwr, lfirst, nHV, nTM, ierr
real*8 detF, Engy, F(3,*), pMAT(*), Sig(*), Cst(*), HV(*)
C
include 'const.h'
real*8 Bmod, Emod, Pratio, Smod
logical incomp, state
C
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE CONS3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iter =', iter
write(iwr,9020) 'lfirst =', lfirst
write(iwr,9020) 'mTYP =', mTYP
write(iwr,9020) 'mVER =', mVER
write(iwr,9020) 'nHV =', nHV
write(iwr,9020) 'nTM =', nTM
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Pmat =', pmat
write(iwr,9030) 'pMAT =', pMAT(1:2)
call rprin0(F , 3, 3, 'F ', iwr)
call rprin0(HV , 1, nHV, 'HV ', iwr)
endif
endif
C
@@ -99,8 +122,9 @@ C
C Compute Shear and Bulk modulus
C
if (mTYP .lt. 0) then
Bmod = pmat(1)
Smod = pmat(2)
Bmod = pmat(1)
Smod = pmat(2)
incomp = Bmod .eq. zero
else
Emod = pmat(1)
Pratio = pmat(2)
@@ -118,14 +142,14 @@ C
C
C Compute Cauchy stresses and spatial constitutive tensor
C
if (abs(mTYP) .eq. 10) then
if (abs(mTYP) .eq. 10) then
C
if (mVER .eq. 0) then
if (mVER .eq. 0) then
C
C Linear elastic isotropic material
C
call liel3d(detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
call liel3d(ipsw, iwr, detF, F, Bmod, Smod, Engy, Sig, Cst,
& ierr)
if (ierr .lt. 0) go to 7000
C
elseif (mVER .gt. 0 .and. mVER .lt. 5) then
@@ -136,16 +160,16 @@ C Compute Lame's parameter for compressibel material
C
if (.not. incomp) Bmod = Bmod - two3 * Smod
C
call stnh3d(mVER, detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
call stnh3d(ipsw, iwr, mVER, detF, F, Bmod, Smod, Engy, Sig,
& Cst, ierr)
if (ierr .lt. 0) go to 7000
C
elseif (mVER .gt. 10 .and. mVER .lt. 15) then
C
C Modified hyperelastic neoHookean model
C
call monh3d(mVER-10, detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
call monh3d(ipsw, iwr, mVER-10, detF, F, Bmod, Smod, Engy,
& Sig, Cst, ierr)
if (ierr .lt. 0) go to 7000
C
elseif (mVER .eq. 20) then
@@ -153,8 +177,8 @@ C
C Isotropic elasticity formulated in terms of
C principal logarithmic stretches
C
call ieps3d(detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
call ieps3d(ipsw, iwr, detF, F, Bmod, Smod, Engy, Sig, Cst,
& ierr)
if (ierr .lt. 0) go to 7000
C
else
@@ -163,6 +187,16 @@ C Illegal material version
C
go to 7100
endif
C
elseif (abs(mTYP) .eq. 40) then
C
C Finite stretch plasticity model
C
call plas3d(ipsw, iwr, iter, lfirst, ntm, pMAT, detF, F(1,1),
& F(1,4), HV(1), HV(7), HV(8), Sig, Cst, state, ierr)
if (ierr .lt. 0) go to 7000
C
Engy = zero
C
else
C
@@ -187,14 +221,20 @@ C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE CONS3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iter =', iter
write(iwr,9020) 'lfirst =', lfirst
write(iwr,9020) 'mTYP =', mTYP
write(iwr,9020) 'mVER =', mVER
write(iwr,9020) 'nHV =', nHV
write(iwr,9020) 'nTM =', nTM
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Pmat =', pmat
write(iwr,9030) 'pMAT =', pMAT(1:2)
call rprin0(F , 3, 3, 'F ', iwr)
call rprin0(HV , 1, nHV, 'HV ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(HV , 1, nHV, 'HV ', iwr)
call rprin0(Sig , 1, 6, 'Sig ', iwr)
call rprin0(Cst , 6, 6, 'Cst ', iwr)
call flush(iwr)
@@ -204,11 +244,12 @@ C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE CONS3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(HV , 1, nHV, 'HV ', iwr)
call rprin0(Sig , 1, 6, 'Sig ', iwr)
call rprin0(Cst , 6, 6, 'Cst ', iwr)
endif

6
Apps/FiniteDefElasticity/Material/const.h
View File

@@ -1,6 +0,0 @@
double precision zero,one,two,four
double precision one2,one3,one4,one6,onep5,two3,four3
parameter ( zero=0.0D0, one=1.0D0, two=2.0D0, four=4.0D0 )
parameter ( one2=0.5D0, one3=one/3.0D0, one4=0.25D0 )
parameter ( one6=one/6.0D0, onep5=1.5D0 )
parameter ( two3=two/3.0D0, four3=four/3.0D0 )

25
Apps/FiniteDefElasticity/Material/eigs3d.f
View File

@@ -1,5 +1,4 @@
subroutine eigs3d (v, d, nrotd,
& ipswb3, iwr, ierr)
subroutine eigs3d (ipsw, iwr, v, d, nrotd, ierr)
C
C ---------------------------------------------------------------------
C
@@ -23,6 +22,8 @@ C v(3,3) - Matrix for which eigenvalues/vectors should be
C computed for
C
C ARGUMENTS OUTPUT:
C ipsw - Print switch
C iwr - Write unit number
C v(3,3) - Matrix of eigenvectors (by column)
C d(3) - Eigenvalues associated with the columns of v
C nrotd - number of rotations to diagonalize
@@ -40,10 +41,10 @@ C
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -63,20 +64,20 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer nrotd, i, its, j, k, ipswb3, iwr, ierr
integer nrotd, i, its, j, k, ipsw, iwr, ierr
C
real*8 g, h, aij, sm, thresh, t, c, s, tau,
* v(3,3), d(3), a(3), b(3), z(3)
C
include 'const.h'
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE EIGS3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
call rprin0(v, 3, 3, 'v ', iwr)
endif
@@ -236,9 +237,9 @@ C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE EIGS3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9020) 'nrotd =', nrotd
call rprin0(v, 3, 3, 'v ', iwr)

52
Apps/FiniteDefElasticity/Material/ieps3d.f
View File

@@ -1,5 +1,5 @@
subroutine ieps3d (detF, F, lambda, mu, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
subroutine ieps3d (ipsw, iwr, detF, F, lambda, mu, Engy, Sig, Cst,
& ierr)
C
C ---------------------------------------------------------------------
C
@@ -14,7 +14,9 @@ C METHOD:
C
C
C
C ARGUMENTS INPUT:
C ARGUMENTS INPUT:
C ipsw - Print switch
C iwr - Write unit number
C detF - Determinant of the deformation gradient
C F(3,3) - Deformation gradient
C lambda - Lame's constant (=K-2*G/3)
@@ -34,14 +36,13 @@ C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C b(6) - Left Cauchy-Green deformation tensor
C
C b(6) - Left Cauchy-Green deformation tensor
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -61,27 +62,27 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer ipswb3, iwr, ierr
integer ipsw, iwr, ierr
integer i, j, k, l, iVF, i3, ierrl, nrotd
C
real*8 detF, Engy, F(3,3), lambda, mu, Sig(6), Cst(6,6)
real*8 detF, Engy, lambda, mu, F(3,3), Sig(6), Cst(6,6)
real*8 aap(6,6), b(6), bpr(6), Cstv, detFi, EngU, EngW,
& epsd(3), Press, q(6,6), qen(3,3), taup(3)
C
include 'const.h'
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE IEPS3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'lambda =', lambda
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
endif
endif
C
@@ -118,7 +119,7 @@ C
C
end do ! i
C
call eigs3d(qen, bpr, nrotd, ipswb3, iwr, ierr)
call eigs3d(ipsw, iwr, qen, bpr, nrotd, ierr)
if (ierrl .lt. 0) go to 7000
C
C Compute transformation matrix from
@@ -133,11 +134,8 @@ C
l = 1 + mod(j,3)
C
q(i,j) = qen(i,j)*qen(i,j)
C
q(i+3,j) = qen(i,j)*qen(k,j)
C
q(j,i+3) = qen(j,i)*qen(j,k)*two
C
q(j+3,i+3) = qen(j,i)*qen(l,k) + qen(j,k)*qen(l,i)
C
end do ! j
@@ -147,13 +145,13 @@ C
C Compute Deviatoric Strains, deviatoric Kirchhoff
C Stresses and Material Moduli in Principal Basis
C
call wlog3d (mu, bpr, epsd, taup, aap, EngW, ipswb3, iwr, ierrl)
call wlog3d (ipsw, iwr, mu, bpr, epsd, taup, aap, EngW, ierrl)
if (ierrl .lt. 0) go to 7000
C
C Transform deviatoric Kirchhoff Stresses and Material Moduli
C to Cauchy stresses and spatial material moduli
C
call amat3d (taup, aap, q, Sig, Cst, ipswb3, iwr, ierrl)
call amat3d (ipsw, iwr, taup, aap, q, Sig, Cst, ierrl)
if (ierrl .lt. 0) go to 7000
C
C Transform deviatoric Kirchhoff Stresses and Material Moduli
@@ -167,8 +165,8 @@ C Compute volumetric stresses (pressure) and material moduli
C
iVF = 1
C
call vfnh3d (iVF, detF, detFi, lambda, EngU, Press, Cstv,
& ipswb3, iwr, ierrl)
call vfnh3d (ipsw, iwr, iVF, detF, detFi, lambda, EngU, Press,
& Cstv, ierrl)
if (ierrl .lt. 0) go to 7000
C
C Accumulate deviatoric and volumetric contribution to
@@ -214,7 +212,7 @@ C ----------
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'lambda =', lambda
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
call rprin0(F , 3, 3, 'F ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
@@ -226,9 +224,9 @@ C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE IEPS3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)

8
Apps/FiniteDefElasticity/Material/include/feninc/const.h Normal file
View File

@@ -0,0 +1,8 @@
double precision zero,one,two,three,four
double precision one2,one3,one4,one6,onep5,two3,four3
double precision sqt23
parameter ( zero=0.0D0, one=1.0D0, two=2.0D0, three=3.0D0 )
parameter ( one2=0.5D0, one3=one/three, one4=0.25D0 )
parameter ( one6=one/6.0D0, onep5=1.5D0, four=4.0D0 )
parameter ( two3=two/three, four3=four/three )
parameter ( sqt23=sqrt(two/three) )

26
Apps/FiniteDefElasticity/Material/invert.f Normal file
View File

@@ -0,0 +1,26 @@
subroutine invert (lda,n,a,ierr)
C
implicit none
C
integer, intent(in) :: lda,n
integer, intent(out) :: ierr
real*8 , intent(inout) :: a(lda,n)
C
integer, allocatable :: ipiv(:)
real*8 , allocatable :: work(:)
real*8 :: rwork
C
allocate(ipiv(n))
call dgetrf (n,n,a,lda,ipiv,ierr)
if (ierr /= 0) then
deallocate(ipiv)
return
end if
C
call dgetri (n,a,lda,ipiv,rwork,-1,ierr)
allocate(work(int(rwork)))
call dgetri (n,a,lda,ipiv,work,size(work),ierr)
deallocate(ipiv,work)
C
return
end

33
Apps/FiniteDefElasticity/Material/liel3d.f
View File

@@ -1,5 +1,5 @@
subroutine liel3d (detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
subroutine liel3d (ipsw, iwr, detF, F, Bmod, Smod, Engy, Sig, Cst,
& ierr)
C
C ---------------------------------------------------------------------
C
@@ -17,7 +17,9 @@ C Stress(3,3) = | Stress(4) Stress(2) Stress(5) |
C | Stress(6) Stress(5) Stress(3) |
C
C
C ARGUMENTS INPUT:
C ARGUMENTS INPUT:
C ipsw - Print switch
C iwr - Write unit number
C detF - Determinant of the deformation gradient
C F - Deformation gradients
C Bmod - Bulk modulus
@@ -42,10 +44,10 @@ C PK2(6) - 2nd Piola-Kirchhoff stresses (material stresses)
C Cmt(6,6) - Material constitutive tensor
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -65,22 +67,22 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer ipswb3, iwr, ierr
integer i, j, i3
integer ipsw, iwr, ierr
integer i, j, i3, iop
C
real*8 detF, Engy, F(9), Bmod, Smod, Sig(6), Cst(6,6)
real*8 c1, c2
real*8 GLe(6), PK2(6), Cmt(6,6)
C
include 'const.h'
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE LIEL3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Bmod =', Bmod
@@ -131,7 +133,8 @@ C
C Push stresses and constitutive tensor forward from
C material to spatial coordinates in current configuration
C
call push3D (detF, F, PK2, Cmt, Sig, Cst, ipswb3, iwr)
iop = 0
call push3D (ipsw, iop, iwr, detF, F, PK2, Cmt, Sig, Cst)
C
go to 8000
C
@@ -162,9 +165,9 @@ C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE LIEL3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)

24
Apps/FiniteDefElasticity/Material/mdma3d.f
View File

@@ -1,4 +1,4 @@
subroutine mdma3d (P_bar, P_mix, Sig, DD, ipswb3, iwr)
subroutine mdma3d (P_bar, P_mix, Sig, DD, ipsw, iwr)
C
C ---------------------------------------------------------------------
C
@@ -36,10 +36,10 @@ C INTERNAL VARIABLES:
C None
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -59,21 +59,21 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer ipswb3, iwr, j
integer ipsw, iwr, j
real*8 DD(7,7), P_bar, P_mix, Sig(6)
real*8 fac1, Sig_D(6)
C
include 'const.h'
include 'include/feninc/const.h'
C
C Entry section
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE MDMA3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'P_bar =', P_bar
write(iwr,9030) 'P_mix =', P_mix
call rprin0(Sig, 1, 6, 'Sig ', iwr)
call rprin0(Sig, 1, 6, 'Sig ', iwr)
call rprin0(DD , 7, 7, 'DD ', iwr)
endif
endif
@@ -115,9 +115,9 @@ C
C Closing section
C
8000 continue
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE MDMA3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
call rprin0(DD , 7, 7, 'DD ', iwr)
endif

75
Apps/FiniteDefElasticity/Material/monh3d.f
View File

@@ -1,5 +1,5 @@
subroutine monh3d (iVF, detF, F, Bmod, Smod, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
subroutine monh3d (ipsw, iwr, iVF, detF, F, Bmod, Smod, Engy,
& Sig, Cst, ierr)
C
C ---------------------------------------------------------------------
C
@@ -29,9 +29,11 @@ C = partial^2_J ( Psi )/JC
C
C
C ARGUMENTS INPUT:
C ipsw - Print switch
C iwr - Write unit number
C iVF - Volumetric function type
C detF - Determinant of the deformation gradient
C F - Deformation gradients
C F - Deformation gradient
C Bmod - Bulk modulus
C Smod - Shear modulus
C
@@ -50,14 +52,15 @@ C None
C
C INTERNAL VARIABLES:
C nSTR - Number of stress/strain components
C b(6) - Left Cauchy-Green deformation tensor
C bm(6) - Modified Left Cauchy-Green deformation tensor
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -77,26 +80,26 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer iVF, ipswb3, iwr, ierr
integer ipsw, iwr, iVF, ierr
integer i, j, ierrl
C
real*8 detF, Engy, F(3,3), Bmod, Smod, Sig(6), Cst(6,6)
real*8 c1, c2, c3, Cstv, bm(6), detFi, Press, trbm3
real*8 detF, Engy, Bmod, Smod, F(3,3), Sig(6), Cst(6,6)
real*8 c1, c2, c3, Cstv, b(6), bm(6), detFi, Press, trbm3
C
include 'const.h'
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE MONH3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Bmod =', Bmod
write(iwr,9030) 'Smod =', Smod
write(iwr,9030) 'Smod =', Smod
call rprin0(F , 3, 3, 'F ', iwr)
endif
endif
@@ -109,27 +112,19 @@ C Compute Left Cauchy-Green deformation tensor
C
C b = F * F^T
C
bm(1) = F(1,1)*F(1,1) + F(1,2)*F(1,2) + F(1,3)*F(1,3)
bm(2) = F(2,1)*F(2,1) + F(2,2)*F(2,2) + F(2,3)*F(2,3)
bm(3) = F(3,1)*F(3,1) + F(3,2)*F(3,2) + F(3,3)*F(3,3)
bm(4) = F(1,1)*F(2,1) + F(1,2)*F(2,2) + F(1,3)*F(2,3)
bm(5) = F(2,1)*F(3,1) + F(2,2)*F(3,2) + F(2,3)*F(3,3)
bm(6) = F(1,1)*F(3,1) + F(1,2)*F(3,2) + F(1,3)*F(3,3)
C
if (ipswb3 .gt. 3) then
call rprin0(bm, 1, 6,'bm ', iwr)
endif
b(1) = F(1,1)*F(1,1) + F(1,2)*F(1,2) + F(1,3)*F(1,3)
b(2) = F(2,1)*F(2,1) + F(2,2)*F(2,2) + F(2,3)*F(2,3)
b(3) = F(3,1)*F(3,1) + F(3,2)*F(3,2) + F(3,3)*F(3,3)
b(4) = F(1,1)*F(2,1) + F(1,2)*F(2,2) + F(1,3)*F(2,3)
b(5) = F(2,1)*F(3,1) + F(2,2)*F(3,2) + F(2,3)*F(3,3)
b(6) = F(1,1)*F(3,1) + F(1,2)*F(3,2) + F(1,3)*F(3,3)
C
C Compute modified Left Cauchy-Green deformation tensor
C _
C bm = b = J^(-2/3) * b
C
detFi = one / detF
bm = bm * (detFi ** two3)
C
if (ipswb3 .gt. 3) then
call rprin0(bm, 1, 6,'bm ', iwr)
endif
bm = b * (detFi ** two3)
C
C Compute one third of the trace of the modified
C Left Cauchy-Green deformation tensor
@@ -140,10 +135,6 @@ C
bm(1) = bm(1) - trbm3
bm(2) = bm(2) - trbm3
bm(3) = bm(3) - trbm3
C
if (ipswb3 .gt. 3) then
call rprin0(bm, 1, 6,'bm ', iwr)
endif
C
C Compute deviatoric part of the constitutive tensor
C
@@ -187,18 +178,18 @@ C
Sig = Sig * detFi
Cst = Cst * detFi
C
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9030) 'Press =', Press
call rprin0(Sig, 1, 6,'Sig ', iwr)
call rprin0(Cst, 6, 6,'Cst ', iwr)
endif
C
C Compute volumetric stresses (presssure) and material moduli
C Compute volumetric stresses (pressure) and material moduli
C according to the volumetric function type
C
call vfnh3d (iVF, detF, detFi, Bmod, Engy, Press, Cstv,
& ipswb3, iwr, ierrl)
if (ierrl .lt. 0) go to 7000
call vfnh3d (ipsw, iwr, iVF, detF, detFi, Bmod, Engy, Press, Cstv,
& ierrl)
if (ierrl .lt. 0) go to 7000
C
Sig(1:3) = Sig(1:3) + Press
C
@@ -240,8 +231,8 @@ C ----------
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Bmod =', Bmod
write(iwr,9030) 'Smod =', Smod
call rprin0(F , 3, 3, 'F ', iwr)
write(iwr,9030) 'Smod =', Smod
call rprin0(F , 3, 3, 'F ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
@@ -253,9 +244,9 @@ C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE MONH3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)

2
Apps/FiniteDefElasticity/Material/pcst3d.f
View File

@@ -61,7 +61,7 @@ C
integer ipswb3, iwr, j
real*8 aa(6,6), dd(7,7)
C
include 'const.h'
include 'include/feninc/const.h'
C
C
C Entry section

866
Apps/FiniteDefElasticity/Material/plas3d.f Normal file
View File

@@ -0,0 +1,866 @@
subroutine plas3d (ipsw, iwr, iter, lfirst, ntm, pMAT, detF,
& Fn1, Fn, be, Epp, Epl, Sig, Cst, state, ierr)
C
C ----------------------------------------------------------------------
C
C F E N R I S ROUTINE: PLAS3D NTNU/SINTEF/DNV
C
C PURPOSE:
C PLAS3D: Compute Cauchy stresses and spatial constitutive tensor
C for Finite Deformation Isotropic I1-J2-J3 Plasticity
C Models in Principal Logarithmic Stretches
C Incremental Lagrangian Formulation
C METHOD:
C Energy function: Compressible Neohookean model
C with J_2/3 regularization
C _ __ _ _ __
C W(J,be) = U(J) + 0.5*mu*(J^(-2/3)be:1 - 3)
C
C Volumetric function type: U_1 = K*[ 0.25*( J^2 - 1 )-0.5*ln(J)]
C U_2 = K*[ 0.5*( J - 1 )^2 ]
C U_3 = K*[ 0.5*( ln(J) )^2 ]
C U_4 = K*[ 2.0*(J - 1 - ln(J)) ]
C Free energy function:
C
C Psi = U(J) -3*alpha*J*U'(J)*[T - T_ref]
C
C up = partial_J ( Psi )
C upp = [ partial_J ( J partial_J Psi ) - up ]/J
C = partial^2_J ( Psi )/JC
C
C
C ARGUMENTS INPUT:
C ipsw - Print switch
C iwr - Write unit
C iter - Iteration number
C lFIRST - Control variable
C = 0 : The number of computed equilibrium configurations
C are greater or equal to one
C = 1 : No equilibrium configuration has been computed
C ntm - Number of stress/strain components pr point
C pMAT - Material properties:
C pMAT(1) = Emod : Elastic Bulk modulus
C pMAT(2) = Poissons ratio
C pMAT(5) = Bmod : Bulk modulus
C pMAT(6) = Smod : Shear modulus
C pMAT(7) = H_iso : Isotropic Hardening Modulus [times sqrt(2/3)]
C pMAT(8) = H_kin : Kinematic Hardening Modulus [times (2/3)]
C pMAT(9) = iYIELD (Yield Function number)
C iYIELD = 1 : von Mises
C pMAT(10) = Y0 : Yield stress (initial)
C pMAT(11) = Y_inf : Yield stress (infinity)
C pMAT(12) = Delta : Hardening exponent
C pMAT(13) = istrt : start state
C = 0 : Elastic
C = 1 : Last solution
C detF - Determinant of the deformation gradient at t_n+1
C Fn1(3,3) - Deformation gradient at t_n+1
C Fn(3,3) - Deformation gradient at t_n
C History variables:
C be(*) - Left Cauchy-Green deformation tensor at t_n
C Epp - Cumulative plastic strain at t_n
C Epl(*) - Plastic Strain for Hardening at t_n
C
C ARGUMENTS OUTPUT:
C History variables:
C be(*) - Left Cauchy-Green deformation tensor at t_n+1
C Epp - Cumulative plastic strain at t_n+1
C Epl(*) - Plastic Strain for Hardening at t_n+1
C Sig(6) - Cauchy stresses (spatial stresses)
C Cst(6,6) - Spatial constitutive tensor
C state -
C ierr - Error flag
C = 1 : Illegal number of element nodes
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C vol - volumetric strain - vol = ( eps : 1 ) / 3
C nn_t(3,3) - principal directions (by columns) tensor form
C ll2(3) - squares of principal stretches = lambda^2
C tau(3) - principal values total Kirchhoff stress
C tt(3) - principal values deviatoric Kirchhoff stress
C pp - pressure volumetric Kirchhoff stress
C dtde(3,3) - Kirchhoff stress derivative
C dtde(a,b) = [ d tau_a / d (lambda_b^2) ] * lambda_b^2
C = [ d tau_a / d ( eps_b ) ] / 2.0d0
C Cstm(6,6) - spatial elastic moduli in principal basis
C yield - yield function value
C yfunct - yield function
C flg = .false. only yield function
C flg = .true. yield function and derivatives
C
C PRINT SWITCH:
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
integer ipsw, iwr, istrt, iter, lfirst, ntm, ierr
real*8 pMAT(*), detF, Fn1(3,3), Fn(3,3),
& Epp, be(*), Epl(*), Sig(*), Cst(6,6)
logical state
C
integer i, j, a, b, rot, it, itmax, p1(6), p2(6), ierrl, iop
real*8 Bmod, Smod, F(3,3), Finv(3,3), detFr, detr
real*8 pp, TwoG, K3inv, G2inv, Ypr, YY, dyld
real*8 ll2_tr(3), ll2(3), tau(3), tt(3)
real*8 eps_tr(3), th_tr, vol_tr, ee_tr(3), alp_n(3)
real*8 eps_e(3), vol_e, ee_e(3) , alp(3), ta(3)
real*8 I1, J2, J3, f1, f2, f3, f11, f22, f33, f12, f13, f23
real*8 gam, Hk, Hkr, aatr, epse, Y0, Yinf, beta
real*8 err, dsol(7), res(7), tres(7,7), fss(3,3)
real*8 dtde(3,3), ff(6,6), Cstm(6,6), Eppn
real*8 nn_t(3,3), nn(3), aa(3), bb(3)
real*8 tolb, tolc, d_el, xx, f112, f123, yield, yfunct
logical conv
C
include 'include/feninc/const.h'
C
data p1 /1,2,3,1,2,3/
data p2 /1,2,3,2,3,1/
data tolb / 1.0d-08 /
C
C Entry section
C
ierr = 0
C
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE PLAS3D'
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iter =', iter
write(iwr,9020) 'lfirst =', lfirst
write(iwr,9020) 'ntm =', ntm
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'Epp =', Epp
call rprin0(pMAT, 1, 13, 'pMAT ', iwr)
call rprin0(Fn1 , 3, 3, 'Fn1 ', iwr)
call rprin0(Fn , 3, 3, 'Fn ', iwr)
call rprin0(be , 1, ntm, 'be ', iwr)
call rprin0(Epl , 1, 3, 'Epl ', iwr)
endif
endif
C
C Initialize Left Cauchy-Green deformation tensor
C
if (lfirst .eq. 1) then
be(1:3) = one
end if
C
Cstm = zero
C
C Fetch material parameters & CONSTANT SETUP
C
Bmod = pMAT(5)
Smod = pMAT(6)
Y0 = sqt23 * pMAT(10) ! N.B. Radius of yield = sqrt(2/3) Sig_y
Hk = two3 * pMAT(8) ! N.B. Kinematic hardening = 2/3 * H_kin
istrt = nint(pMAT(13))
C
if (Hk .gt. zero) then
Hkr = one / Hk
else
Hkr = zero
endif
C
TwoG = two * Smod
C
itmax = 50 ! Iteration maximum #
tolc = 1.d-09 ! Tolerance for convergence
C
C Check state for iterations
C
if (iter .eq. 0) then ! First iteration in step
C
if(istrt.eq.0) then ! Elastic state requested
state = .false.
dyld = zero
else ! Last state requested
state = .true.
dyld = 1.0d-08*Y0
endif
C
else ! Not first iteration in step
C
state = .true.
dyld = zero
C
endif
C
C KINEMATIC COMPUTATIONS (Incremental Lagrangian formulation)
C
C Inverse of Fn(3,3) * J_n
C
Finv(1,1) = Fn(2,2)*Fn(3,3) - Fn(2,3)*Fn(3,2)
Finv(2,1) = Fn(2,3)*Fn(3,1) - Fn(2,1)*Fn(3,3)
Finv(3,1) = Fn(2,1)*Fn(3,2) - Fn(2,2)*Fn(3,1)
C
Finv(1,2) = Fn(3,2)*Fn(1,3) - Fn(3,3)*Fn(1,2)
Finv(2,2) = Fn(3,3)*Fn(1,1) - Fn(3,1)*Fn(1,3)
Finv(3,2) = Fn(3,1)*Fn(1,2) - Fn(3,2)*Fn(1,1)
C
Finv(1,3) = Fn(1,2)*Fn(2,3) - Fn(1,3)*Fn(2,2)
Finv(2,3) = Fn(1,3)*Fn(2,1) - Fn(1,1)*Fn(2,3)
Finv(3,3) = Fn(1,1)*Fn(2,2) - Fn(1,2)*Fn(2,1)
C
detFr = one /(Fn(1,1)*Finv(1,1)+
& Fn(1,2)*Finv(2,1)+
& Fn(1,3)*Finv(3,1))
C
do j = 1,3
do i = 1,3
F(i,j) = (Fn1(i,1)*Finv(1,j)
& + Fn1(i,2)*Finv(2,j)
& + Fn1(i,3)*Finv(3,j))*detFr
end do ! i
end do ! j
C
C Compute Elastic left Cauchy-Green tensor: be = F * Cp * F^T
C in current configuration
C
iop = 1
call push3d (ipsw, iop, iwr, one, F, be, Cst, be, Cst)
C call pushr2(F,be,be,one) ! Update configuration
C
C Compute principal stretches and directions
C
nn_t(1,1) = be(1)
nn_t(2,2) = be(2)
nn_t(3,3) = be(3)
C
nn_t(1,2) = be(4)
nn_t(2,1) = be(4)
C
nn_t(2,3) = be(5)
nn_t(3,2) = be(5)
C
nn_t(1,3) = be(6)
nn_t(3,1) = be(6)
C
call eigs3d (ipsw, iwr, nn_t, ll2_tr, rot, ierrl)
if (ierrl .lt. 0) go to 7000
C call eig3(nn_t, ll2_tr, rot)
C
C COMPUTE TRIAL KIRCHHOFF STRESS ( pressure and deviator )
C
do a = 1, 3
ee_tr(a) = log( sqrt(ll2_tr(a)) ) ! log ( lambda(a)^TR )
end do ! a
C
th_tr = ee_tr(1) + ee_tr(2) + ee_tr(3)
C
do a = 1,3
ee_tr(a) = ee_tr(a) - one3*th_tr ! Trial deviators
end do ! a
C
pp = Bmod * th_tr ! Pressure: K*th_tr
Eppn = Epp
C
do a = 1, 3
tt(a) = TwoG * ee_tr(a) ! Trial deviatoric stress
tau(a) = tt(a) + pp
alp(a) = Hk * Epl(a)
alp_n(a) = alp(a)
end do ! a
C
C Deviatoric: ta = tt - alp_dev
C
aatr = (alp(1) + alp(2) + alp(3))*one3
C
do a = 1,3
ta(a) = tt(a) - alp(a) + aatr ! Trial SigMA = ss - alp
end do ! a
C
C CHECK ELASTIC / PLASTIC STEP
C
C Compute stress invariant and yield function
C
I1 = three * (pp - aatr)
J2 = ( ta(1)*ta(1) + ta(2)*ta(2) + ta(3)*ta(3) ) * one2
J3 = ( ta(1)*ta(1)*ta(1) + ta(2)*ta(2)*ta(2) +
& ta(3)*ta(3)*ta(3) ) * one3
C
gam = zero
yield = yfunct(ipsw,iwr,pMAT,
& Epp,I1,J2,J3, f1,f2,f3,
& f11,f22,f33,f12,f13,f23,Ypr,YY,gam,.false.)
C
C Check yield
C
if ( (yield+dyld .gt. zero) .and. state ) then
C
C PLASTIC STEP --> RETURN MAP
C
conv = .false.
it = 1
C
K3inv = one3 / Bmod
G2inv = one2 / Smod
C
vol_tr = th_tr * one3
vol_e = K3inv * pp
C
do a = 1, 3
eps_tr(a) = ee_tr(a) + vol_tr
ee_e(a) = G2inv * tt(a)
eps_e(a) = ee_e(a) + vol_e
end do ! a
C
epse = 1.d0/(abs(eps_tr(1)) + abs(eps_tr(2)) + abs(eps_tr(3)))
C
d_el = ( K3inv - G2inv ) * one3
C
do while ((.not.conv).and.(it.le.itmax))
yield = yfunct(ipsw,iwr,pMAT,
& Epp,I1,J2,J3, f1,f2,f3,
& f11,f22,f33,f12,f13,f23,Ypr,YY,gam,.true.)
C
xx = f1 - f3 * two3 * J2
C
do a = 1, 3
nn(a) = xx + ( f2 + f3 * ta(a) ) * ta(a)
res(a) = eps_e(a) - eps_tr(a) + gam * nn(a)
res(a+3) = (alp_n(a) - alp(a))*Hkr + gam * nn(a)
end do ! a
C
res(7) = yield
C
err = (abs(res(1)) + abs(res(2)) + abs(res(3)))*epse
& + abs(res(7))/Y0
conv = err .lt. tolc
C
C Construct local tangent matrix
C
do a = 1, 3
aa(a) = f23 * ta(a) + f13
bb(a) = ta(a) * ta(a) - two3 * J2
end do ! a
C
f112 = f11 - one3 * f2
f123 = f12 - two3 * f3
C
do a = 1, 3
C
do b = a, 3
fss(a,b) = ( f112 + f22 * ta(a) * ta(b)
& + f33 * bb(a) * bb(b)
& + f123 * ( ta(b) + ta(a) )
& + aa(a) * bb(b) + bb(a) * aa(b) ) * gam
fss(b,a) = fss(a,b)
end do ! b
C
fss(a,a) = fss(a,a) + gam * ( f2 + two * f3 * ta(a) )
C
end do ! a
C
do a = 1, 3
C
do b = 1, 3
tres(a,b) = fss(a,b) + d_el
end do ! b
C
tres(a,a) = tres(a,a) + G2inv
C
end do ! a
C
C Kinematic and Isotropic Hardening
C
if (Hk .gt. zero) then
C
do a = 1,3
C
do b = 1,3
tres(a,b+3) = -fss(a,b)
tres(a+3,b) = -fss(a,b)
tres(a+3,b+3) = fss(a,b)
end do ! b
C
tres(a+3,a+3) = tres(a+3,a+3) + one/Hk
tres(a ,7) = nn(a)
tres(a+3,7) = -nn(a)
tres(7,a ) = nn(a)
tres(7,a+3) = -nn(a)
C
end do ! a
C
tres(7,7) = -Ypr
C
call invert (7, 7, tres, ierrl)
if (ierrl .lt. 0) go to 7000
C call invert(7,7,tres,ierrl)
C
do i = 1, 7
dsol(i) = - tres(i,1)*res(1)
& - tres(i,2)*res(2)
& - tres(i,3)*res(3)
& - tres(i,4)*res(4)
& - tres(i,5)*res(5)
& - tres(i,6)*res(6)
& - tres(i,7)*res(7)
end do ! i
C
C Isotropic hardening only
C
else
C
do a = 1,3
tres(a,4) = nn(a)
tres(4,a) = nn(a)
end do ! a
C
tres(4,4) = -Ypr
C
call invert (7, 4, tres, ierrl)
if (ierrl .lt. 0) go to 7000
C call invert(4,7,tres,ierrl)
C
do i = 1, 4
dsol(i) = - tres(i,1)*res(1)
& - tres(i,2)*res(2)
& - tres(i,3)*res(3)
& - tres(i,4)*res(7)
end do ! i
C
dsol(7) = dsol(4)
dsol(4) = zero
dsol(5) = zero
dsol(6) = zero
C
endif
C
C Update Kirchhoff stress and plastic flow
C
tau(1) = tau(1) + dsol(1)
tau(2) = tau(2) + dsol(2)
tau(3) = tau(3) + dsol(3)
gam = gam + dsol(7)
C
C Update accumulated plastic strain
C
Epp = Eppn + sqt23*gam
C
C Update Back Stress
C
alp(1) = alp(1) + dsol(4)
alp(2) = alp(2) + dsol(5)
alp(3) = alp(3) + dsol(6)
C
C Update vol.-dev. Kirchhoff stress and stress invariants
C
pp = ( tau(1) + tau(2) + tau(3) ) * one3
tt(1) = tau(1) - pp
tt(2) = tau(2) - pp
tt(3) = tau(3) - pp
C
aatr = (alp(1) + alp(2) + alp(3))*one3
C
do a = 1,3
ta(a) = tt(a) - alp(a) + aatr
end do ! a
C
I1 = three * (pp - aatr)
J2 = ( ta(1)*ta(1) + ta(2)*ta(2) +
& ta(3)*ta(3) ) * one2
J3 = ( ta(1)*ta(1)*ta(1) + ta(2)*ta(2)*ta(2) +
& ta(3)*ta(3)*ta(3) ) * one3
C
C Update vol.-dev. logarithmic strain
C
vol_e = K3inv * pp
C
do a = 1, 3
ee_e(a) = G2inv * tt(a)
eps_e(a) = ee_e(a) + vol_e
end do ! a
C
it = it + 1
C
end do ! while
C
C Warning: check convergence
C
if (.not.conv .and. iter.gt.0) then
write( *,*) ' *WARNING* No convergence in PLASFD',err,tolc
write(iwr,*) ' *WARNING* No convergence in PLASFD',err,tolc
C call plstop()
endif
C
C Update elastic left Cauchy-Green tensor and plastic acc. strain
C
do a = 1, 3
ll2 (a) = exp( two * eps_e(a) )
end do ! a
C
do i = 1, 6
be(i) = ll2(1) * nn_t(p1(i),1) * nn_t(p2(i),1)
& + ll2(2) * nn_t(p1(i),2) * nn_t(p2(i),2)
& + ll2(3) * nn_t(p1(i),3) * nn_t(p2(i),3)
end do ! i
C
C Update plastic strains
C
Epl(1) = Epl(1) + gam * nn(1)
Epl(2) = Epl(2) + gam * nn(2)
Epl(3) = Epl(3) + gam * nn(3)
C
C Compute elasto-plastic tangent
C
do b = 1, 3
do a = 1, 3
dtde(a,b) = one2 * tres(a,b)
end do ! a
end do ! b
C
else
C
C ELASTIC STEP ( only tangent computation )
C
state = .false. ! Indicate elastic on return
C
dtde(1,1) = one2 * Bmod + two3 * Smod
dtde(1,2) = one2 * Bmod - one3 * Smod
dtde(1,3) = dtde(1,2)
dtde(2,1) = dtde(1,2)
dtde(2,2) = dtde(1,1)
dtde(2,3) = dtde(1,2)
dtde(3,1) = dtde(1,2)
dtde(3,2) = dtde(1,2)
dtde(3,3) = dtde(1,1)
C
end if
C
C COMPUTE CAUCHY STRESS
C
detr = one / detF
C
do i = 1, ntm
Sig(i) = (tau(1) * nn_t(p1(i),1)*nn_t(p2(i),1)
& + tau(2) * nn_t(p1(i),2)*nn_t(p2(i),2)
& + tau(3) * nn_t(p1(i),3)*nn_t(p2(i),3)) * detr
end do ! i
C
Sig( 9) = Epp ! Plastic strain for output
Sig(10) = sqrt(onep5)*(YY + yield) ! Yield stress value
C
C TANGENT TRANSFORMATION
C
C Material tangent (computation in the principal basis)
C
do a = 1, 3
C
C Upper 3x3 block of Cstm()
C
do b = 1, 3
Cstm(b,a) = two * dtde(b,a)
end do ! b
C
Cstm(a,a) = Cstm(a,a) - two * tau(a)
C
C Lower 3x3 block of Cstm() [ diagonal block ]
C
b = mod(a,3) + 1
C
if (abs(ll2_tr(a)-ll2_tr(b)).gt.tolb) then
Cstm(a+3,a+3) = ( ll2_tr(a)*tau(b) - ll2_tr(b)*tau(a) ) /
& ( ll2_tr(b) - ll2_tr(a))
else
Cstm(a+3,a+3) = dtde(a,a) - dtde(b,a) - tau(a)
end if
C
end do ! a
C
C Transform matrix to standard basis
C
iop = 2
call push3d (ipsw, iop, iwr, detF, nn_t, Sig, Cstm, Sig, Cst)
C call tranr4(nn_t,nn_t,ff)
C call pushr4(ff,ff,dd_l,dd,detf)
go to 8000
C
C Error section
C
7000 continue
ierr =-1
go to 7900
7100 continue
ierr =-2
go to 7900
C ----------
7900 continue
write(iwr,9020) '*** ERROR IN SUBROUTINE PLAS3D *** IERR = ', ierr
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iter =', iter
write(iwr,9020) 'lfirst =', lfirst
write(iwr,9020) 'ntm =', ntm
write(iwr,9030) 'detF =', detF
call rprin0(pMAT, 1, 13, 'pMAT ', iwr)
call rprin0(Fn1 , 3, 3, 'Fn1 ', iwr)
call rprin0(Fn , 3, 3, 'Fn ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Epp =', Epp
call rprin0(be , 1, ntm, 'be ', iwr)
call rprin0(Epl , 1, 3, 'Epl ', iwr)
call rprin0(Sig , 1, 10, 'Sig ', iwr)
call rprin0(Cst , 6, 6, 'Cst ', iwr)
go to 8010
C
C Closing section
C
8000 continue
C
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE PLAS3D'
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Epp =', Epp
call rprin0(be , 1, ntm, 'be ', iwr)
call rprin0(Epl , 1, 3, 'Epl ', iwr)
call rprin0(Sig , 1, 10, 'Sig ', iwr)
call rprin0(Cst , 6, 6, 'Cst ', iwr)
endif
endif
C
8010 continue
C
return
C
9010 format(3X,A)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end
function yfunct(ipsw, iwr, pMAT, Epp, I1, J2, J3, f1, f2, f3,
& f11, f22, f33, f12, f13, f23, Ypr, YY, gam, flg)
C
C ---------------------------------------------------------------------
C
C F E N R I S ROUTINE: YFUNCT NTNU/SINTEF/DNV
C
C PURPOSE:
C YFUNCT: Compute yield function value
C
C METHOD:
C
C
C ARGUMENTS INPUT:
C ipsw - Print switch
C iwr - Write unit
C pMAT - Material properties:
C pMAT(1) = Emod : Elastic Bulk modulus
C pMAT(2) = Poissons ratio
C pMAT(5) = Bmod : Bulk modulus
C pMAT(6) = Smod : Shear modulus
C pMAT(7) = H_iso : Isotropic Hardening Modulus [times sqrt(2/3)]
C pMAT(8) = H_kin : Kinematic Hardening Modulus [times (2/3)]
C pMAT(9) = iYIELD (Yield Function number)
C iYIELD = 1 : von Mises
C pMAT(10) = Y0 : Yield stress (initial)
C pMAT(11) = Y_inf : Yield stress (infinity)
C pMAT(12) = Delta : Hardening exponent
C pMAT(13) = istrt : start state
C Epp - Cumulative plastic strain
C
C ARGUMENTS OUTPUT:
C Yield - Yield function value
C
C COMMON BLOCKS:
C const
C
C PERIPHERAL UNITS:
C None
C
C INTERNAL VARIABLES:
C vol - volumetric strain - vol = ( eps : 1 ) / 3
C nn_t(3,3) - principal directions (by columns) tensor form
C ll2(3) - squares of principal stretches = lambda^2
C tau(3) - principal values total Kirchhoff stress
C tt(3) - principal values deviatoric Kirchhoff stress
C pp - pressure volumetric Kirchhoff stress
C dtde(3,3) - Kirchhoff stress derivative
C dtde(a,b) = [ d tau_a / d (lambda_b^2) ] * lambda_b^2
C = [ d tau_a / d ( eps_b ) ] / 2.0d0
C Cstm(6,6) - spatial elastic moduli in principal basis
C yield - yield function value
C yfunct - yield function
C flg = .false. only yield function
C flg = .true. yield function and derivatives
C
C PRINT SWITCH:
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
C CODING:
C STANDART ADHERED TO: Computas report no 78-949
C AUTHOR: Kjell Magne Mathisen
C NTNU, Department of structural engineering
C September 14, 2010
C
C INSPECTED:
C
C REVISIONS:
C xx-OCT-2010 K.M.Mathisen NTNU
C Changed XG1-XG8 to XG and UG1-UG8 to UG.
C
C ---------------------------------------------------------------------
C
implicit none
C
include 'include/feninc/const.h'
C
integer ipsw, iwr, iYIELD
real*8 yfunct, yield, gam
C
real*8 pMAT(*), Epp
real*8 I1, J2, J3
real*8 f1, f2, f3, f11, f22, f33, f12, f13, f23, Ypr
real*8 J2_12, J2_1, J2_2
real*8 YY, aa, bb, ss, c1, c2
real*8 beta, Hi, Y0, Yinf
logical flg
C
f1 = zero
f2 = zero
f3 = zero
f11 = zero
f22 = zero
f33 = zero
f12 = zero
f13 = zero
f23 = zero
C
iYIELD = nint(pMAT(9)) ! Type of yield function
Y0 = sqt23 * pMAT(10) ! Radius of yield = sqrt(2/3) Sig_y
Yinf = sqt23 * pMAT(11) ! Radius of yield = sqrt(2/3) Sig_inf
beta = pMAT(12)
Hi = sqt23 * pMAT(7) ! Isotropic hardening = sqrt(2/3)*H_iso
YY = Y0 + Hi*Epp ! Sig_y(t_n) yield stress
Ypr = two3 * pMAT(7) ! Isotropic yield = 2/3 H_iso
yield =-YY ! Default return = no yield
C
C VON MISES YIELD FUNCTION
C
if (iYIELD .eq. 1) then
C
aa = (Y0 - Yinf)*exp(-beta*Epp)
YY = Yinf + aa + Hi*Epp
Ypr = Ypr - sqt23*beta*aa
ss = sqrt( two * J2 )
yield = ss - YY
C
C Compute yield function derivatives
C
if (flg) then
C
if (ss .ne. zero) then
f2 = one / ss
f22 = - f2**3
end if
C
end if
C
C DRUCKER PRAGER YIELD FUNCTION
C
else if (iYIELD .eq. 2) then
C
aa = Yinf
ss = sqrt( two * J2 )
yield = ss + one3 * aa * I1 - YY
C
C Compute yield function derivatives
C
if (flg) then
C
f1 = one3 * aa
C
if (ss .ne. zero) then
f2 = one / ss
f22 = - f2**3
end if
C
end if
C
C PRAGER-LODE YIELD FUNCTION
C
else if (iYIELD .eq. 3) then
C
if (J2 .gt. zero) then
C
if (flg) then
C
C Read in material parameters
C
bb = sqt23*pMAT(11)! Radius of yield = sqrt(2/3) Sig_inf
C
C Compute yield function
C
c1 = one / sqrt(two)
c2 = sqrt(13.5d0) * bb
C
C yield = J2 * ( one + bb * J3 * J2 **(-onep5)) - YY
yield = ( sqrt(two*J2) + c2 * J3 / J2 ) - YY
C
C Compute yield function derivatives
C
J2_1 = one / J2
J2_2 = J2_1**2
J2_12 = sqrt(J2_1)
C
f2 = c1 * J2_12 - c2 * J3 * J2_2
f3 = c2 * J2_1
f22 = - c1 * J2_12**3 * one2 + c2 * J3 * J2_1**3 * two
f23 = - c2 * J2_2
C
end if
C
else
C
yield = -YY
C
end if
C
else
C
C NO YIELD FUNCTION SPECIFIED
C
yield = zero
write (iwr,9000)
C
end if
C
yfunct = yield
C
C Format
C
9000 format(' *** NO YIELD FUNCTION SPECIFIED ***')
C
end

156
Apps/FiniteDefElasticity/Material/push3d.f
View File

@@ -1,4 +1,4 @@
subroutine push3d (detF, F, PK2, Cmt, Sig, Cst, ipswb3, iwr)
subroutine push3d (ipsw, iop, iwr, detF, F, PK2, Cmt, Sig, Cst)
C
C ---------------------------------------------------------------------
C
@@ -16,6 +16,12 @@ C Stress(3,3) = | Stress(4) Stress(2) Stress(5) |
C | Stress(6) Stress(5) Stress(3) |
C
C ARGUMENTS INPUT:
C ipsw - Prints switch
C iop - Optional parameter
C = 0 : Push stresses and constitutive tensor
C = 1 : Push stresses only
C = 2 : Push constitutive tensor only
C iwr - Write unit number
C detF - Determinant of deformation gradient
C F(3,3) - Deformation gradient
C PK2(6) - 2nd Piola-Kirchhoff stresses (material stresses)
@@ -35,10 +41,10 @@ C INTERNAL VARIABLES:
C Tmat(6,6) - Transformation matrix
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -58,28 +64,64 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer ipswb3, iwr
integer iop, ipsw, iwr
integer i, j, i1(6), i2(6)
real*8 detF, F(3,3), PK2(6), Cmt(6,6), Sig(6), Cst(6,6)
real*8 detFi, Tmat(6,6)
C
include 'const.h'
include 'include/feninc/const.h'
C
data i1 /1,2,3,1,2,3/
data i2 /1,2,3,2,3,1/
C
C Entry section
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE PUSH3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'detF =', detF
write(iwr,9020) 'iop =', iop
write(iwr,9030) 'detF =', detF
call rprin0(F , 3, 3,'F ',iwr)
call rprin0(PK2, 1, 6,'PK2 ',iwr)
call rprin0(Cmt, 6, 6,'Cmt ',iwr)
endif
endif
if (iop .eq. 0 .or. iop .eq. 1) then
call rprin0(PK2, 1, 6,'PK2 ',iwr)
end if
if (iop .eq. 0 .or. iop .eq. 2) then
call rprin0(Cmt, 6, 6,'Cmt ',iwr)
end if
end if
end if
C
detFi = one / detF
C
if (iop .eq. 0 .or. iop .eq. 1) then
C
C Push forward material stress (PK2) to current configuration
C
C Sig(i,j) = F(i,k)*PK2(k,l)*F(j,l)/detf
C
C Tmat = F^T * PK2
C
do i = 1,3
Tmat(i,1) = F(i,1)*PK2(1) + F(i,2)*PK2(4) + F(i,3)*PK2(6)
Tmat(i,2) = F(i,1)*PK2(4) + F(i,2)*PK2(2) + F(i,3)*PK2(5)
Tmat(i,3) = F(i,1)*PK2(6) + F(i,2)*PK2(5) + F(i,3)*PK2(3)
end do ! i
C
C Sig = Tmat*F
C
Sig(1) = Tmat(1,1)*F(1,1) + Tmat(1,2)*F(1,2) + Tmat(1,3)*F(1,3)
Sig(2) = Tmat(2,1)*F(2,1) + Tmat(2,2)*F(2,2) + Tmat(2,3)*F(2,3)
Sig(3) = Tmat(3,1)*F(3,1) + Tmat(3,2)*F(3,2) + Tmat(3,3)*F(3,3)
Sig(4) = Tmat(1,1)*F(2,1) + Tmat(1,2)*F(2,2) + Tmat(1,3)*F(2,3)
Sig(5) = Tmat(2,1)*F(3,1) + Tmat(2,2)*F(3,2) + Tmat(2,3)*F(3,3)
Sig(6) = Tmat(3,1)*F(1,1) + Tmat(3,2)*F(1,2) + Tmat(3,3)*F(1,3)
C
Sig = detFi * Sig
C
end if
C
if (iop .eq. 0 .or. iop .eq. 2) then
C
C Form transformation array for a 4th rank tensor in matrix form
C
@@ -89,74 +131,56 @@ C ------+-----------------------------
C (I,J) | 1,1 2,2 3,3 1,2 2,3 3,1
C or (i,j) | 2,1 3,2 1,3
C
do i = 1,3
do j = 1,3
Tmat(i,j) = F(i1(j),i1(i))*F(i2(j),i2(i))
end do ! j
do j = 4,6
Tmat(i,j) = (F(i1(j),i1(i))*F(i2(j),i2(i))
& + F(i2(j),i2(i))*F(i1(j),i1(i)))*one2
end do ! j
end do ! i
do i = 4,6
do j = 1,3
Tmat(i,j) = F(i1(j),i1(i))*F(i2(j),i2(i))
& + F(i2(j),i2(i))*F(i1(j),i1(i))
end do ! j
do j = 4,6
Tmat(i,j) = (F(i1(j),i1(i))*F(i2(j),i2(i))
& + F(i2(j),i1(i))*F(i1(j),i2(i))
& + F(i1(j),i2(i))*F(i2(j),i1(i))
& + F(i2(j),i2(i))*F(i1(j),i1(i)))*one2
end do ! j
end do ! i
do i = 1,3
do j = 1,3
Tmat(i,j) = F(i1(j),i1(i))*F(i2(j),i2(i))
end do ! j
do j = 4,6
Tmat(i,j) = (F(i1(j),i1(i))*F(i2(j),i2(i))
& + F(i2(j),i2(i))*F(i1(j),i1(i)))*one2
end do ! j
end do ! i
C
do i = 4,6
do j = 1,3
Tmat(i,j) = F(i1(j),i1(i))*F(i2(j),i2(i))
& + F(i2(j),i2(i))*F(i1(j),i1(i))
end do ! j
do j = 4,6
Tmat(i,j) = (F(i1(j),i1(i))*F(i2(j),i2(i))
& + F(i2(j),i1(i))*F(i1(j),i2(i))
& + F(i1(j),i2(i))*F(i2(j),i1(i))
& + F(i2(j),i2(i))*F(i1(j),i1(i)))*one2
end do ! j
end do ! i
C
C Compute matrix product: Cst = Tmat^t * Cmt * Tmat
C
Cst = matmul(transpose(Tmat),matmul(Cmt,Tmat))
Cst = matmul(transpose(Tmat),matmul(Cmt,Tmat))
C
C Push forward material stress (PK2) to current configuration
Cst = detFi * Cst
C
C Sig(i,j) = F(i,k)*PK2(k,l)*F(j,l)/detf
C
C
C Tmat = F^T * PK2
C
do i = 1,3
Tmat(i,1) = F(i,1)*PK2(1) + F(i,2)*PK2(4) + F(i,3)*PK2(6)
Tmat(i,2) = F(i,1)*PK2(4) + F(i,2)*PK2(2) + F(i,3)*PK2(5)
Tmat(i,3) = F(i,1)*PK2(6) + F(i,2)*PK2(5) + F(i,3)*PK2(3)
end do ! i
C
C Sig = Tmat*F
C
Sig(1) = Tmat(1,1)*F(1,1) + Tmat(1,2)*F(1,2) + Tmat(1,3)*F(1,3)
Sig(2) = Tmat(2,1)*F(2,1) + Tmat(2,2)*F(2,2) + Tmat(2,3)*F(2,3)
Sig(3) = Tmat(3,1)*F(3,1) + Tmat(3,2)*F(3,2) + Tmat(3,3)*F(3,3)
Sig(4) = Tmat(1,1)*F(2,1) + Tmat(1,2)*F(2,2) + Tmat(1,3)*F(2,3)
Sig(5) = Tmat(2,1)*F(3,1) + Tmat(2,2)*F(3,2) + Tmat(2,3)*F(3,3)
Sig(6) = Tmat(3,1)*F(1,1) + Tmat(3,2)*F(1,2) + Tmat(3,3)*F(1,3)
C
detFi = one / detf
Cst = detFi * Cst
Sig = detFi * Sig
end if
C
C Closing section
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE PUSH3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
call rprin0(Sig, 1, 6,'Sig ',iwr)
call rprin0(Cst, 6, 6,'Cst ',iwr)
if (iop .eq. 0 .or. iop .eq. 1) then
call rprin0(Sig, 1, 6,'Sig ',iwr)
end if
if (iop .eq. 0 .or. iop .eq. 2) then
call rprin0(Cmt, 6, 6,'Cmt ',iwr)
end if
endif
endif
C
return
C
9010 format(3X,A)
9020 format(3X,A,E15.6)
9020 format(3X,A,I6)
9030 format(3X,A,E15.6)
C
end

48
Apps/FiniteDefElasticity/Material/stnh3d.f
View File

@@ -1,5 +1,5 @@
subroutine stnh3d (iVF, detF, F, lambda, mu, Engy, Sig, Cst,
& ipswb3, iwr, ierr)
subroutine stnh3d (ipsw, iwr, iVF, detF, F, lambda, mu, Engy, Sig,
& Cst, ierr)
C
C ---------------------------------------------------------------------
C
@@ -31,9 +31,11 @@ C = partial^2_J ( Psi )/JC
C
C
C ARGUMENTS INPUT:
C ipsw - Print switch
C iwr - Write unit number
C iVF - Volumetric function type
C detF - Determinant of the deformation gradient
C F - Deformation gradients
C F - Deformation gradient
C lambda - Lame's constant (=K-2*G/3)
C mu - Shear modulus
C
@@ -52,13 +54,12 @@ C None
C
C INTERNAL VARIABLES:
C b(6) - Left Cauchy-Green deformation tensor
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -78,27 +79,27 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer iVF, ipswb3, iwr, ierr
integer iVF, ipsw, iwr, ierr
integer i, i3, ierrl
C
real*8 detF, Engy, F(3,3), lambda, mu, Sig(6), Cst(6,6)
real*8 c1, c2, b(6), Cstv, detFi, Press
real*8 detF, Engy, lambda, mu, b(6), F(3,3), Sig(6), Cst(6,6)
real*8 c1, c2, Cstv, detFi, Press
C
include 'const.h'
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE STNH3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'lambda =', lambda
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
endif
endif
C
@@ -118,9 +119,9 @@ C according to the volumetric function type
C
detFi = one / detF
C
call vfnh3d (iVF, detF, detFi, lambda, Engy, Press, Cstv,
& ipswb3, iwr, ierrl)
if (ierrl .lt. 0) go to 7000
call vfnh3d (ipsw, iwr, iVF, detF, detFi, lambda, Engy, Press,
& Cstv, ierrl)
if (ierrl .lt. 0) go to 7000
C
C Initialize spatial constitutive tensor
C
@@ -130,7 +131,6 @@ C Accumulate deviatoric and volumetric contribution to
C spatial stresses and material moduli
C
c1 = mu * detFi
C c2 = two * mu
c2 = two * c1
Cstv = Cstv * detF
C
@@ -176,8 +176,8 @@ C ----------
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
write(iwr,9030) 'lambda =', lambda
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
write(iwr,9030) 'mu =', mu
call rprin0(F , 3, 3, 'F ', iwr)
C
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
@@ -189,9 +189,9 @@ C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE STNH3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'Engy =', Engy
call rprin0(Sig, 1, 6,'Sig ', iwr)

34
Apps/FiniteDefElasticity/Material/vfnh3d.f
View File

@@ -1,5 +1,5 @@
subroutine vfnh3d (iVF, detF, detFi, lambda, U, Up, Upp,
& ipswb3, iwr, ierr)
subroutine vfnh3d (ipsw, iwr, iVF, detF, detFi, lambda, U, Up,
& Upp, ierr)
C
C ---------------------------------------------------------------------
C
@@ -22,6 +22,8 @@ C = partial^2_J ( Psi )/JC
C
C
C ARGUMENTS INPUT:
C ipsw - Prints switch
C iwr - Write unit number
C iVF - Volumetric function type
C detF - Determinant of the deformation gradient
C detFi - Inverse of the determinant of the Jacobian gradient
@@ -45,10 +47,10 @@ C U = U(J)
C Up = partial_J ( Psi )
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -68,20 +70,20 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer iVF, ipswb3, iwr, ierr
integer ipsw, iwr, iVF, ierr
C
real*8 detF, detFi, lambda,
& detF2, detFi2, detFm1, U, Up, Upp
C
include 'const.h'
include 'include/feninc/const.h'
C
C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE VFNH3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9020) 'iVF =', iVF
write(iwr,9030) 'detF =', detF
@@ -93,7 +95,7 @@ C
C Compute volumetric stresses (pressure) and material moduli
C according to the volumetric function type
C
if (iVF . eq. 1) then
if (iVF . eq. 1) then
C
C Volumetric function type 1:
C U(J) = lambda*0.25*(J^2 - 1 - 2*(log J))
@@ -104,7 +106,7 @@ C
Up = one2 * lambda * ( detF - detFi )
Upp = one2 * lambda * ( one + detFi2 )
C
else if (iVF .eq. 2) then
else if (iVF .eq. 2) then
C
C Volumetric function type 2:
C U(J) = lambda*0.5*(J-1)^2
@@ -114,7 +116,7 @@ C
Up = lambda * detFm1
Upp = lambda
C
else if (IVF .eq. 3) then
else if (IVF .eq. 3) then
C
C Volumetric function type 3:
C U(J) = lambda*0.5*(log J)^2
@@ -123,7 +125,7 @@ C
Up = lambda * log(abs(detF)) * detFi
Upp = ( lambda * detFi - Up ) * detFi
C
else if (iVF .eq. 4) then
else if (iVF .eq. 4) then
C
C Volumetric function type 4:
C U(J) = lambda*2.0*(J - 1 - log J)
@@ -166,9 +168,9 @@ C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE VFNH3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'U =', U
write(iwr,9030) 'Up =', Up

27
Apps/FiniteDefElasticity/Material/wlog3d.f
View File

@@ -1,5 +1,4 @@
subroutine wlog3d (mu, bpr, epsd, taup, aap, w,
& ipswb3, iwr, ierr)
subroutine wlog3d (ipsw, iwr, mu, bpr, epsd, taup, aap, w, ierr)
C
C ---------------------------------------------------------------------
C
@@ -17,7 +16,9 @@ C
C Logarithmic strain energy function
C w = mu * [log(lambda)]**2
C
C ARGUMENTS INPUT:
C ARGUMENTS INPUT:
C ipsw - Print switch
C iwr - Write unit number
C mu - Shear modulus
C bpr(3) - Principal stretch (squared)
C
@@ -39,10 +40,10 @@ C
C
C
C PRINT SWITCH:
C ipswb3 = 0 Gives no print
C ipswb3 = 2 Gives enter and leave
C ipswb3 = 3 Gives in addition parameters on input
C ipswb3 = 5 Gives in addition parameters on output
C ipsw = 0 Gives no print
C ipsw = 2 Gives enter and leave
C ipsw = 3 Gives in addition parameters on input
C ipsw = 5 Gives in addition parameters on output
C
C LIMITS:
C
@@ -62,12 +63,12 @@ C ---------------------------------------------------------------------
C
implicit none
C
integer i, j, ipswb3, iwr, ierr
integer i, j, ipsw, iwr, ierr
C
real*8 jthird, mu, twomu, fourmu, w, tol,
& bpr(3), lamt(3), epsd(3), taup(3), aap(6,6)
C
include 'const.h'
include 'include/feninc/const.h'
C
data tol / 1.0d-08 /
C
@@ -75,9 +76,9 @@ C Entry section
C
ierr = 0
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'ENTERING SUBROUTINE WLOG3D'
if (ipswb3 .gt. 2) then
if (ipsw .gt. 2) then
write(iwr,9010) 'WITH INPUT ARGUMENTS'
write(iwr,9030) 'mu =', mu
call rprin0(bpr, 1, 3, 'bpr ', iwr)
@@ -162,9 +163,9 @@ C Closing section
C
8000 continue
C
if (ipswb3 .gt. 0) then
if (ipsw .gt. 0) then
write(iwr,9010) 'LEAVING SUBROUTINE WLOG3D'
if (ipswb3 .gt. 3) then
if (ipsw .gt. 3) then
write(iwr,9010) 'WITH OUTPUT ARGUMENTS'
write(iwr,9030) 'w =', w
call rprin0(epsd, 1, 3, 'epsd ', iwr)

4
Apps/FiniteDefElasticity/NonlinearElasticity.C
View File

@@ -1,4 +1,4 @@
// $Id: NonlinearElasticity.C,v 1.3 2011-02-08 09:06:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonlinearElasticity.C
@@ -181,7 +181,7 @@ bool NonlinearElasticity::evalSol (Vector& s, const Vector&,
int ierr = 0;
if (eV && !primsol.front().empty())
if (ierr = utl::gather(MNPC,nsd,primsol.front(),*eV))
if ((ierr = utl::gather(MNPC,nsd,primsol.front(),*eV)))
{
std::cerr <<" *** NonlinearElasticity::evalSol: Detected "
<< ierr <<" node numbers out of range."<< std::endl;

25
Apps/FiniteDefElasticity/NonlinearElasticityUL.C
View File

@@ -20,15 +20,17 @@
#ifdef USE_FTNMAT
extern "C" {
//! \brief Interface to 2D nonlinear material routines (FORTRAN-77 code).
void cons2d_(const int& mTYP, const int& mVER, const double& detF,
const double* F, const double* pmat, double& Engy,
const double* Sig, double* Cst,
const int& ipsw, const int& iwr, int& ierr);
void cons2d_(const int& ipsw, const int& iter, const int& iwr,
const int& lfirst, const int& mTYP, const int& mVER,
const int& nHV, const int& nTM, const double& detF,
const double* F, const double* pMAT, double* HV, double& Engy,
const double* Sig, double* Cst, int& ierr);
//! \brief Interface to 3D nonlinear material routines (FORTRAN-77 code).
void cons3d_(const int& mTYP, const int& mVER, const double& detF,
const double* F, const double* pmat, double& Engy,
const double* Sig, double* Cst,
const int& ipsw, const int& iwr, int& ierr);
void cons3d_(const int& ipsw, const int& iter, const int& iwr,
const int& lfirst, const int& mTYP, const int& mVER,
const int& nHV, const int& nTM, const double& detF,
const double* F, const double* pMAT, double* HV, double& Engy,
const double* Sig, double* Cst, int& ierr);
}
#ifndef INT_DEBUG
#define INT_DEBUG 0
@@ -280,7 +282,6 @@ bool NonlinearElasticityUL::kinematics (const Matrix& dNdX, SymmTensor& E) const
}
const size_t nenod = dNdX.rows();
const size_t nstrc = nsd*(nsd+1)/2;
if (eV->size() != nenod*nsd || dNdX.cols() < nsd)
{
std::cerr <<" *** NonlinearElasticityUL::kinematics: Invalid dimension,"
@@ -348,9 +349,11 @@ bool NonlinearElasticityUL::constitutive (Matrix& C, SymmTensor& sigma,
int ipsw = INT_DEBUG > 1 ? 9 : 0;
int ierr = 0;
if (ndim == 2)
cons2d_(mTYP,mVER,J,F.ptr(),pmat,U,sigma.ptr(),C.ptr(),ipsw,6,ierr);
cons2d_(ipsw,0,6,0,mTYP,mVER,0,3,J,F.ptr(),pmat,
&U,U,sigma.ptr(),C.ptr(),ierr);
else
cons3d_(mTYP,mVER,J,F.ptr(),pmat,U,sigma.ptr(),C.ptr(),ipsw,6,ierr);
cons3d_(ipsw,0,6,0,mTYP,mVER,0,6,J,F.ptr(),pmat,
&U,U,sigma.ptr(),C.ptr(),ierr);
if (calcStress == 2)
{
// Transform to 2nd Piola-Kirchhoff stresses,

2
Apps/FiniteDefElasticity/NonlinearElasticityULMX.C
View File

@@ -135,7 +135,7 @@ bool NonlinearElasticityULMX::initElement (const std::vector<int>& MNPC,
// Extract the element level displacements of previous increment
int ierr = 0;
if (!primsol[1].empty())
if (ierr = utl::gather(MNPC,nsd,primsol[1],myMats->b[2]))
if ((ierr = utl::gather(MNPC,nsd,primsol[1],myMats->b[2])))
std::cerr <<" *** NonlinearElasticityULMX::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;

6
Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C
View File

@@ -1,4 +1,4 @@
// $Id: NonlinearElasticityULMixed.C,v 1.4 2011-02-08 09:06:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonlinearElasticityULMixed.C
@@ -261,7 +261,7 @@ bool NonlinearElasticityULMixed::initElement (const std::vector<int>& MNPC1,
// Extract the element level volumetric change and pressure vectors
int ierr = 0;
if (!primsol.front().empty())
if (ierr = (utl::gather(MNPC2,0,2,primsol.front(),myMats->b[T],nsd*n1,n1) +
if ((ierr = utl::gather(MNPC2,0,2,primsol.front(),myMats->b[T],nsd*n1,n1) +
utl::gather(MNPC2,1,2,primsol.front(),myMats->b[P],nsd*n1,n1)))
std::cerr <<" *** NonlinearElasticityULMixed::initElement: Detected "
<< ierr/2 <<" node numbers out of range."<< std::endl;
@@ -285,7 +285,7 @@ bool NonlinearElasticityULMixed::initElementBou (const std::vector<int>& MNPC1,
// Extract the element level displacement vector
int ierr = 0;
if (!primsol.front().empty())
if (ierr = utl::gather(MNPC1,nsd,primsol.front(),myMats->b[U]))
if ((ierr = utl::gather(MNPC1,nsd,primsol.front(),myMats->b[U])))
std::cerr <<" *** NonlinearElasticityULMixed::initElementBou: Detected "
<< ierr <<" node numbers out of range."<< std::endl;

4
Apps/SIMLinEl3D.C
View File

@@ -306,7 +306,7 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
double p = atof(strtok(NULL," "));
std::cout <<"\tPressure on P"<< press.patch
<<" F"<< (int)press.lindx <<" direction "<< pdir <<": ";
if (cline = strtok(NULL," "))
if ((cline = strtok(NULL," ")))
myTracs[1+i] = new PressureField(utl::parseRealFunc(cline,p),pdir);
else
{
@@ -330,7 +330,7 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
double E = atof(strtok(cline," "));
double nu = atof(strtok(NULL," "));
double rho = atof(strtok(NULL," "));
while (cline = strtok(NULL," "))
while ((cline = strtok(NULL," ")))
if (!strncasecmp(cline,"ALL",3))
{
std::cout <<"\tMaterial for all patches: "