适合HCP mg合金的黄永刚Vumat的晶体塑性（之前推文对应的vumat版本，同时附上对应的材料参数）

晶体塑性有限元

2025年5月20日 11:49

subroutine vumat(nblock, ndir, nshr, nstatev, nfieldv, nprops,

* lanneal, steptime, totaltime, dt, cmname, coordmp, charlength,

* props, density, straininc, relspininc, tempold, stretchold,

* defgradold, fieldold, stressold, stateold, enerinternold,

* enerinelasold, tempnew, stretchnew, defgradnew, fieldnew,

* stressnew, statenew, enerinternnew, enerinelasnew)

include 'vaba_param.inc'

c modified from harewood vumat- jg:20/07/12

dimension stressold(nblock,ndir+nshr),

1 stressnew(nblock,ndir+nshr),

1 stateold(nblock,nstatev),statenew(nblock,nstatev),

1 straininc(nblock, ndir+nshr)

dimension slpdir(3,18),slpnor(3,18),slpdef(6,18),slpspn(3,18),

1 dspdir(3,18),dspnor(3,18),d(6,6),fslip(18),dfdxsp(18),

1 ddemsd(6,18),h(18,18),dgamma(18),dtausp(18),dgslip(18),

1 dstres(6),delats(6),dvgrad(3,3),dspin(3),workst(18,18),

1 indx(18),term(3,3),trm0(3,3),props(nprops),itrm(3),

1 rotate(3,3),rwkdir(3,24),rwknor(3,24),indxL(3),termd(3),

1 termn(3),tauslpA(18)

do km=1,nblock

C----- As the VUMAT passes in tensor shear strain and this subroutine

C----- uses engineering strain --> STRAININC(shr) x 2

do i=1,nshr

straininc(km,i+3)=straininc(km,i+3)*2.d0

end do

if (nshr.gt.1) then

save=straininc(km,5)

straininc(km,5)=straininc(km,6)

straininc(km,6)=save

save=stressold(km,5)

stressold(km,5)=stressold(km,6)

stressold(km,6)=save

end if

C Elastic Matrix {D}

gshear = props(1)/(2.d0*(1.d0+props(2)))

e11 = 2.d0*gshear*(1.d0-props(2))/(1.d0-2.d0*props(2))

e12 = 2.d0*gshear*props(2)/(1.d0-2.d0*props(2))

d = 0.d0

do j = 1,3

d(j,j) = e11

do i = 1,3

if(i.ne.j) d(i,j) = e12

end do

d(j+3,j+3) = gshear

end do

c Lin Elastic Response for Exp/Packager

if(totaltime==0.)then

C Calculation of Stress Inc

dstres=0.d0

do i=1,ndir+nshr

do j=1,ndir+nshr

dstres(i)=dstres(i)+d(i,j)*straininc(km,j)

end do

C Calculation of Stress New

do i=1,ndir+nshr

stressnew(km,i)=stressold(km,i)+dstres(i)

end do

cycle

endif

c------ Crystal Type:

ictype=nint(props(9))

if(ictype == 2)then

nslptl = 6

elseif(ictype == 3)then

nslptl = 12

else

nslptl = 18

endif

C----- Integration parameter: THETA

theta = 0.5d0

term = 0.d0

trm0 = 0.d0

do i=1,3

term(i,i)=2.d0

end do

call ludcmp (term, 3, 3, itrm, ddcmp)

do j=1,3

call lubksb (term, 3, 3, itrm, trm0(1,j))

end do

dspin(1)=trm0(2,1)-trm0(1,2)

dspin(2)=trm0(1,3)-trm0(3,1)

dspin(3)=trm0(3,2)-trm0(2,3)

C----- Increment of dilatational strain: DEV

dev=0.d0

do i=1,ndir

dev=dev+straininc(km,i)

end do

C----- Check whether the current stress state is the initial state

if (stateold(km,1).eq.0.) then

c ##### basal and prismatic#####

c----- generating slip directions and normals for hcp-basal

rwkdir = 0.

rwknor = 0.

angle = acos(-1.)/3.

cangle = cos(angle)

sangle = sin(angle)

rwkdir(1,1) = 1.

rwkdir(2,1) = 0.

rwkdir(1,2) = cangle

rwkdir(2,2) = sangle

rwkdir(1,3) = -cangle

rwkdir(2,3) = sangle

rwknor(3,1) = 1.

do i = 1,3

do k = 1,3

slpdir(k,i) = rwkdir(k,i)

