黄永刚单晶塑性有限元umat子程序

SUBROUTINE UMAT(stress,statev,ddsdde,sse,spd,scd,

1 rpl, ddsddt, drplde, drpldt,

2 stran,dstran,time,dtime,temp,dtemp,predef,dpred,cmname,

3 ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt,

4 celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)

c WRITE (6,*) '

c NOTE: MODIFICATIONS TO *UMAT FOR ABAQUS VERSION 5.3 (14 APR '94) c

c (1) The list of variables above defining the *UMAT subroutine,

c an

d th

e first (standard) block o

f variables dimensioned below,

c have variable names adde

d compared to earlier ABAQUS versions.

c

c (2) The statement: include 'aba_param.inc' must be adde

d as below.

c

c (3) As of version 5.3, ABAQUS files use double precision only.

c The file aba_param.inc has a line "implicit real*8" and, since

c it is include

d in th

e main subroutine, it will define the variables

c there as double precision. But other subroutines still nee

d the

c definition "implicit real*8" since there may be variables that are

c not passe

d to them through th

e list or common block.

c

c (4) This is current as of version 5.6 of ABAQUS.

c

c (5) Note adde

d by J. W. Kysar (4 November 1997). This UMAT has been

c modifie

d to keep track of th

e cumulative shear strain in each

c individual slip system. This information is neede

d to correct an

c error in the implementation of the Bassani an

d Wu hardening law.

c Any line of code which has been adde

d or modified is preceded

c immediately by a line beginning CFIXA an

d succeeded by a line

c beginning CFIXB. Any comment line adde

d or modified will begin

c with CFIX.

c

c The hardening law by Bassani an

d Wu was implemented incorrectly.

c This law is a function of both hyperbolic secant square

d and hyperbolic

c tangent. However, the arguments of sech an

d tanh ar

e related to the *total*

c slip on individual slip systems. Formerly, the UMAT implemente

d this

c hardening law by using the *current* slip on each slip system. Therein

c lay the problem. The UMAT di

d not restrict th

e current slip to be a

c positive value. So when a slip with a negative sign was encountered, the

c term containing tanh le

d to a negativ

e hardening rate (since tanh is an

c od

d function).

c The UMA T has been fixe

d by adding stat

e variables to keep track o

f the

c *total* slip on each slip system by integrating up the absolute value

c of slip rates for each individual slip system. These "solution dependent

c variables" are available for postprocessing. The only require

d change

c in the input file is that the DEPV AR comman

d must b

e changed.

c

C----- Use single precision on Cray by

C (1) deleting the statement "IMPLICIT*8 (A-H,O-Z)";

C (2) changing "REAL*8 FUNCTION" to "FUNCTION";

C (3) changing double precision functions DSIGN to SIGN.

C

C----- Subroutines:

C

C ROTATION -- forming rotation matrix, i.e. the direction

C cosines of cubic crystal [100], [010] and [001]

C directions in global system at the initial

C state

C

C SLIPSYS -- calculating number of slip systems, unit

C vectors in slip directions and unit normals to

C slip planes in a cubic crystal at the initial

C state

C

C GSLPINIT -- calculating initial value of current strengths

C at initial state

C

C STRAINRA TE -- based on current values of resolved shear

C stresses and current strength, calculating

C shear strain-rates in slip systems

C

C LATENTHARDEN -- forming self- and latent-hardening matrix

C

C ITERATION -- generating arrays for the Newton-Rhapson

C iteration

C

C LUDCMP -- LU decomposition

C

C LUBKSB -- linear equation solver based on LU

C decomposition method (must call LUDCMP first) C----- Function subprogram:

C F -- shear strain-rates in slip systems

C----- Variables:

C

C STRESS -- stresses (INPUT & OUTPUT)

C Cauchy stresses for finite deformation

C STATEV -- solution dependent state variables (INPUT & OUTPUT) C DDSDDE -- Jacobian matrix (OUTPUT)

C----- Variables passed in for information:

C

C STRAN -- strains

C logarithmic strain for finite deformation

C (actually, integral of the symmetric part of velocity

C gradient with respect to time)

C DSTRAN -- increments of strains

C CMNAME -- name given in the *MATERIAL option

C NDI -- number of direct stress components

C NSHR -- number of engineering shear stress components

C NTENS -- NDI+NSHR

C NSTATV -- number of solution dependent state variables (as

C defined in the *DEPVAR option)

C PROPS -- material constants entered in the *USER MA TERIAL C option

C NPROPS -- number of material constants

C

C----- This subroutine provides the plastic constitutive relation of

C single crystals for finite element code ABAQUS. The plastic slip

C of single crystal obeys the Schmid law. The program gives the

C choice of small deformation theory and theory of finite rotation

C and finite strain.

C The strain increment is composed of elastic part and plastic

C part. The elastic strain increment corresponds to lattice

C stretching, the plastic part is the sum over all slip systems of

C plastic slip. The shear strain increment for each slip system is

C assumed a function of the ratio of corresponding resolved shear

C stress over current strength, and of the time step. The resolved

C shear stress is the double product of stress tensor with the slip

C deformation tensor (Schmid factor), and the increment of current

C strength is related to shear strain increments over all slip

C systems through self- and latent-hardening functions.

C----- The implicit integration method proposed by Peirce, Shih and

C Needleman (1984) is used here. The subroutine provides an option C of iteration to solve stresses and solution dependent state

C variables within each increment.

C----- The present program is for a single CUBIC crystal. However,

C this code can be generalized for other crystals (e.g. HCP,

C Tetragonal, Orthotropic, etc.). Only subroutines ROTATION and

C SLIPSYS need to be modified to include the effect of crystal

C aspect ratio.

C

C----- Important notice:

C

C (1) The number of state variables NSTATV must be larger than (or CFIX equal to) TEN (10) times the total number of slip systems in

C all sets, NSLPTL, plus FIVE (5)

CFIX NSTATV >= 10 * NSLPTL + 5

C Denote s as a slip direction and m as normal to a slip plane.

C Here (s,-m), (-s,m) and (-s,-m) are NOT considered

C independent of (s,m). The number of slip systems in each set

C could be either 6, 12, 24 or 48 for a cubic crystal, e.g. 12

C for {110}<111>.

C

C Users who need more parameters to characterize the

C constitutive law of single crystal, e.g. the framework

C proposed by Zarka, should make NSTATV larger than (or equal C to) the number of those parameters NPARMT plus nine times

C the total number of slip systems, NSLPTL, plus five

CFIX NSTATV >= NPARMT + 10 * NSLPTL + 5

C

C (2) The tangent stiffness matrix in general is not symmetric if

C latent hardening is considered. Users must declare "UNSYMM"

C in the input file, at the *USER MATERIAL card.

C

PARAMETER (ND=150)

C----- The parameter ND determines the dimensions of the arrays in

C this subroutine. The current choice 150 is a upper bound for a

C cubic crystal with up to three sets of slip systems activated.

C Users may reduce the parameter N

D to any number as long as larger C than or equal to the total number of slip systems in all sets.

C For example, if {110}<111> is the only set of slip system

C potentially activated, N

D could be taken as twelve (12).

c

include 'aba_param.inc'

c

CHARACTER*8 CMNAME

EXTERNAL F

dimension stress(ntens),statev(nstatv),

1 ddsdde(ntens,ntens),ddsddt(ntens),drplde(ntens),

2 stran(ntens),dstran(ntens),time(2),predef(1),dpred(1),

3 props(nprops),coords(3),drot(3,3),dfgrd0(3,3),dfgrd1(3,3)

DIMENSION ISPDIR(3), ISPNOR(3), NSLIP(3),

2 SLPDIR(3,ND), SLPNOR(3,ND), SLPDEF(6,ND),

3 SLPSPN(3,ND), DSPDIR(3,ND), DSPNOR(3,ND),

4 DLOCAL(6,6), D(6,6), ROTD(6,6), ROTATE(3,3),

5 FSLIP(ND), DFDXSP(ND), DDEMSD(6,ND),

6 H(ND,ND), DDGDDE(ND,6),

7 DSTRES(6), DELATS(6), DSPIN(3), DVGRAD(3,3),

8 DGAMMA(ND), DTAUSP(ND), DGSLIP(ND),

9 WORKST(ND,ND), INDX(ND), TERM(3,3), TRM0(3,3), ITRM(3)

DIMENSION FSLIP1(ND), STRES1(6), GAMMA1(ND), TAUSP1(ND),

2 GSLP1(ND), SPNOR1(3,ND), SPDIR1(3,ND), DDSDE1(6,6),

3 DSOLD(6), DGAMOD(ND), DTAUOD(ND), DGSPOD(ND),

4 DSPNRO(3,ND), DSPDRO(3,ND),

5 DHDGDG(ND,ND)

C----- NSLIP -- number of slip systems in each set

C----- SLPDIR -- slip directions (unit vectors in the initial state)

C----- SLPNOR -- normals to slip planes (unit normals in the initial

C state)

C----- SLPDEF -- slip deformation tensors (Schmid factors)

C SLPDEF(1,i) -- SLPDIR(1,i)*SLPNOR(1,i)

C SLPDEF(2,i) -- SLPDIR(2,i)*SLPNOR(2,i)

C SLPDEF(3,i) -- SLPDIR(3,i)*SLPNOR(3,i)

C SLPDEF(4,i) -- SLPDIR(1,i)*SLPNOR(2,i)+

C SLPDIR(2,i)*SLPNOR(1,i)

C SLPDEF(5,i) -- SLPDIR(1,i)*SLPNOR(3,i)+

C SLPDIR(3,i)*SLPNOR(1,i)

C SLPDEF(6,i) -- SLPDIR(2,i)*SLPNOR(3,i)+

C SLPDIR(3,i)*SLPNOR(2,i)

C where index i corresponds to the ith slip system

C----- SLPSPN -- slip spin tensors (only needed for finite rotation)

C SLPSPN(1,i) -- [SLPDIR(1,i)*SLPNOR(2,i)-

C SLPDIR(2,i)*SLPNOR(1,i)]/2

C SLPSPN(2,i) -- [SLPDIR(3,i)*SLPNOR(1,i)-

C SLPDIR(1,i)*SLPNOR(3,i)]/2

C SLPSPN(3,i) -- [SLPDIR(2,i)*SLPNOR(3,i)-

C SLPDIR(3,i)*SLPNOR(2,i)]/2

C where index i corresponds to the ith slip system

C----- DSPDIR -- increments of slip directions

C----- DSPNOR -- increments of normals to slip planes

C

C----- DLOCAL -- elastic matrix in local cubic crystal system

C----- D -- elastic matrix in global system

C----- ROTD -- rotation matrix transforming DLOCAL to D

C

C----- ROTATE -- rotation matrix, direction cosines of [100], [010]

C and [001] of cubic crystal in global system

C

C----- FSLIP -- shear strain-rates in slip systems

C----- DFDXSP -- derivatives of FSLIP w.r.t x=TAUSLP/GSLIP, where

C TAUSLP is the resolved shear stress and GSLIP is the C current strength

C

C----- DDEMSD -- double dot product of the elastic moduli tensor with

C the slip deformation tensor plus, only for finite

C rotation, the dot product of slip spin tensor with

C the stress

C

C----- H -- self- and latent-hardening matrix

C H(i,i) -- self hardening modulus of the ith slip

C system (no sum over i)

C H(i,j) -- latent hardening molulus of the ith slip

C system due to a slip in the jth slip system C (i not equal j)

C

C----- DDGDDE -- derivatice of the shear strain increments in slip

C systems w.r.t. the increment of strains

C

C----- DSTRES -- Jaumann increments of stresses, i.e. corotational

C stress-increments formed on axes spinning with the

C material

C----- DELATS -- strain-increments associated with lattice stretching

C DELATS(1) - DELATS(3) -- normal strain increments C DELATS(4) - DELATS(6) -- engineering shear strain C increments

C----- DSPIN -- spin-increments associated with the material element

C DSPIN(1) -- component 12 of the spin tensor

C DSPIN(2) -- component 31 of the spin tensor

C DSPIN(3) -- component 23 of the spin tensor

C

C----- DVGRAD -- increments of deformation gradient in the current

C state, i.e. velocity gradient times the increment of

C time

C

C----- DGAMMA -- increment of shear strains in slip systems

C----- DTAUSP -- increment of resolved shear stresses in slip systems

C----- DGSLIP -- increment of current strengths in slip systems

C

C

C----- Arrays for iteration:

C

C FSLIP1, STRES1, GAMMA1, TAUSP1, GSLP1 , SPNOR1, SPDIR1,

C DDSDE1, DSOL

D , DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO, C DHDGDG

C

C

C----- Solution dependent state variable STATEV:

C Denote the number of total slip systems by NSLPTL, which

C will be calculated in this code.

C

C Array STATEV:

C 1 - NSLPTL : current strength in slip systems

C NSLPTL+1 - 2*NSLPTL : shear strain in slip systems

C 2*NSLPTL+1 - 3*NSLPTL : resolved shear stress in slip systems

C

C 3*NSLPTL+1 - 6*NSLPTL : current components of normals to slip

C slip planes

C 6*NSLPTL+1 - 9*NSLPTL : current components of slip directions

C

CFIX 9*NSLPTL+1 - 10*NSLPTL : total cumulative shear strain on each

CFIX slip system (sum of the absolute

CFIX values of shear strains in each slip

CFIX system individually)

CFIX

CFIX 10*NSLPTL+1 : total cumulative shear strain on all

C slip systems (sum of the absolute

C values of shear strains in all slip

C systems)

C

CFIX 10*NSLPTL+2 - NSTA TV-4 : additional parameters users may need

C to characterize the constitutive law

C of a single crystal (if there are

C any).

C

C NSTATV-3 : number of slip systems in the 1st set

C NSTATV-2 : number of slip systems in the 2nd set

C NSTATV-1 : number of slip systems in the 3rd set

C NSTATV : total number of slip systems in all

C sets

C

C

C----- Material constants PROPS:

C

C PROPS(1) - PROPS(21) -- elastic constants for a general elastic

C anisotropic material

C

C isotropic : PROPS(i)=0 for i>2

C PROPS(1) -- Young's modulus

C PROPS(2) -- Poisson's ratio

C

C cubic : PROPS(i)=0 for i>3

C PROPS(1) -- c11

C PROPS(2) -- c12

C PROPS(3) -- c44

C

C orthotropic : PORPS(i)=0 for i>9

C PROPS(1) - PROPS(9) are D1111, D1122, D2222, C D1133, D2233, D3333, D1212, D1313, D2323,

C respectively, which has the same definition

C as ABAQUS for orthotropic materials

C (see *ELASTIC card)

C

C anisotropic : PROPS(1) - PROPS(21) are D1111, D1122,

C D2222, D1133, D2233, D3333, D1112, D2212,

C D3312, D1212, D1113, D2213, D3313, D1213,

C D1313, D1123, D2223, D3323, D1223, D1323,

C D2323, respectively, which has the same

C definition as ABAQUS for anisotropic

C materials (see *ELASTIC card)

C

C

C PROPS(25) - PROPS(56) -- parameters characterizing all slip

C systems to be activated in a cubic

C crystal

C

C PROPS(25) -- number of sets of slip systems (maximum 3),

C e.g. (110)[1-11] and (101)[11-1] are in the

C same set of slip systems, (110)[1-11] and

C (121)[1-11] belong to different sets of slip

C systems

C (It must be a real number, e.g. 3., not 3 !)

C

C PROPS(33) - PROPS(35) -- normal to a typical slip plane in

C the first set of slip systems,

C e.g. (1 1 0)

C (They must be real numbers, e.g.

C 1. 1. 0., not 1 1 0 !)

C PROPS(36) - PROPS(38) -- a typical slip direction in the

C first set of slip systems, e.g.

C [1 1 1]

C (They must be real numbers, e.g.

C 1. 1. 1., not 1 1 1 !)

C

C PROPS(41) - PROPS(43) -- normal to a typical slip plane in

C the second set of slip systems

C (real numbers)

C PROPS(44) - PROPS(46) -- a typical slip direction in the

C second set of slip systems

C (real numbers)

C

C PROPS(49) - PROPS(51) -- normal to a typical slip plane in

C the third set of slip systems

C (real numbers)

C PROPS(52) - PROPS(54) -- a typical slip direction in the

C third set of slip systems

C (real numbers)

C

C

C PROPS(57) - PROPS(72) -- parameters characterizing the initial

C orientation of a single crystal in

C global system

C The directions in global system and directions in local

C cubic crystal system of two nonparallel vectors are needed

C to determine the crystal orientation.

C

C PROPS(57) - PROPS(59) -- [p1 p2 p3], direction of first

C vector in local cubic crystal

C system, e.g. [1 1 0]

C (They must be real numbers, e.g.

C 1. 1. 0., not 1 1 0 !)

C PROPS(60) - PROPS(62) -- [P1 P2 P3], direction of first

C vector in global system, e.g.

C [2. 1. 0.]

C (It does not have to be a unit

C vector)

C

C PROPS(65) - PROPS(67) -- direction of second vector in

C local cubic crystal system (real C numbers)

C PROPS(68) - PROPS(70) -- direction of second vector in

C global system

C

C

C PROPS(73) - PROPS(96) -- parameters characterizing the visco-

C plastic constitutive law (shear

C strain-rate vs. resolved shear

C stress), e.g. a power-law relation

C

C PROPS(73) - PROPS(80) -- parameters for the first set of

C slip systems

C PROPS(81) - PROPS(88) -- parameters for the second set of C slip systems

C PROPS(89) - PROPS(96) -- parameters for the third set of

C slip systems

C

C

C PROPS(97) - PROPS(144)-- parameters characterizing the self-

C and latent-hardening laws of slip

C systems

C

C PROPS(97) - PROPS(104)-- self-hardening parameters for the C first set of slip systems

C PROPS(105)- PROPS(112)-- latent-hardening parameters for C the first set of slip systems and C interaction with other sets of

C slip systems

C

C PROPS(113)- PROPS(120)-- self-hardening parameters for the C second set of slip systems

C PROPS(121)- PROPS(128)-- latent-hardening parameters for C the second set of slip systems C and interaction with other sets C of slip systems

C

C PROPS(129)- PROPS(136)-- self-hardening parameters for the C third set of slip systems

C PROPS(137)- PROPS(144)-- latent-hardening parameters for C the third set of slip systems and C interaction with other sets of

C slip systems

C

C

C PROPS(145)- PROPS(152)-- parameters characterizing forward time C integration scheme and finite

C deformation

C

C PROPS(145) -- parameter theta controlling the implicit

C integration, which is between 0 and 1

C 0. : explicit integration

C 0.5 : recommended value

C 1. : fully implicit integration

C

C PROPS(146) -- parameter NLGEOM controlling whether the C effect of finite rotation and finite strain

C of crystal is considered,

C 0. : small deformation theory

C otherwise : theory of finite rotation and

C finite strain

C

C

C PROPS(153)- PROPS(160)-- parameters characterizing iteration

C method

C

C PROPS(153) -- parameter ITRATN controlling whether the

C iteration method is used,

C 0. : no iteration

C otherwise : iteration

C

C PROPS(154) -- maximum number of iteration ITRMAX

C

C PROPS(155) -- absolute error of shear strains in slip

C systems GAMERR

C

C----- Elastic matrix in local cubic crystal system: DLOCAL

DO J=1,6

DO I=1,6

DLOCAL(I,J)=0.

END DO

END DO

CHECK=0.

DO J=10,21

CHECK=CHECK+ABS(PROPS(J))

END DO

IF (CHECK.EQ.0.) THEN

DO J=4,9

CHECK=CHECK+ABS(PROPS(J))

END DO

IF (CHECK.EQ.0.) THEN

IF (PROPS(3).EQ.0.) THEN

C----- Isotropic material

GSHEAR=PROPS(1)/2./(1.+PROPS(2))

E11=2.*GSHEAR*(1.-PROPS(2))/(1.-2.*PROPS(2))

E12=2.*GSHEAR*PROPS(2)/(1.-2.*PROPS(2))

DO J=1,3

DLOCAL(J,J)=E11

DO I=1,3

IF (I.NE.J) DLOCAL(I,J)=E12

END DO

DLOCAL(J+3,J+3)=GSHEAR

END DO

ELSE

C----- Cubic material

DO J=1,3

DLOCAL(J,J)=PROPS(1)

DO I=1,3

IF (I.NE.J) DLOCAL(I,J)=PROPS(2)

END DO

DLOCAL(J+3,J+3)=PROPS(3)

END DO

END IF

ELSE

C----- Orthotropic metarial

DLOCAL(1,1)=PROPS(1)

DLOCAL(1,2)=PROPS(2)

DLOCAL(2,1)=PROPS(2)

DLOCAL(2,2)=PROPS(3)

DLOCAL(1,3)=PROPS(4)

DLOCAL(3,1)=PROPS(4)

DLOCAL(2,3)=PROPS(5)

DLOCAL(3,2)=PROPS(5)

DLOCAL(3,3)=PROPS(6)

DLOCAL(4,4)=PROPS(7)

DLOCAL(5,5)=PROPS(8)

DLOCAL(6,6)=PROPS(9)

END IF

ELSE

C----- General anisotropic material

ID=0

DO J=1,6

DO I=1,J

ID=ID+1

DLOCAL(I,J)=PROPS(ID)

DLOCAL(J,I)=DLOCAL(I,J)

END DO

END DO

END IF

C----- Rotation matrix: ROTA TE, i.e. direction cosines of [100], [010]

C and [001] of a cubic crystal in global system

C

CALL ROTATION (PROPS(57), ROTA TE)

C----- Rotation matrix: ROTD to transform local elastic matrix DLOCAL C to global elastic matrix D

C

DO J=1,3

J1=1+J/3

J2=2+J/2

DO I=1,3

I1=1+I/3

I2=2+I/2

ROTD(I,J)=ROTATE(I,J)**2

ROTD(I,J+3)=2.*ROTATE(I,J1)*ROTA TE(I,J2)

ROTD(I+3,J)=ROTA TE(I1,J)*ROTATE(I2,J)

ROTD(I+3,J+3)=ROTA TE(I1,J1)*ROTA TE(I2,J2)+

2 ROTA TE(I1,J2)*ROTA TE(I2,J1)

END DO

END DO

C----- Elastic matrix in global system: D

C {D} = {ROTD} * {DLOCAL} * {ROTD}transpose

C

DO J=1,6

DO I=1,6

D(I,J)=0.

END DO

END DO

DO J=1,6

DO I=1,J

DO K=1,6

DO L=1,6

D(I,J)=D(I,J)+DLOCAL(K,L)*ROTD(I,K)*ROTD(J,L) END DO

END DO

D(J,I)=D(I,J)

END DO

END DO

C----- Total number of sets of slip systems: NSET

NSET=NINT(PROPS(25))

IF (NSET.LT.1) THEN

WRITE (6,*) '***ERROR - zero sets of slip systems'

STOP

ELSE IF (NSET.GT.3) THEN

WRITE (6,*)

2 '***ERROR - more than three sets of slip systems'

STOP

END IF

C----- Implicit integration parameter: THETA

THETA=PROPS(145)

C----- Finite deformation ?

C----- NLGEOM = 0, small deformation theory

C otherwise, theory of finite rotation and finite strain, Users C must declare "NLGEOM" in the input file, at the *STEP card C

IF (PROPS(146).EQ.0.) THEN

NLGEOM=0

ELSE

NLGEOM=1

END IF

C----- Iteration?

C----- ITRATN = 0, no iteration

C otherwise, iteration (solving increments of stresses and

C solution dependent state variables)

C

IF (PROPS(153).EQ.0.) THEN

ITRATN=0

ELSE

ITRATN=1

END IF

ITRMAX=NINT(PROPS(154))

GAMERR=PROPS(155)

NITRTN=-1

DO I=1,NTENS

DSOLD(I)=0.

END DO

DO J=1,ND

DGAMOD(J)=0.

DTAUOD(J)=0.

DGSPOD(J)=0.

DO I=1,3

DSPNRO(I,J)=0.

DSPDRO(I,J)=0.

END DO

END DO

C----- Increment of spin associated with the material element: DSPIN C (only needed for finite rotation)

C

IF (NLGEOM.NE.0) THEN

DO J=1,3

DO I=1,3

TERM(I,J)=DROT(J,I)

TRM0(I,J)=DROT(J,I)

END DO

TERM(J,J)=TERM(J,J)+1.D0

TRM0(J,J)=TRM0(J,J)-1.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)