slpnor(k,i) = rwknor(k,1)

enddo

c----- generating slip directions and normals for hcp-prismatic

rwknor = 0.

rwknor(1,1) = 0.

rwknor(2,1) = -1.

rwknor(1,2) = sangle

rwknor(2,2) = -cangle

rwknor(1,3) = -sangle

rwknor(2,3) = -cangle

do i = 4,6

do k = 1,3

slpdir(k,i) = rwkdir(k,i-3)

slpnor(k,i) = rwknor(k,i-3)

enddo

if(ictype >= 3)then

c ##### 2nd order pyramidal <a+c> #####

aspect = 1.623

c slip directions

do j = 1,6

rwkdir(3,j) = -aspect

enddo

rwkdir(1,1) = cangle

rwkdir(2,1) = sangle

rwkdir(1,2) = -cangle

rwkdir(2,2) = sangle

rwkdir(1,3) = 2.*cangle

rwkdir(2,3) = 0.

do j = 4,6

rwkdir(1,j) = -rwkdir(1,j-3)

rwkdir(2,j) = -rwkdir(2,j-3)

enddo

rlength=sqrt(1.+aspect*aspect)

do j = 1,6

do i = 1,3

rwkdir(i,j) = rwkdir(i,j)/rlength

enddo

c slip normals

rwknor(1,1) = aspect*sangle

rwknor(2,1) = 3.*aspect*cangle

rwknor(3,1) = 4.*sangle*cangle

rwknor(1,2) = aspect*sangle

rwknor(2,2) = -3.*aspect*cangle

rwknor(3,2) = 4.*sangle*cangle

rwknor(1,3) = 2.*aspect*sangle

rwknor(2,3) = 0.

rwknor(3,3) = 4.*sangle*cangle

do j = 4,6

rwknor(1,j) = rwknor(1,j-3)

rwknor(2,j) = rwknor(2,j-3)

rwknor(3,j) = -rwknor(3,j-3)

enddo

rlength=sqrt(3.*(1.+aspect*aspect))

do j = 1,6

do i = 1,3

rwknor(i,j) = rwknor(i,j)/rlength

enddo

nslip = 6

do j = 1,6

nslip = nslip+1

do i = 1,3

slpdir(i,nslip) = rwkdir(i,j)

slpnor(i,nslip) = rwknor(i,j)

enddo

if(ictype == 4)then

c ###### twinning planes #####

c slip directions

do j = 1,6

rwkdir(3,j) = -aspect

enddo

rwkdir(1,1) = 0.

rwkdir(2,1) = 2.*sangle

rwkdir(1,2) = 3.*cangle

rwkdir(2,2) = 1.*sangle

rwkdir(1,3) = 3.*cangle

rwkdir(2,3) = -1.*sangle

do j = 4,6

rwkdir(1,j) = -rwkdir(1,j-3)

rwkdir(2,j) = -rwkdir(2,j-3)

enddo

rlength=sqrt(3.+aspect*aspect)

do j = 1,6

do i = 1,3

rwkdir(i,j) = rwkdir(i,j)/rlength

enddo

c slip normals

rwknor(1,1) = 0.

rwknor(2,1) = -2.*aspect*cangle

rwknor(3,1) = -4.*sangle*cangle

rwknor(1,2) = -aspect*sangle

rwknor(2,2) = -aspect*cangle

rwknor(3,2) = -4.*sangle*cangle

rwknor(1,3) = -aspect*sangle

rwknor(2,3) = aspect*cangle

rwknor(3,3) = -4.*sangle*cangle

do j = 4,6

rwknor(1,j) = -rwknor(1,j-3)

rwknor(2,j) = -rwknor(2,j-3)

rwknor(3,j) = rwknor(3,j-3)

enddo

do j = 1,6

do i = 1,3

rwknor(i,j) = rwknor(i,j)/rlength

enddo

do j = 1,6

nslip = nslip+1

do i = 1,3

slpdir(i,nslip) = rwkdir(i,j)

slpnor(i,nslip) = rwknor(i,j)

enddo

endif

C----- Unit vectors in slip dirs and unit norms-Global system

c----- Generate rotation matrix

do i = 1,3

term(i,1) = props(i+2)

term(i,2) = props(i+5)

enddo

term(1,3) = term(2,1)*term(3,2)-term(3,1)*term(2,2)

term(2,3) = term(3,1)*term(1,2)-term(1,1)*term(3,2)

term(3,3) = term(1,1)*term(2,2)-term(2,1)*term(1,2)

call ludcmp (term, 3, 3, indxL, dcmp)

rotate = 0.

do j = 1,3

rotate(j,j) = 1.

end do

do j = 1,3

call lubksb (term, 3, 3, indxL, rotate(1,j))

end do

c--- Rotate slip normals and directions to global system

do j = 1,nslptl

do i = 1,3

termd(i) = 0.

termn(i) = 0.

do k = 1,3

termd(i) = termd(i)+rotate(i,k)*slpdir(k,j)

termn(i) = termn(i)+rotate(i,k)*slpnor(k,j)

end do

do i = 1,3

slpdir(i,j) = termd(i)

slpnor(i,j) = termn(i)

end do

C----- Slip deformation tensor: SLPDEF (Schmid factors)

do j=1,nslptl

slpdef(1,j)=slpdir(1,j)*slpnor(1,j)

slpdef(2,j)=slpdir(2,j)*slpnor(2,j)

slpdef(3,j)=slpdir(3,j)*slpnor(3,j)

slpdef(4,j)=slpdir(1,j)*slpnor(2,j)

1 +slpdir(2,j)*slpnor(1,j)

slpdef(5,j)=slpdir(1,j)*slpnor(3,j)

1 +slpdir(3,j)*slpnor(1,j)

slpdef(6,j)=slpdir(2,j)*slpnor(3,j)

1 +slpdir(3,j)*slpnor(2,j)

end do

C----- Store normals and directions

stateold(km,nstatev) = nslptl

idnor=3*nslptl

iddir=6*nslptl

do j=1,nslptl

do i=1,3

idnor=idnor+1

stateold(km,idnor)=slpnor(i,j)

iddir=iddir+1

stateold(km,iddir)=slpdir(i,j)

end do

C----- Initial value of the current strength for all slip systems

do j=1,nslptl

if(j<=3)then

stateold(km,j)=props(10)

elseif(j<=6)then

stateold(km,j)=props(13)

elseif(j<=12)then

stateold(km,j)=props(16)

else

stateold(km,j)=props(19)

endif

enddo

C----- Initial value of shear strain in slip systems

do i=1,nslptl

stateold(km,nslptl+i)=0.d0

end do

stateold(km,9*nslptl+1)=0.d0

C----- Initial value of the resolved shear stress in slip systems

do i=1,nslptl

term1=0.

do j=1,ndir+nshr

if (j.le.ndir) then

term1=term1+slpdef(j,i)*stressold(km,j)

else

term1=term1+slpdef(j-ndir+3,i)*stressold(km,j)

end if

end do

stateold(km,2*nslptl+i)=term1

end do

else

C----- Current stress state

idnor=3*nslptl

iddir=6*nslptl

do j=1,nslptl

do i=1,3

idnor=idnor+1

slpnor(i,j)=stateold(km,idnor)

iddir=iddir+1

slpdir(i,j)=stateold(km,iddir)

end do

C----- Slip deformation tensor: SLPDEF (Schmid factors)

do j=1,nslptl

slpdef(1,j)=slpdir(1,j)*slpnor(1,j)

slpdef(2,j)=slpdir(2,j)*slpnor(2,j)

slpdef(3,j)=slpdir(3,j)*slpnor(3,j)

slpdef(4,j)=slpdir(1,j)*slpnor(2,j)

1 +slpdir(2,j)*slpnor(1,j)

slpdef(5,j)=slpdir(1,j)*slpnor(3,j)

1 +slpdir(3,j)*slpnor(1,j)

slpdef(6,j)=slpdir(2,j)*slpnor(3,j)

1 +slpdir(3,j)*slpnor(2,j)

end do

end if

C----- Slip spin tensor: SLPSPN

do j=1,nslptl

slpspn(1,j)=0.5d0*(slpdir(1,j)*slpnor(2,j)-

2 slpdir(2,j)*slpnor(1,j))

slpspn(2,j)=0.5d0*(slpdir(3,j)*slpnor(1,j)-

2 slpdir(1,j)*slpnor(3,j))

slpspn(3,j)=0.5d0*(slpdir(2,j)*slpnor(3,j)-

2 slpdir(3,j)*slpnor(2,j))

end do

C----- Double dot product of elastic moduli tensor with the slip

C deformation tensor

do j=1,nslptl

do i=1,6

ddemsd(i,j)=0.d0

do k=1,6

ddemsd(i,j)=ddemsd(i,j)+d(k,i)*slpdef(k,j)