END IF

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

DEV=0.D0

DO I=1,NDI

DEV=DEV+DSTRAN(I)

END DO

C----- Iteration starts (only when iteration method is used)

1000 CONTINUE

C----- Parameter NITRTN: number of iterations

C NITRTN = 0 --- no-iteration solution

C

NITRTN=NITRTN+1

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

IF (STATEV(1).EQ.0.) THEN

C----- Initial state

C

C----- Generating the following parameters and variables at initial

C state:

C Total number of slip systems in all the sets NSLPTL

C Number of slip systems in each set NSLIP

C Unit vectors in initial slip directions SLPDIR

C Unit normals to initial slip planes SLPNOR

C

NSLPTL=0

DO I=1,NSET

ISPNOR(1)=NINT(PROPS(25+8*I))

ISPNOR(2)=NINT(PROPS(26+8*I))

ISPNOR(3)=NINT(PROPS(27+8*I))

ISPDIR(1)=NINT(PROPS(28+8*I))

ISPDIR(2)=NINT(PROPS(29+8*I))

ISPDIR(3)=NINT(PROPS(30+8*I))

CALL SLIPSYS (ISPDIR, ISPNOR, NSLIP(I), SLPDIR(1,NSLPTL+1),

2 SLPNOR(1,NSLPTL+1), ROTATE)

NSLPTL=NSLPTL+NSLIP(I)

END DO

IF (ND.LT.NSLPTL) THEN

WRITE (6,*)

2 '***ERROR - parameter ND chosen by the present user is less than

3 the total number of slip systems NSLPTL'

STOP

END IF

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)+SLPDIR(2,J)*SLPNOR(1,J)

SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)

SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J) END DO

C----- Initial value of state variables: unit normal to a slip plane

C and unit vector in a slip direction

C

STATEV(NSTATV)=FLOAT(NSLPTL)

DO I=1,NSET

STA TEV(NSTA TV-4+I)=FLOAT(NSLIP(I))

END DO

IDNOR=3*NSLPTL

IDDIR=6*NSLPTL

DO J=1,NSLPTL

DO I=1,3

IDNOR=IDNOR+1

STA TEV(IDNOR)=SLPNOR(I,J)

IDDIR=IDDIR+1

STA TEV(IDDIR)=SLPDIR(I,J)

END DO

END DO

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

C

CALL GSLPINIT (STATEV(1), NSLIP, NSLPTL, NSET, PROPS(97))

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

CFIX-- Initial value of cumulative shear strain in each slip systems

DO I=1,NSLPTL

STA TEV(NSLPTL+I)=0.

CFIXA

STA TEV(9*NSLPTL+I)=0.

CFIXB

END DO

CFIXA

STATEV(10*NSLPTL+1)=0.

CFIXB

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

DO I=1,NSLPTL

TERM1=0.

DO J=1,NTENS

IF (J.LE.NDI) THEN

TERM1=TERM1+SLPDEF(J,I)*STRESS(J)

ELSE

TERM1=TERM1+SLPDEF(J-NDI+3,I)*STRESS(J) END IF

END DO

STA TEV(2*NSLPTL+I)=TERM1

END DO

ELSE

C----- Current stress state

C

C----- Copying from the array of state variables STA TVE the following C parameters and variables at current stress state:

C Total number of slip systems in all the sets NSLPTL

C Number of slip systems in each set NSLIP