end do

do j=1,nslptl

ddemsd(4,j)=ddemsd(4,j)-slpspn(1,j)*stressold(km,1)

ddemsd(5,j)=ddemsd(5,j)+slpspn(2,j)*stressold(km,1)

if (ndir.gt.1) then

ddemsd(4,j)=ddemsd(4,j)+slpspn(1,j)*stressold(km,2)

ddemsd(6,j)=ddemsd(6,j)-slpspn(3,j)*stressold(km,2)

end if

if (ndir.gt.2) then

ddemsd(5,j)=ddemsd(5,j)-slpspn(2,j)*stressold(km,3)

ddemsd(6,j)=ddemsd(6,j)+slpspn(3,j)*stressold(km,3)

end if

if (nshr.ge.1) then

ddemsd(1,j)=ddemsd(1,j)+slpspn(1,j)

2 *stressold(km,ndir+1)

ddemsd(2,j)=ddemsd(2,j)-slpspn(1,j)

2 *stressold(km,ndir+1)

ddemsd(5,j)=ddemsd(5,j)-slpspn(3,j)

2 *stressold(km,ndir+1)

ddemsd(6,j)=ddemsd(6,j)+slpspn(2,j)

2 *stressold(km,ndir+1)

end if

if (nshr.ge.2) then

ddemsd(1,j)=ddemsd(1,j)-slpspn(2,j)

2 *stressold(km,ndir+2)

ddemsd(3,j)=ddemsd(3,j)+slpspn(2,j)

2 *stressold(km,ndir+2)

ddemsd(4,j)=ddemsd(4,j)+slpspn(3,j)

2 *stressold(km,ndir+2)

ddemsd(6,j)=ddemsd(6,j)-slpspn(1,j)

2 *stressold(km,ndir+2)

end if

if (nshr.eq.3) then

ddemsd(2,j)=ddemsd(2,j)+slpspn(3,j)

2 *stressold(km,ndir+3)

ddemsd(3,j)=ddemsd(3,j)-slpspn(3,j)

2 *stressold(km,ndir+3)

ddemsd(4,j)=ddemsd(4,j)-slpspn(2,j)

2 *stressold(km,ndir+3)

ddemsd(5,j)=ddemsd(5,j)+slpspn(1,j)

2 *stressold(km,ndir+3)

end if

end do

C----- Shear strain-rate in a slip system at the start of increment:

do i=1,nslptl

tauslp=stateold(km,2*nslptl+i)

if(i>=13.and.tauslp<=0)then

yield=1.e6

else

yield=stateold(km,i)

endif

x=tauslp/yield

fslip(i)=0.001d0*((abs(x))**50.)*dsign(1.d0,x)

dfdxsp(i)=50.d0*0.001d0*(abs(x))**(50.-1.0)

end do

C----- Self- and latent-hardening laws

gamtol=stateold(km,9*nslptl+1)

do iself = 1,nslptl

do latent = 1,nslptl

C BASAL

if(iself <= 3)then

if(latent <= 3)then

q = 0.2

else

q = 0.5

endif

if(iself == latent)q = 1.

hlatnt = q*props(12)

C PRISM

elseif(iself <= 6)then

if(latent <= 12)then

q = 0.2

else

q = 0.5

endif

if(iself == latent)q = 1.

hlatnt = q*props(15)*(1.d0-(props(13)/props(14)))

* *exp(-props(15)*(gamtol/props(14)))

C PYRM

elseif(iself <= 12)then

if(latent <= 6)then

q = 1.

elseif(latent <= 12)then

q = 0.2

else

q = 0.25

endif

if(iself == latent)q = 1.

hlatnt = q*props(18)*(1.d0-props(16)/props(17))

* *exp(-props(18)*gamtol/props(17))

C TWIN

else

if(latent <= 6)then

q = 1.

else

q = 0.2

endif

if(iself == latent)q = 1.

if(gamtol <= props(21))then

hlatnt = q*props(20)

else

hlatnt = q*props(20)*(gamtol/props(21))

* **(props(22)-1.)

endif

h(iself,latent) = hlatnt

enddo

end do

C----- Solve the increment of shear strain in a slip system

term1=theta*dt

do i=1,nslptl

tauslp=stateold(km,2*nslptl+i)

if(i>=13.and.tauslp<=0)then

yield=1.e6

else

yield=stateold(km,i)

endif

x=tauslp/yield

term2=term1*dfdxsp(i)/yield

term3=term1*x*dfdxsp(i)/yield