C Current slip directions SLPDIR

C Normals to current slip planes SLPNOR

C

NSLPTL=NINT(STA TEV(NSTATV))

DO I=1,NSET

NSLIP(I)=NINT(STATEV(NSTATV-4+I))

END DO

IDNOR=3*NSLPTL

IDDIR=6*NSLPTL

DO J=1,NSLPTL

DO I=1,3

IDNOR=IDNOR+1

SLPNOR(I,J)=STATEV(IDNOR)

IDDIR=IDDIR+1

SLPDIR(I,J)=STA TEV(IDDIR)

END DO

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)+SLPDIR(2,J)*SLPNOR(1,J)

SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)

SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J) END DO

END IF

C----- Slip spin tensor: SLPSPN (only needed for finite rotation)

IF (NLGEOM.NE.0) THEN

DO J=1,NSLPTL

SLPSPN(1,J)=0.5*(SLPDIR(1,J)*SLPNOR(2,J)-

2 SLPDIR(2,J)*SLPNOR(1,J))

SLPSPN(2,J)=0.5*(SLPDIR(3,J)*SLPNOR(1,J)-

2 SLPDIR(1,J)*SLPNOR(3,J))

SLPSPN(3,J)=0.5*(SLPDIR(2,J)*SLPNOR(3,J)-

2 SLPDIR(3,J)*SLPNOR(2,J))

END DO

END IF

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

C deformation tensor (Schmid factors) plus, only for finite

C rotation, the dot product of slip spin tensor with the stress:

C DDEMSD

C

DO J=1,NSLPTL

DO I=1,6

DDEMSD(I,J)=0.

DO K=1,6

DDEMSD(I,J)=DDEMSD(I,J)+D(K,I)*SLPDEF(K,J)

END DO

END DO

END DO

IF (NLGEOM.NE.0) THEN

DO J=1,NSLPTL

DDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(1,J)*STRESS(1)

DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(2,J)*STRESS(1)

IF (NDI.GT.1) THEN

DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(1,J)*STRESS(2)

纯钛塑性变形行为的晶体塑性有限元模拟

纯钛塑性变形行为的晶体塑性有限元模拟本文采用晶体塑性有限元(Crystal Plasty Finite Element Method,CP-FEM)为基础,结合Voronoi图表技术建立了包含纯钛材料基本参数的单晶体和多晶体塑性变形有限元模型,从细观角度研究了晶体的变形行为。首先利用晶体塑性有限元方法,建立了纯钛单晶塑性变形模型,对钛单晶体的拉伸变形行为进行了研究,模拟结果表明在边界约束条件下,单晶体拉伸变形的同时发生了取向的旋转,其滑移系的滑移方向会逐步趋向于与拉伸方向平行。 利用建立的单晶体模型研究了晶体初始取向与晶体滑移的宏观力学响关系,结果表明,软位向方向拉伸的拉伸应力小于硬位向方向拉升应力,模拟变形过程符合Schmid定律。通过Voronoi图表信息建立了TA2多晶体模型,并模拟了拉伸变形,对照真实材料拉伸实验,分析了晶体变形过程的外形、组织形貌及应力应变关系,模拟与实验结果接近。 分段表征了TA2多晶体模型变形的应力应变演化过程,描述了多晶体塑性变形弹性阶段、屈服阶段、硬化阶段与颈缩阶段应力应变分布,详细分析了晶粒变形的不均匀性的微观机制,说明了桔皮组织的形成过程。分析了不同晶粒尺寸对多晶体塑性变形行为的影响,研究结果表明,随着晶粒度减小,晶体屈服强度增大,硬化曲线变陡,结果符合Hall-Petch关系,预测TA2纯钛的Hall-Petch常数为σ0=64.9MPa,K=473.3MPaμm1/2。 利用设定同取向晶粒建立织构模型的方法研究了初始织构对多晶体变形各向异性及应力应变关系的影响,结果表明随着织构强度的升高,多晶体变形不均匀性增大,多晶体屈服强度升高,硬化曲线变陡。

有限元与数值方法-讲稿19 弹塑性增量有限元分析课件

材料非线性问题有限元方法 教学要求和内容 1.掌握弹塑性本构关系和塑性力学的基本法则； 2.掌握弹塑性增量分析的有限元格式； 3.学习常用非线性方程组的求解方法： （1）直接迭代法； （2） Newton-Raphson 方法，修正的N－R 方法； （3）增量法等。 请大家预习，争取对相关内容有大概的了解和把握。

弹塑性增量有限元分析 一．材料弹塑性行为的描述 弹塑性材料进入塑性的特点：存在 不可恢复的塑性变形； 卸载时：非线性弹性材料按原路径 卸载； 弹塑性材料按不同的路径卸载，并 且有残余应变，称为塑性应变。

1．单向加载 1） 弹性阶段: 卸载时不留下残余变形; 2） 初始屈服：s σσ= 3） 强化阶段：超过初始屈服之后,按弹性规律卸载，再加载弹性范 围扩大：ss σσ'>，s σ'为相继屈服应力。

4） 鲍氏现象（Bauschinger ）： 二．塑性力学的基本法则 1．初始屈服准则： 00(,)0ij F k σ= 已经建立了多种屈服准则： （1） V . Mises 准则：000(,)()0ij ij F k f k σσ=-= 2 2 001 1 ()(),()2 3ij ij ij s f s s J k σσ===第二应力不变量1122221 ,() 3 ij ij ij m m s σδσσσσσ=-=++偏应力张量：平均应力： （2） Tresca 准则（最大剪应力准则）： 0max ()0ij s F S ττ=-=

2．流动法则 V . Mises 流动法则： 0(,)()ij ij p ij ij ij F k f d d d σσελ λ σσ??==??， 0d λ> 待定有限量 塑性应变增量 p ij d ε 沿屈服面当前应力点的法线方向增加。 因此，称为法向流动法则。 3．硬化法则： （1）各向同性硬化：(,)()0ij ij F k f k σσ=-=

弹塑性及有限元题目整理

一、应力 1. 什么是偏应力状态？什么是静水压力状态？举例说明？ 静水压力状态时指微六面体的每个面只有正应力作用，偏应力状态是从应力状态中扣除静水压力后剩下的部分。 2. 应力边界条件所描述的物理本质是什么？ 物体边界点的平衡条件。 3. 为什么定义物体内部应力状态的时候要采取在一点的领域取极限的方法？ 不规则，内部受力不一样。 4.Pie平面上的点所代表的应力状态有何特点？ 该平面上任意一点的所代表值的应力状态1+2+3=0，为偏应力状态，且该平面上任一法线所代表的应力状态其应力解不唯一。 5.固体力学解答必须满足的三个条件是什么？可否用其他条件代替？ 可以。能量原理处于整个系统。 6. 解释应力空间中为什么应力状态不能位于加载面之外？ 保证位移单值连续。连续体的形变分量、、不是互相独立的，而是相关，否则导致位移不单值，不连续。 二、应变 1．从数学和物理的不同角度，阐述相容方程的意义。 从数学角度看，由于几何方程是6个，而待求的位移分量是3个，方程数目多于未知函数的数目，求解出的位移不单值。从物理角度看，物体各点可以想象成微小六面体，微单元体之间就会出现"裂缝"或者相互"嵌入"，即产生不连续。 2.两个材料不同、但几何形状、边界条件及体积力（且体积力为常数）等都完全相同的线弹性平面问题，它们的应力分布是否相同？为什么？ 相同。应力分布受到平衡方程、变形协调方程及力边界条件，未涉及本构方程，与材料性质无关。 3.应力状态是否可以位于加载面外？为什么？ 不可以。保证位移单值连续。连续体的形变分量、、不是互相独立的，而是相关，否则导致位移不单值，不连续。 4．给定单值连续的位移函数，通过几何方程可求出应变分量，问这些应变分量是否满足变形协调方程？为什么？ 满足。根据几何方程求出各应变分量，则变形协调方程自然满足，因为变形协调方程本身是从几何方程中推导出来的。

刚塑性有限元数值模拟中产生误差的原因及改进方法(精)

刚塑性有限元数值模拟中产生误差的原因及改进方法 1 引言 塑性加工过程的有限元数值模拟，可以获得金属变形的详细规律，如网格变形、速度场、应力和应变场的分布规律，以及载荷-行程曲线。通过对模拟结果的可视化分析，可以在现有的模具设计上预测金属的流动规律，包括缺陷的产生（如角部充不满、折叠、回流和断裂等）。利用得到的力边界条件对模具进行结构分析，从而改进模具设计，提高模具设计的合理性和模具的使用寿命，减少模具重新试制的次数。在制造技术飞速发展、市场竞争日益加剧的今天，塑性加工过程的计算机模拟可在模具虚拟设计、制造阶段就能充分检验模具设计的合理性，减少新产品模具的开发研制时间，对用户需求做出快速响应，提高市场竞争能力。由此可见，金属成型过程的有限元模拟已是模具计算机集成制造系统中必不可少的模具设计检验环节。 金属成形工艺分体积成形和板料成形两大类，相应地，用于分析其流动规律的有限元法也分为两类，即：刚塑性、刚粘塑性有限元和弹塑性有限元。体积成形中的挤压成形和锻造成形在实际生产中应用很广，中外学者在这方面进行了很多研究，其中二维模拟技术已相当成熟，三维模拟是目前的世界研究热点。刚塑性、刚粘塑性有限元模拟能否对模具设计的合理性做出可靠校验，取决于模拟的精度和效率。作者结合从事二维塑性有限元模拟的经验和当前的三维塑性有限元模拟系统开发的实践，对刚塑性、刚粘塑性有限元模拟过程中产生误差的原因进行了全面的详细分析，并提出相应的解决方法，同时以具体实例说明。 2 刚塑性、刚粘塑性有限元模拟中产生误差的原因及改进方法 2.1 刚塑性有限元法求解的数学基础 刚塑性有限元法是假设材料具有刚塑性的特点，把实际的加工过程定义为边值问题，从刚塑性材料的变分原理或上界定理出发，接有限元模式把能耗率表示为节点速度的非线性函数，利用数学上的最优化原理，在给定变形体某些表面的力边界条件和速度边界条件的情况下，求满足平衡方程、本构方程和体积不变条件的速度场和应力场。速度场的真实解使以动可容速度场建立的能量泛函取极小值。但所得到的塑性力学的微分方程组一般不能用解析法求解，常采用数值解近似，而采用数值解，则会出现各种误差。误差取决于所用的数值方法。下述处理方式易引起系统误差。 2.1.1时间和空间的离散化

纯钛塑性变形行为的晶体塑性有限元模拟

硕士学位论文 纯钛塑性变形行为的晶体塑性 有限元模拟 THE FINITE ELEMENT SIMULATION ON PLASTIC BEHAVIOR OF PURE TITANIUM 黄晓华 2010年7月

国内图书分类号：TG113.12 学校代码：10213国际图书分类号: 669.017.3 密级：公开 工学硕士学位论文 纯钛塑性变形行为的晶体塑性 有限元模拟 硕 士 研 究 生： 黄晓华 导 师： 来忠红 副教授 申 请 学 位： 工学硕士 学 科、专 业： 材 料 学 所 在 单 位： 材料科学与工程学院 答 辩 日 期： 2010年7月 授予学位单位： 哈尔滨工业大学

Classified Index: TG113.12 U.D.C: 669.017.3 Dissertation for the Master Degree in Engineering THE FINITE ELEMENT SIMULATION ON PLASTIC BEHAVIOR OF PURE TITANIUM Candidate：Huang Xiaohua Supervisor：Associate Prof. Lai Zhonghong Academic Degree Applied for：Master of Engineering Speciality：Materials Science Affiliation：School of Materials Science and Engineering Date of Defence:July, 2010 Degree-Conferring-Institution:Harbin Institute of Technology

弹塑性力学有限单元法

中南大学2014年博士研究生入学考试 《弹塑性力学有限单元法》考试大纲 本考试大纲由交通运输工程学院教授委员会于2013年7月通过。 I.考试性质 弹塑性力学有限单元法是我校“载运工具运用工程”专业博士生入学考试的专业基础课，它是为我校招收本专业博士生而实施的具有选拔功能的水平考试；其目的是科学、公平、有效地测试考生掌握弹性力学、塑性力学及有限单元数值方法课程的基本知识、基本理论，以及相关理论和方法分析解决实际问题的能力；评价的标准是高等学校优秀硕士毕业生能达到的及格或及格以上水平，以保证被录取者能较好的掌握了本专业必备的基础知识。 II.考查目标 弹塑性力学有限单元法课程考试弹性力学、塑性力学及有限单元数值方法等内容，重点在检查力学基本概念与基本方法的掌握和应用，难度适中，覆盖主要章节，能区分学生优劣层次。要求考生：（1）掌握弹塑性力学的基本知识、结构有限元分析的基本方法和过程，要求学生具备使用有限元方法进行车辆结构强度分析的能力。 Ⅲ.考试形式和试卷结构 1、试卷满分及考试时间 本试卷满分为100 分，考试时间为180 分钟 2、答题方式 答题方式为闭卷，笔试。 3、试卷内容结构 弹性力学约30 % 30 有限单元法约50 % 50

塑性力学基本理论约20 % 20 Ⅳ.考查内容 1. 弹性力学 （1）掌握弹性力学问题基本方程及边界条件。 （2）掌握应力理论及变形理论、二阶张量的坐标转换； （3）掌握使用位移法和应力法求解弹性力学问题； （4）掌握使用半逆解法求解简单平面问题； 2. 有限单元法 （1）掌握有限元方法的基本概念； （2）掌握平面、空间及等参单元分析的过程 （3）掌握有限单元位移模式的选取、刚度矩阵数值积分方法；（4）掌握结构刚度矩阵性质、边界条件处理； （5）掌握薄板弯曲问题有限元分析方法； （6）掌握车辆典型结构有限元分析的步骤和处理技巧； 3. 塑性力学 （1）掌握塑性力学的基本概念； （2）掌握Tresca和Mises屈服条件； （3）掌握几种常用的弹塑性力学模型； （4）掌握应力空间和屈服曲面的概念、加载曲面和塑性流动法则；

弹塑性有限元方法

第三章 弹塑性有限元方法的实施 §3.1 增量平衡方程和切线刚度矩阵 1、 分段线性化的求解思想 塑性变形的特点决定了塑性本构关系的非线性和多值性，上面由塑性增量理论给 出了塑性应力—应变关系{}{}ep d D d σε=???? 其中 [][] {}{}[]{}[]{} T ep T F F D D D D F F A D σσ σ σ ????=- ??+ ?????? 说明当前应力状态不仅与当前应变有关，而且和达到这一变形状态的路径（加载历史）有关。这里包含了屈服准则、强化条件和加卸载准则。 由此，对物理非线性问题，通常采用分段线性化的纯增量法和逐次迭代的方法求解。即将加载过程分成若干个增量步，选择其中任意一个增量步建立它的增量平衡方程并求解，对整个过程的求解有普遍意义。 2、 增量平衡方程和切线刚度矩阵 设t 时刻（加载至i -1步终），结构（单元）在当前载荷（广义体力{}v f 和表面力{}s f ） 的作用下处于平衡状态，此时物体内一点的应力、应变状态为{}{}σε、。在此基础上，施加一个载荷增量{}{}v s f f ??和，即从t t t →+?时刻，则在体内必然引起一个位移增量{}u ?和相应的{}σ?、{}ε?，只要{}{}v s f f ??和足够小，就有{}{}ep D σε?=?????。 倘若初始状态{}σ已知，加载过程已知，则ep D ????可以确定（即p ij d ε?可以确定，然后 可在硬化曲线上得到1p ε所对应的硬化系数）于是上面的方程成为线性的。在t t t →+?这一增量过程中，应用于虚功原理可得到如下虚功方程： ()()()0e e T T T V V s s V S f f u dV f f u dS σσδεδδ??+?-+??-+??=?? ?? （1） 根据小变形几何关系u N q B q ε?=??=?和，再由虚位移()q δ?的任意性，并设 ()()e e T T v v s s V S P P N f f dV N f f dS +?= +?+ +?? ? ，展开后，其中单元在t 时刻载荷等效节点 力：e e T T v s V S P N f dV N f dS = + ? ? ；t ?内增量载荷的等效力e e T T v s V S P N f dV N f dS ?= ?+ ?? ? 。

塑性成形过程中的有限元法

塑性成形过程中的有限元法 金属塑性成形技术是现代化制造业中金属加工的重要方法之一。它是金属材料在模具和锻压设备作用下发生变形，获得所需要求的形状、尺寸和性能的制件的加工过程。金属成形件在汽车、飞机仪表、机械设备等产品的零部件中占有相当大的比例。由于其具有生产效率高，生产费用低的特点，适合于大批量生产，是现代高速发展的制造业的重要成形工艺。据统计，在发达国家中，金属塑性成形件的产值在国民经济中的比重居行业之首，在我国也占有相当大的比例。 随着现代制造业的高速发展，对塑性成形工艺分析和模具设计方面提出了更高的要求。若工艺分析不完善、模具设计不合理或材料选择不当，则会造成产品达不到质量要求，造成大量的次品和废品，增加了模具的设计制造时间和费用。为了防止缺陷的产生，以提高产品质量，降低产品成本，国内外许多大公司企业及大专院校和研究机构对塑性成形件的性能、成形过程中的应力应变分布及变化规律进行了大量的理论分析、实验研究与数值计算，力图发现各种制件、产品成形工艺所遵循的共同规律以及力学失效所反映的共同特征。由于塑性成形工艺影响因素甚多，有些因素如摩擦与润滑、变形过程中材料的本构关系等机理尚未被人们完全认识和掌握，因而到目前为止还未能对各种材料各种形状的制件成形过程作出准确的定量判定。正因为大变形机理非常复杂，使得塑性成形研究领域一直成为一个充满挑战和机遇的领域。 一般来说，产品研究与开发的目标之一就是确定生产高质量产品的优化准则，而不同的产品要求不同的优化准则，建立适当的优化准则需要对产品制造过程的全面了解。如果不掌握诸如摩擦条件、材料性能、工件几何形状、成形力等工艺参数对成形过程的影响，就不可能正确地设计模具和选择加工设备，更无法预测和防止缺陷的生成。在传统工艺分析和模具设计中，主要还是依靠工程类比和设计经验，经过反复试模修模，调整工艺参数以期望消除成形过程中的产品缺陷如失稳起皱、充填不满、局部破裂等。仅仅依靠类比和传统的经验工艺分析和模具设计方法已无法满足高速发展的现代金属加工工业的要求。因此，现代金属成形工艺分析过程中，建立适当的“过程模拟”非常重要。随着计算机技术的发展，人们已经认识到数值模拟在金属成形工程中的重要价值，这一领域已成为现代国内外学者的研究热点。 应用塑性成形的数值模拟方法主要有上限法（Upper Bound Method）、边界元法(Boundary Element Method)和有限元法（Finite Element Method）。上限元法常用于分析较为简单的准稳态变形问题；而边界元法主要用于模具设计分析和温度计算。对于大变形的体积成形和板料成形，变形过程常呈非稳态，形状、边界、材料性质等都会发生很大的变化，有限元法可由实验和理论方法给出的本构关系、边界条件、摩擦关系式，按变分原理推导出场方程，根据离散技术建立计算模型，从而实现对复杂成形问题进行数值模拟。分析成形过程中的应力应变分布及其变化规律，由此提供较为可靠的主要成形工艺参数。因此基于有限元法的塑性成形数值模拟技术是当前国际上极具发展潜力的成形技术前沿研究课题之一。 正确设计和控制金属塑性成形过程的前提条件是充分掌握金属流动、应力应变状态、热传导、润滑、加热与冷却及模具结构设计等方面的知识。任何分析方法都是为工程技术人员服务的，其目的是帮助工程技术人员掌握金属流动过程中应力应变状态等方面知识，一个好的分析方法至少应包括以下几个功能： （1）在未变形体（毛坯）与变形体（产品）之间建立运动学关系，预测金属塑性成形过程中的金属流动规律，其中包括应力应变场量变化、温度变化及热传导等。 （2）计算金属塑性成形极限，即保证金属材料在塑性变形过程中不产生任何表面及内部缺陷的最大变形量可能性。 （3）预测金属塑性成形过程得以顺利进行所需的成形力及能量，为正确选择加工设备和进行模具设计提供依据。 当前，有限元法已成为分析和研究金属塑性成形问题的最重要的数值分析方法之一，它具有以下优点：（1）由于单元形状具有多样性，有限元法使用与任何材料模型，任意的边界条件，任意的结构形状，在原则上一般不会发生处理上的困难。金属材料的塑性加工过程，均可以利用有限元法进行分析，而其它的数值