do j=1,nslptl

term4=0.d0

do k=1,6

term4=term4+ddemsd(k,i)*slpdef(k,j)

end do

workst(i,j)=term2*term4+h(i,j)*term3

2 *dsign(1.d0,fslip(j))

end do

workst(i,i)=workst(i,i)+1.d0

end do

call ludcmp (workst, nslptl, 18, indx, ddcmp)

C----- Increment of shear strain in a slip system

term1=theta*dt

do i=1,nslptl

tauslp=stateold(km,2*nslptl+i)

if(i>=13.and.tauslp<=0)then

yield=1.e6

else

yield=stateold(km,i)

endif

term2=term1*dfdxsp(i)/yield

dgamma(i)=0.

do j=1,ndir

dgamma(i)=dgamma(i)+ddemsd(j,i)*straininc(km,j)

end do

if (nshr.gt.0) then

do j=1,nshr

if (j.eq.1) then

dgamma(i)=dgamma(i)+ddemsd(4,i)

2 *straininc(km,4)

elseif (j.eq.2) then

dgamma(i)=dgamma(i)+ddemsd(6,i)

2 *straininc(km,5)

elseif (j.eq.3) then

dgamma(i)=dgamma(i)+ddemsd(5,i)

2 *straininc(km,6)

end if

end do

end if

dgamma(i)=dgamma(i)*term2+fslip(i)*dt

end do

call lubksb (workst, nslptl, 18, indx, dgamma)

C----- Update the shear strain in a slip system:

do i=1,nslptl

stateold(km,nslptl+i)=stateold(km,nslptl+i)+dgamma(i)

end do

C----- Increment of current strength in a slip system: DGSLIP

do i=1,nslptl

dgslip(i)=0.

do j=1,nslptl

dgslip(i)=dgslip(i)+h(i,j)*abs(dgamma(j))

end do

C----- Update the current strength in a slip system:

do i=1,nslptl

stateold(km,i)=stateold(km,i)+dgslip(i)

end do

C----- Increment of strain associated with lattice stretching:

delats=0.

do j=1,3

if (j.le.ndir) delats(j)=straininc(km,j)

do i=1,nslptl

delats(j)=delats(j)-slpdef(j,i)*dgamma(i)

end do

do j=1,3

if (j.le.nshr) delats(j+3)=straininc(km,j+ndir)

do i=1,nslptl

delats(j+3)=delats(j+3)-slpdef(j+3,i)*dgamma(i)

end do

C----- Increment of deformation gradient associated with lattice stretching

do j=1,3

do i=1,3

if (i.eq.j) then

dvgrad(i,j)=delats(i)

else

dvgrad(i,j)=delats(i+j+1)

end if

end do

do j=1,3

do i=1,j

if (j.gt.i) then

ij2=i+j-2

if (mod(ij2,2).eq.1) then

term1=1.

else

term1=-1.

end if

do k=1,nslptl

dvgrad(i,j)=dvgrad(i,j)-term1*dgamma(k)*

2 slpspn(ij2,k)

dvgrad(j,i)=dvgrad(j,i)+term1*dgamma(k)*

2 slpspn(ij2,k)

end do

end if

end do

C----- Increment of resolved shear stress in a slip system: DTAUSP

do i=1,nslptl

dtausp(i)=0.

do j=1,6

dtausp(i)=dtausp(i)+ddemsd(j,i)*delats(j)

end do

C----- Update the resolved shear stress in a slip system:

do i=1,nslptl

stateold(km,2*nslptl+i)=stateold(km,2*nslptl+i)

2 +dtausp(i)

end do

C----- Increment of stress: DSTRES

do i=1,ndir+nshr

dstres(i)=-stressold(km,i)*dev

end do

do i=1,ndir

do j=1,ndir

dstres(i)=dstres(i)+d(i,j)*straininc(km,j)

end do

if (nshr.gt.0)then

do j=1,nshr

dstres(i)=dstres(i)+d(i,j+3)*

2 straininc(km,j+ndir)

end do

end if

do j=1,nslptl

dstres(i)=dstres(i)-ddemsd(i,j)*dgamma(j)

end do

if (nshr.gt.0) then

do i=1,nshr

do j=1,ndir

dstres(i+ndir)=dstres(i+ndir)+d(i+3,j)

2 *straininc(km,j)

end do

do j=1,nshr

dstres(i+ndir)=dstres(i+ndir)+d(i+3,j+3)*

2 straininc(km,j+ndir)

end do

do j=1,nslptl

dstres(i+ndir)=dstres(i+ndir)-ddemsd(i+3,j)

2 *dgamma(j)

end do

end if

C----- Update the stress: STRESSOLD

do i=1,ndir+nshr

stressold(km,i)=stressold(km,i)+dstres(i)

end do

C----- Increment of normal to a slip plane and a slip direction

do j=1,nslptl

do i=1,3

dspnor(i,j)=0.

dspdir(i,j)=0.

do k=1,3

dspnor(i,j)=dspnor(i,j)-slpnor(k,j)*dvgrad(k,i)

dspdir(i,j)=dspdir(i,j)+slpdir(k,j)*dvgrad(i,k)

end do

C----- Update the normal to a slip plane and a slip direction

idnor=3*nslptl

iddir=6*nslptl

do j=1,nslptl

do i=1,3

idnor=idnor+1

stateold(km,idnor)=stateold(km,idnor)+dspnor(i,j)

iddir=iddir+1

stateold(km,iddir)=stateold(km,iddir)+dspdir(i,j)

end do

C----- Total cumulative shear strains on all slip systems

do i=1,nslptl

stateold(km,9*nslptl+1)=stateold(km,9*nslptl+1)

2 +abs(dgamma(i))

end do

c----- update stressold to stressnew

do i=1,ndir+nshr

stressnew(km,i)=stressold(km,i)

end do

c----- update stateold to statenew for 1 - nstatev

do i=1,nstatev

statenew(km,i)=stateold(km,i)

end do

if (nshr.gt.1) then

save=straininc(km,5)

straininc(km,5)=straininc(km,6)

straininc(km,6)=save

save=stressnew(km,5)

stressnew(km,5)=stressnew(km,6)

stressnew(km,6)=save

end if

enddo

return

end

c----------------------------------------------------------------------

subroutine ludcmp (a, n, np, indx, d)

include 'vaba_param.inc'

parameter (nmax=200, tiny=1.0e-20)

dimension a(np,np), indx(n), vv(nmax)

d = 1.d0

do i = 1,n

aamax = 0.

do j = 1,n

if (abs(a(i,j)).gt.aamax) aamax = abs(a(i,j))

end do

if (aamax.eq.0.) pause 'singular matrix.'

vv(i) = 1./aamax

end do

do j = 1,n

do i = 1,j-1

sum = a(i,j)

do k = 1,i-1

sum = sum-a(i,k)*a(k,j)

end do

a(i,j) = sum

end do

aamax = 0.

do i = j,n

sum = a(i,j)

do k = 1,j-1

sum = sum-a(i,k)*a(k,j)

end do

a(i,j) = sum

dum = vv(i)*abs(sum)

if (dum.ge.aamax) then

imax = i

aamax = dum

end if

end do

if (j.ne.imax) then

do k = 1,n

dum = a(imax,k)

a(imax,k) = a(j,k)

a(j,k) = dum

end do

d = -d

vv(imax) = vv(j)

end if

indx(j) = imax

if (a(j,j).eq.0.) a(j,j) = tiny

if (j.ne.n) then

dum = 1./a(j,j)

do i = j+1,n

a(i,j) = a(i,j)*dum

end do

end if

end do

return

end

C----------------------------------------------------------------------

subroutine lubksb (a, n, np, indx, b)

include 'vaba_param.inc'

dimension a(np,np), indx(n), b(n)

ii = 0

do i = 1,n

ll = indx(i)

sum = b(ll)

b(ll) = b(i)

if (ii.ne.0) then

do j = ii,i-1

sum = sum-a(i,j)*b(j)

end do

else if (sum.ne.0.) then

ii = i

end if

b(i) = sum

end do

do i = n,1,-1

sum = b(i)

if (i.lt.n) then

do j = i+1,n

sum = sum-a(i,j)*b(j)

end do

end if

b(i) = sum/a(i,i)

end do

return

end

C----------------------------------------------------------------------

适合HCP mg合金的黄永刚Vumat的晶体塑性（之前推文对应的vumat版本，同时附上对应的材料参数）的图1

相关文件上传到知识星球，直接再知识星球下载对应的程序和材料文件，参考文献查看之前的推文

登录后免费查看全文

立即登录

技术邻APP
工程师必备

项目客服
培训客服
平台客服

TOP

适合HCP mg合金的黄永刚Vumat的晶体塑性（之前推文对应的vumat版本，同时附上对应的材料参数）

推荐阅读