基于D_P准则的三维弹塑性有限元增量计算的有效算法

基于D-P准则的三维弹塑性有限元 增量计算的有效算法 A practical3D ela sto2pla stic incremental method in FEM ba sed on D-P yield criteria 杨　强,陈　新,周维垣 (清华大学水利系,北京　100084) 摘　要:针对岩土材料常用的D-P准则,提出了一种新的增量分析方法,不用形成弹塑性增量矩阵,直接导出了符合正交流动法则的转移应力的解析解。该方法无论是对小步长还是大步长加载均有良好的收敛性。当采用精细的步长划分时,它就是严格意义上的理想弹塑性增量计算。在大步长情况下,在收敛域内最大载荷低于结构真实的极限承载力;对应的应力场是一个静力容许应力场;同时由于正交流动法则在平均意义下得到满足,收敛域内最大载荷接近结构真实的极限承载力。按此法所得结果接近真解且偏于安全。将整个计算模型装入三维非线性有限元程序TFI NE中,对某拱坝进行了超载分析。 关键词:转移应力;极限载荷;点安全度 中图分类号:T U452 文献标识码:A 文章编号:1000-4548(2002)01-0016-05 作者简介:杨　强(1964-),男,云南人。1988年在清华大学获硕士学位,1996年在奥地利Innsbruck大学获博士学位,现为清华大学水利系教授。主要从事水工结构及岩石力学方面的研究工作。 Y ANG Qiang,CHE N X in,ZH OU Wei2yuan (Department of Hydraulic Engineering,Tsinghua University,Beijing100084,China) Abstract:In this paper,focused on popularly used D-P yield criteria in geomaterials,a new incremental method in which the stresses to be trans2 ferred according to normal flow rule are directly derived without forming elasto2plastic increment matrix,was proposed.This method converges for either small load steps or large load steps.When very small load steps are used,the method is equivalent to standard elasto2plastic incremental method.When large load steps are used,the maximum load applied is lower than limit load in structure,the calculated stress field is an static ad2 missible one.As normal flow rule is satisfied in average,the maximum load is close to limit load.The soltion calculated by the method is on the safe side and close to real solution.The method was embedded into a3D nonlinear FEM software named TFINE,and overloading analysis was performed on an arch dam. K ey words:the stresses to be transferred;normal flow rule;limit load 1　引 言Ξ 岩土材料具有很复杂的本构特性,如各向异性、硬化、软化等,目前描述岩土材料的本构模型非常多。但在实际工程三维有限元计算分析中,尤其是在岩体工程里,大量使用的仍是最简单D-P准则及理想弹塑性分析。其主要原因是参数选取不易。如在二滩高拱坝建设中,做了大量坝肩岩体现场大型抗剪试验,但具体到某一岩级,试验点数仍然很有限,且离散性很大。很难完全依赖试验确定参数,一般都要进行工程类比,对中、小工程工程类比更是参数确定的主要手段。最终一般只能给出岩体的抗剪参数f,c值。在这种情况下,从工程实用角度来说,追求本构关系的精致、完备并无太多实用意义。 相对而言,在岩土工程三维非线性有限元分析里,计算收敛性是一个较大的问题。弹塑性增量计算要采用精细的步长划分,才能确保计算收敛到正确解。在岩土工程,尤其是岩体工程里,荷载量级都很大,如高拱坝对水荷载的极限承载力可达上亿吨,而这对两岸高陡边坡所承受的的自重荷载来说,还只是一个小数,又如高地应力区大型地下洞室、高边坡(如三峡船闸高边坡)开挖过程中的释放荷载量级也十分巨大。若采用精细的步长划分,计算量将很大。岩体地质构造复杂,三维网格划分时经常会有畸形单元。由于地址缺陷或加固措施导致相邻单元材料性质差异过大,再加上高水平的荷载,各种因素交互影响,使得在计算过程中,经常出现局部发散现象,使得增量计算难以进行下去,最终结果可信度低,也难以从计算结果判断何时结构丧失稳定性。而对岩土工程来说,往往更关注结构的稳定性和极限承载力,而非应力和位移分布。 Ξ基金项目:国家自然科学基金资助项目(59879005);清华大学基础研究基金资助项目 收稿日期:2001-04-12 　第24卷　第1期岩　土　工　程　学　报V ol.24　N o.1 2002年 1月Chinese Journal of G eotechnical Engineering Jan.,　2002

弹塑性有限元法与刚塑性有限元法

弹塑性有限元法与刚塑性有限元法 板料成形数值模拟涉及到连续介质力学中材料非线性、几何非线性、边界条件非线性三非线性问题的计算，难度很大。随着非线性连续介质力学理论、有限元方法和计算机技术的发展，通过高精度的数值计算来模拟板料成形过程已成为可能。从70年代后期开始，经过近二十年的发展，板料成形数值模拟逐渐走向成熟，并开始在汽车、飞机等工业领域得到实际应用。 本文评述了板料成形数值模拟的发展历史和最新进展，并指出了该领域的发展趋势。 1、板料成形的典型成形过程、物理过程与力学模型 典型成形过程 板料成形的具体过程多种多样，在模拟分析时，可归纳成如图1所示的典型成形过程。成形时，冲头在压力机的作用下向下运动，给板料一个作用压力，板料因此产生运动与变形。同时，冲头、压力圈和凹模按一定方式共同约束板料的运动与变形，从而获得所要求的形状与尺寸。 物理过程 板料成形的物理过程包括模具与板料间的接触与摩擦；由于金属的塑性变形而导致的加工硬化和各向异性化；加工中可能产生的皱曲、微裂纹与破裂及由于卸载而在零件中产生回弹。 力学模型 板料成形过程可归纳成如下的力学问题：

给定冲头位移、凹模位移及压边圈历程函数，求出板料的位移历程函数，使其满足运动方程、初始条件、边界条件、本构关系及接触摩擦条件。 2板料成形数值模拟的发展历史 塑性有限元方法的发展 根据材料的本构关系，用于板料成形分析的非线性有限元法大体上分为刚-（粘）塑性与弹-（粘）塑性两类。 粘塑性有限元法很早就在板料成形分析中应用过，只是未能推广。事实上，粘塑性有限元法适用于热加工。在热加工时，应变硬化效应不显著，材料形变对变形速率有较大敏感性。

金属塑性成形综述

金属塑性成形 摘要：金属塑性成形技术是机械冶金、汽车拖拉机、电工仪表、宇航军工、五金日用品等制造业最基本，最古老，亦是极重要的加工手段之一。文章主要对塑性成形的基本方法、主要研究内容，发展趋势做了综合介绍。 一、引言 塑性成形技术具有高产、优质、低耗等显著特点,已成为当今先进制造技术的重要发展方向。据国际生产技术协会预测,21世纪,机械制造工业零件粗加工的75%和精加工的50%都采用塑性成形的方式实现。【1】 在现代制造技术中，人们广泛的利用金属材料生产各种零件和产品。金属加工方法多种多样，包括成型、切削等。金属塑性成形是其中一种重要的加工方法，是利用金属在外力作用下产生的塑性变形来获得具有一定形状、尺寸和力学性能的原材料、毛坯或零件的生产方法，因此也称为金属塑性加工或金属压力加工。 图1 传统金属塑性成形工艺 二、金属塑性成形的主要形式 金属塑性成形工艺的种类有很多，包括轧制、挤压、拉拔、锻造和冲压等基本工艺类型。随着技术的发展，也有很多新的成型方式出现，它们具备精密、高效、节能、节材、清洁等优点，得到广泛关注。

2.1 体积成型 金属体积成型是指对金属块料、棒料或厚板在高温或室温下进行成形加工的方法，主要分为热态金属体积成型和冷温态金属体积成型。热态金属变形过程可分为热锻、轧制、挤压、拉拔、辗压等工艺技术；冷温态变形过程可分为冷锻、冷精轧、冷挤压、冷拔、冷辗扩等工艺。 2.2 板材成型 所谓板材成型是指用板材、薄壁管、薄型材等作为原材料进行塑性加工的成形方法。在忽略板厚的变化时，可视为平面变形问题来处理，板材成型可分为:冲裁、弯曲、拉延、胀形、翻边、扩孔、辊压等工艺技术。 2.3 粉末态金属成形 随着制粉技术的发展，其应用领域不断扩展，对于复杂形状的机械零件来说，它具有高效、精密成形的特点，但成本较高，机械性能不如整体金属材料。粉末态金属成形的工艺过程为制粉、造型、压实、烧结、精锻。 2.4半固态金属材料成形 70年代开发研究的新技术，原金属材料作过特殊前处理，当材料加热到一定温度时可使30%的金属材料处于融溶状态，其余70%的金属材料呈均匀细颗粒组织的固态。在此状态加压变形，其流动性特好，可成形结构形状特别复杂的零件，而变形杭力很小。 2.5 复合成形技术 现代的科学越来越相互交叉、渗透，出现许多边缘学科、交叉学科一样，材料成形技术也逐渐突破原有铸、锻、焊、粉末冶金等技术相互独立的格局，相互融合、渗透，产生了种类繁多的“复合成形技术”。【2】金属塑性的复合成型技术主要有两个方面 （1）各种成形工艺的组合优化达到优化工艺和产品的目的。 （2）铸、锻、焊、热处理等不同加工方法的组合。 三、金属塑性成形技术主要研究内容 由于压力加工中,少、无切屑的特点和精密加工技术的发展,使金属塑性成型理论的研究受到日益广泛的重视而进入工程应用的前列.一般认为,研究金属塑性科学的历史开始于Tresa在1864年提出的屈服准则,至今不过100多年,而

有限元法的发展及在塑性加工研究中的应用

有限元法的发展及在塑性加工研究中的应用 金属成形过程的研究方法大致可分为三类[33]：第一类是基于经典塑性理论的解析方法，其中包括：①精确地联立求解塑性理论基本方程的数学解析法。②将平衡方程和塑性条件简化后联立求解塑性理论基本方程的主应力法。③针对平面问题提出的滑移线法。④基于能量守恒原理的能量法和上限法；第二类是以实验数据为分析基础的实验力学研究方法，如视塑性法、网格法、密栅云纹法等；第三类是随着塑性理论和计算机应用的发展，由传统的方法演化出来的数值方法，主要包括：①上限元法(UBET)，②矩阵算子法，③有限差分法(FDM)，④加权余数法(WRM)，⑤边界元法(BEM)，⑥有限元法(FEM)等。 近年来，随着计算机技术的发展和普遍应用，对塑性成形问题的求解起到了很大的促进作用。金属成形过程伴随着很大的塑性变形，既存在材料非线性（应力与应变之间的非线性），又有几何非线性（应变与位移之间的非线性），变形机制十分复杂。加上复杂的边界条件以及数学上的困难，使得成形过程的求解十分复杂。有限元法是20世纪60年代提出的一种分析方法，已能在考虑变形热效应以及工件与模具、周围介质热交换的情况下，确定变形体内的应力、应变和温度情况。成形过程的设计一般依赖于图表、经验公式等传统的经验设计，这对于工艺设计复杂而工艺参数控制要求严格的精确塑性成形而言，其浪费会更大，有限元数值模拟方法在计算机上再现了零件的制造过程，使人们能够直观、全面地了解成形过程。精确塑性成形需要高质量的控形和控性来保证，通过模拟可以在零件实际制造前，发现成形条件和工艺中可能存在的问题及缺陷，优化成形工艺，提高生产效率，降低生产成本，是一种非常有效的研究方法。有限元法所以能获得如此广泛应用，与计算机技术的发展、特别是模拟的实践指导与应用密切相关。有限元模拟有其独特的优越性，通过模拟计算可以作为试验设计依据，从而简化试验，降低试验成本、避免多次试验的尝试和反复修改。计算机数值模拟技术已成为研究和发展先进塑性成形理论和技术的强有力工具，为高质量、低成本、短周期地实现塑性成形产品的开发创造了有利条件。 有限元法是对塑性成型过程进行数值模拟的最有效方法。可以比较精确地求解变形过程中的各种场变量，如速度（位移）场、应变场和应力场等，从而为变形工艺分析提供科学依据和具体实践的指导。 Marcal和King[34]于1967年最早提出了弹塑性有限元，它同时考虑弹性变形和塑性变形，弹性区采用Hook定律，塑性区采用Prandtl-Reuss方程和Mises屈服准则。1973年S.Kobayashi和C.H.Lee[35]针对弹塑性有限元法存在的问题提出了所谓“矩阵法”的刚塑性有限元法，用来分析金属塑性成形问题。这类有限元法不计弹性变形，采用Levy-Mises 方程作为本构方程，满足体积不变条件，并采用率方程描述，变形后物体的形状通过在离散空间上对速度积分而获得，从而避开了有限变形中几何非线性问题。近年来，刚塑性有限元法己被广泛应用并解决了许多金属塑性成形问题。与此同时，针对刚塑性有限元法在求解过程中存在的一些问题，提出了各种相应的解决方法，如处理材料不可压缩条件的Lagrange乘子法[35]、罚函数法[36]和体积可压缩法[37]等。O.C.Zienkiewicz[38]等将刚塑性材料看作是非牛顿不可压缩粘性流体，导出了刚粘塑性有限元求解列式，并于1972年提出了刚粘塑性有限元法。 针对以上提出的模型，后人根据具体材料的不同情况进行了大量的模拟试验，使这些有限元方法得到了广泛的应用，同时结合实际对模型进行修正，促进有限元模拟在生产的指导作用[39~43]。N.L.Dung等[39]采用刚塑性有限元法，对涡轮叶片精锻过程进行了二维有限元模拟，在模拟过程中采用三角形和四边形的混合单元，并采用Lagrange乘子法处理体积不可压缩条件，对金属在对称和不对称的上下模的作用下的流动情况进行了对比分析。王梦寒等[42]用刚粘塑性有限元法对油泵定了温挤成形过程的模拟分析，优化了成形工艺及模型参数，通过与物理试验比较，验证了数值模拟分析的准确性。 对钛合金塑性成形过程进行数值模拟是实现计算机仿真的基础。由于有限元法可以全面地考虑变形过程中材料的动态特性、各种边界条件和初始条件的影响，即使对于复杂边界仍可达到满意的模拟精度。因此，目前对塑性成形过程进行数值模拟的方法主要是有限元法。 有限元法的一般解题步骤是： （1）连续体的离散化。把求解的连续体离散成有限数目的单元，单元的类型有多种，二维中有三边形、四边形，三维中有四面体、六面体等。模拟中根据实际模拟对象和模拟的需要，选择合理的单元类型、大小和数目。 （2）选择插值函数。选择满足某些需要（单元连续性、边界协调性等）的、联系单元节点和单元内部各点位移（速度）的插值函数，保证计算结果的精确性。插值函数通常为多项式，以便于微分和积分。 （3）建立单元的刚度矩阵或能量泛函。根据变分原理，对弹性和弹塑性有限元，推导单元的刚度矩阵[K]e，用此矩阵把单元节点位移{u}e和节点力{P}e联系起来,即e [K]{u}e={P}e(1-1)对刚塑性有限元，建立以节点位移为自变函数的单元能量泛函 （4）建立整体方程。对弹性和弹塑性有限元，这个过程包括由各单元的刚度矩阵集合成整个变形体的总刚度矩阵[K]，以及由单元节点力列阵集合成的总载荷列阵{P}，从而建立表示整个变形体的节点位移和总载荷关系的方程组，即[K]{u}={P}( 对于刚塑性有限元，则建立整个变形体的能量泛函变分方程组，

金属塑性

1、什么是金属塑性？什么是塑性成型？塑性成型有何特点？ 塑性：在外力作用下使金属材料发生塑性变形而不破坏其完整性的能力称为塑性。利用金属在一定的外力作用下产生塑性变形，并获得具有一定形状、尺寸和机械性能的材料、毛坯或零件的加工方法，称为金属的塑性成形（也称压力加工）。塑性成型特点：1组织、性能好2材料利用率高3尺寸精度高4生产率高，易实现连续化、自动化、高速、大批量生产 不足：设备较庞大，相对能耗较高，成本较高2试述塑性成型的一般分类？ 一、板料成型：1、一次加工：1）轧制2）挤压3）拉拔2、二次加工：1）自由锻2)模锻二、块料成型：1、分离工序:1)冲裁2）落料2、成型工序：1）弯曲2）拉深三、按温度分：热成型、冷成型、温成型3、试简述滑移和孪生两种变形机理的主要区别？滑移与孪生的比较相同点：通过位错运动实现；两者都不改变晶体结构类型区别：1）晶体中的位向滑移：晶体中已滑移部分与未滑移部分的位向相同孪生：已孪生部分（孪晶）和未孪生部分(基体)的位向不同，两部分之间具有特定的位向关系（镜面对称）2）变形机制：滑移是全位错运动的结果；孪生是部分位错3）对塑性变形的贡献：总变形量大；孪生（小）4）变形应力：近似临界分切应力；高于临界分切应力5）变形条件：一般情况下，先发生滑移变形；滑移变形难以进行时，或晶体对称度很低、变形温度较低、加载速率较高，发生孪生变形4、试分析多晶体塑性变形的特点？（1）各晶粒变形的不同时性首先在位向有利、滑移系上切应力分量已优先达到临界值的晶粒内发生2）各晶粒变形的相互协调性晶粒的变形需要相互协调配合，才能保持晶粒之间的连续性，即变形不是孤立和任意的。（3）变形的不均匀性软位向的晶粒先变形，硬位向的晶粒后变形，其结果必然是各晶粒变形量的差异，这是由多晶体的结构特点所决定的。5、什么是加工硬化？加工硬化产生的原因？加工硬化对塑性加工有何利弊？1）加工硬化：塑性变形时，随着内部组织结构变化，金属金属强度、硬度增加，而塑性、韧性降低的现象。2）加工硬化是位错与交互作用有关，随着塑性变形的进行，位错密度不断增加，位错反应和相互交割加剧，结果产生固定割阶、位错纠缠等障碍。以致形成细胞亚状结构，是位错难以越过这些障碍而被限制在一定的范围内运动。金属要继续变形，就要不断外力，才能克服强大的交互作用。3）有利的方面：1、是金属强化的重要途径2、对不能用热处理方法强化的材料，借助冷塑性变形来提高其力学性能。3、对改善板料成型性能有积极的意义。不利的一面：金属塑性下降、变形抗力升高、继续变形越来越困难；对高硬化速率的多道次成形，需增加中间退火来消除加工硬化，降低了生产效率、提高成本6、什么是动态再结晶？影响动态再结晶的因素有哪些？1）动态在结晶：是在热塑性变形过程发生的再结晶。2）影响因素：位错能的高低，晶界迁移的难易程度、应变速率、变形温度等有关。7、什么是动态回复？为什么说动态回复是热塑变形的主要软化机制？动态回复：在热塑性变形过程中发生的回复。原因: 层错能高，变形时扩展位错的宽度窄、集束容易，位错的交滑移和攀移容易进行，位错容易在滑移面间转移，而使异号位错互相抵消，结果使位错密度下降，畸变能降低，不足以达到动态再结晶所需的能量水平8、与常规塑性变形相比，超塑性变形具有哪些主要特征？1、大伸长率，高达百分之几千2、无缩颈，拉伸时变现均匀的截面缩小，断面收缩率甚至可接近100%3、低流动应力，仅（几个—几十个）N/mm，对应变速率非常敏感4、具有极好的流动性和充填性，加工复杂精确的零件。9、什么是细晶超塑性？什么是相变超塑性？细晶超塑性：一定恒温，应变速率和晶粒度都满足要求，呈现的超塑性相变超塑性：要求具有相变或同素异构转变，一定的外力下，金属或合金在相变温度附近反复加热和冷却，经过一定的循环次数后，获得很大的伸长率。10、超塑性变形的力学方程中的m的物理意义是什么？m值指的是应变速率敏感指数：反应材料抗局部收缩或产生均匀拉深变形的能力。 11、什么是塑性？什么是塑性指标？为什么说塑性指标只具有相对意义？1）塑性：是金属在外力作用下产生永久变形而不破坏其完整性的能力。2）塑性指标：金属在破坏前产生的

有限元法小结

Elements in ABAQUS

A1.2
Elements in ABAQUS
? ABAQUS单元库中提供广泛的单元类型，适应不同的结构和几何特征 The wide range of elements in the ABAQUS element library provides flexibility in modeling different geometries and structures. – Each element can be characterized by considering the following: 单元特性： ? Family 单元类型 ? Number of nodes 节点数 ? Degrees of freedom 自由度数 ? Formulation 公式 ? Integration 积分

A1.3
Elements in ABAQUS
?单元类型（Family） –A family of finite elements is the broadest category used to classify elements. –同类型单元有很多相 同的基本特征。 Elements in the same family share many basic features. –同种类单元又有很多 变化：There are many variations within a family.
continuum (solid elements)
shell elements
beam elements
rigid elements
membrane elements
infinite elements
special-purpose elements like springs, dashpots, and masses
truss elements

面心纯铝轧制织构的晶体塑性有限元模拟_皮华春

面心纯铝轧制织构的晶体塑性有限元模拟 皮华春1) 韩静涛1) 章传国1) A.K.Tieu 2) 姜正义2) 1)北京科技大学材料科学与工程学院,北京100083 2)伍伦贡大学机械材料与机电学院,伍伦贡NSW 2522,澳大利亚 摘 要 基于率相关晶体塑性本构模型,分别将Taylor 模型和有限单元模型两种多晶模型嵌入大型有限元程序ABAQUS,实现了晶体塑性学有限元模拟.直接将电子背散射衍射(EBSD)获取的晶粒初始取向输入晶体塑性有限元模型,预测了两种不同应变情况下面心1050纯铝轧制织构的演化.模拟结果与EBSD 实验测得的织构演化结果有较好的一致性,随着变形程度的增加,预测织构与实测织构变得更加锋锐.经过比较,T aylor 型模型预测出了{4411}3111184的Dillamore 取向,而有限单元模型预测出了铜型织构取向,比T aylor 模型预测结果更接近实验验证结果.两种模型并不能预测出{011}32114黄铜取向、{123}35234S 取向、{011}31004Goss 取向及其他理想取向.关键词 轧制织构;晶体塑性;有限元;电子背散射衍射(EBSD)分类号 TG302 收稿日期:2006204217 修回日期:2006206226 基金项目:澳大利亚研究院(ARC)国际合作项目(No.DP0451197)作者简介:皮华春(1978)),男,博士研究生;韩静涛(1957)),男,教授,博士生导师 大多数工程使用材料为多晶材料,并具有择优取向也就是织构.织构不仅与材料的性能相关,而且能够揭示材料的变形历史[1].研究人员很早就开始了对变形织构的模拟[2-4].近年来,他们中的一部分人将晶体塑性理论与力学界经常运用的有限元法有机结合在一起,在介观尺度形成了一种新的织构模拟方法)))晶体塑性有限元法(CPFEM).作为一个强大的模拟工具,晶体塑性有限元法不仅用于模拟变形织构的演化,而且能够用于与材料织构相关的各种性能响应[5-11]. 通常由X 射线测定的晶体学织构是从统计学的角度出发观察多晶体取向分布状况,因而是从宏观的角度分析问题.宏观织构往往不能确定微观结构和晶粒的单个取向.以往的晶体塑性模拟一般是将X 射线获得的织构极图通过离散化方法得到表征该材料晶粒取向的欧拉角,再将这些数据输入模拟程序中.近年来,人们在扫描电子分析技术基础上开发出了背散射电子衍射分析技术(EBSD),该技术可以在观测微观组织结构的同时快速、统计性地获取多晶体各个晶粒的取向信息[12] .这就使得晶体塑性模拟能够直接运用EBSD 获得的晶粒取向进行模拟与验证,从而不仅简化了晶体塑性模拟过程而且保证了模拟结果的准确性. 本文直接运用EBSD 获取的晶粒取向数据通过率相关多晶体塑性模型模拟了面心纯铝轧制织构的 演化,并通过EBSD 实验验证了模拟的结果. 1 实验过程 具有工业纯度的1050纯铝热轧后取样,化学成分(质量分数)为:Al 9917%,Fe 0116%,Si 0104%,Ca 0101%,Mg 0106%,Cu 0103%.采用毛卫民[13]和Raabe 等[14]的方法获得晶粒细小的初始任意织构试样.首先将纯铝试样在三个相互垂直的方向进行三次循环锻造,三次锻造变形量分别为25%,15%和5%依次降低.然后对锻后样品进行500e 、015h 的退火.冷轧实验在北京科技大学冷轧实验中心的二辊轧机上进行,轧辊直径为250mm ,试样原始尺寸为40mm @25mm @20mm (RD @TD @ND ),以煤油作为润滑剂.经过5道次轧制第1个试样的最终厚度为617mm ,计算压下量为67%(真应变约111);经过12道次轧制第2个试样的最终厚度为112mm ,计算压下量为94%(真应变约218).初始取向和轧后织构的EBSD 结果是在北京工业大学的Jeol JSM 6500F 场发射扫描电镜上测得,该设备的加速电压015~30kV ,分辨率115nm .为了得到清晰的菊池线,EBSD 试样经过了打磨、机械抛光和电解抛光.观察面为试样在厚度方向的中心面. 2 晶体塑性有限元模型 晶体塑性理论源于Taylor [15] 开创性的工作,在 他的研究工作中引入了滑移和晶格转动的思想.随后,H ill 和Rice [16]、Asaro [17]在此基础上给出了一套 第29卷第9期2007年9月北京科技大学学报 Jour nal of Univer sity of Science and Technology Beijing Vol.29No.9Sep.2007