求VUMAT 塑性各向同性强化模型

小侯

求VUMAT塑性各向同性强化模型  感谢了 急用！！！    小弟对塑性力学有了一些基本认识  愿意和大侠们讨论  向大侠们学习

RathLusia

here is the VUMAT for Isotropic Hardening:

1. C
   
2. C Coding for Isotropic Hardening Plasticity VUMAT
   
3. C
   
4. 
   
5.        SUBROUTINE VUMAT(
   
6. C Read only -
   
7.      1 NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, LANNEAL,
   
8.      2 STEPTIME, TOTALTIME, DT, CMNAME, COORDMP, CHARLENGTH,
   
9.      3 PROPS, DENSITY, STRAININC, RELSPININC,
   
10.      4 TEMPOLD, STRETCHOLD, DEFGRADOLD, FIELDOLD,
    
11.      5 STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,
    
12.      6 TEMPNEW, STRETCHNEW, DEFGRADNEW, FIELDNEW,
    
13. C Write only -
    
14.      7 STRESSNEW, STATENEW, ENERINTERNNEW, ENERINELASNEW)
    
15. C
    
16.        INCLUDE 'VABA_PARAM.INC'
    
17. C
    
18.        DIMENSION PROPS(NPROPS), DENSITY(NBLOCK), COORDMP(NBLOCK),
    
19.      1 CHARLENGTH(NBLOCK), STRAININC(NBLOCK, NDIR+NSHR),
    
20.      2 RELSPININC(NBLOCK, NSHR), TEMPOLD(NBLOCK),
    
21.      3 STRETCHOLD(NBLOCK, NDIR+NSHR),DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR),
    
22.      4 FIELDOLD(NBLOCK, NFIELDV), STRESSOLD(NBLOCK, NDIR+NSHR),
    
23.      5 STATEOLD(NBLOCK, NSTATEV), ENERINTERNOLD(NBLOCK),
    
24.      6 ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),
    
25.      7 STRETCHNEW(NBLOCK, NDIR+NSHR),DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),
    
26.      8 FIELDNEW(NBLOCK, NFIELDV), STRESSNEW(NBLOCK,NDIR+NSHR),
    
27.      9 STATENEW(NBLOCK, NSTATEV), ENERINTERNNEW(NBLOCK),
    
28.      1 ENERINELASNEW(NBLOCK)
    
29. 
    
30.        CHARACTER*80 CMNAME
    
31. 
    
32. C LOCAL ARRAYS
    
33. C ----------------------------------------------------------------
    
34. C EELAS - ELASTIC STRAINS
    
35. C EPLAS - PLASTIC STRAINS
    
36. C FLOW - DIRECTION OF PLASTIC FLOW
    
37. C ----------------------------------------------------------------
    
38. C
    
39.        DIMENSION EELAS(6),EPLAS(6),FLOW(6), HARD(3)
    
40. C
    
41.        PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, SIX=6.D0,
    
42.      1 ENUMAX=.4999D0, NEWTON=10, TOLER=1.0D-6,HALF=0.5D0)
    
43. 
    
44. C ----------------------------------------------------------------
    
45. C UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC MISES PLASTICITY
    
46. C CANNOT BE USED FOR PLANE STRESS
    
47. C ----------------------------------------------------------------
    
48. C PROPS(1) - E
    
49. C PROPS(2) - NU
    
50. C PROPS(3..) - SYIELD AN HARDENING DATA
    
51. C CALLS UHARD FOR CURVE OF YIELD STRESS VS. PLASTIC STRAIN
    
52. C ----------------------------------------------------------------
    
53. C
    
54. C ELASTIC PROPERTIES
    
55. C
    
56.        EMOD=PROPS(1)
    
57.        ENU=MIN(PROPS(2), ENUMAX)
    
58.        EBULK3=EMOD/(ONE-TWO*ENU)
    
59.        EG2=EMOD/(ONE+ENU)
    
60.        EG=EG2/TWO
    
61.        EG3=THREE*EG
    
62.        ELAM=(EBULK3-EG2)/THREE
    
63.        NVALUE=NPROPS/2-1
    
64. C
    
65. C CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN
    
66. C
    
67.        IF ( STEPTIME .EQ. ZERO ) THEN
    
68.        DO K = 1, NBLOCK
    
69.        TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
    
70.        STRESSNEW(K,1) = STRESSOLD(K,1)
    
71.      *                  + EG2 * STRAININC(K,1) + ELAM * TRACE
    
72.        STRESSNEW(K,2) = STRESSOLD(K,2)
    
73.      *                  + EG2 * STRAININC(K,2) + ELAM * TRACE
    
74.        STRESSNEW(K,3) = STRESSOLD(K,3)
    
75.      *                  + EG2 * STRAININC(K,3) + ELAM * TRACE
    
76.        STRESSNEW(K,4)=STRESSOLD(K,4) + EG2 * STRAININC(K,4)
    
77.        IF ( NSHR .GT. 1 ) THEN
    
78.        STRESSNEW(K,5)=STRESSOLD(K,5) + EG2 * STRAININC(K,5)
    
79.        STRESSNEW(K,6)=STRESSOLD(K,6) + EG2 * STRAININC(K,6)
    
80.        END IF
    
81.        END DO
    
82.        ELSE
    
83. C
    
84. C111111111111
    
85. C
    
86.        DO  K = 1, NBLOCK
    
87. C
    
88.        TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
    
89.        STRESSNEW(K,1) = STRESSOLD(K,1)
    
90.      *                  + EG2 * STRAININC(K,1) + ELAM * TRACE
    
91.        STRESSNEW(K,2) = STRESSOLD(K,2)
    
92.      *                  + EG2 * STRAININC(K,2) + ELAM * TRACE
    
93.        STRESSNEW(K,3) = STRESSOLD(K,3)
    
94.      *                  + EG2 * STRAININC(K,3) + ELAM * TRACE
    
95.        STRESSNEW(K,4)=STRESSOLD(K,4) + EG2 * STRAININC(K,4)
    
96.        IF ( NSHR .GT. 1 ) THEN
    
97.        STRESSNEW(K,5)=STRESSOLD(K,5) + EG2 * STRAININC(K,5)
    
98.        STRESSNEW(K,6)=STRESSOLD(K,6) + EG2 * STRAININC(K,6)
    
99.        END IF
    
100. C
     
101.        SMISES=(STRESSNEW(K,1)-STRESSNEW(K,2))**2 +
     
102.      1        (STRESSNEW(K,2)-STRESSNEW(K,3))**2 +
     
103.      2        (STRESSNEW(K,3)-STRESSNEW(K,1))**2
     
104.        DO 90 K1=NDIR+1,NDIR+NSHR
     
105.          SMISES=SMISES+SIX*STRESSNEW(K,K1)**2
     
106. 90     CONTINUE
     
107.        SMISES=SQRT(SMISES/TWO)
     
108. C
     
109. C RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE FORWARD
     
110. C ALSO RECOVER EQUIVALENT PLASTIC STRAIN
     
111. C
     
112. C      CALL ROTSIG(STATEOLD(K,1), DROT, EELAS, 2, NDIR, NSHR)   
     
113. C      CALL ROTSIG(STATEOLD(K,1+NDIR+NSHR), DROT, EPLAS, 2, NDIR, NSHR)
     
114. 
     
115.        DO K1=1, NDIR+NSHR
     
116.        EELAS(K1)=STATEOLD(K,K1)
     
117.        EPLAS(K1)=STATEOLD(K,K1+NDIR+NSHR)
     
118.        END DO
     
119.        EQPLAS=STATEOLD(K,1+2*(NDIR+NSHR))!!!
     
120. C
     
121. C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
     
122. C
     
123. 
     
124.        CALL VUHARD(SYIEL0, HARD, EQPLAS, EQPLASRT,TOTALTIME,DT,TEMP,
     
125.      1             DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,
     
126.      2             NUMFIELDV,NVALUE,PROPS(3))
     
127. C
     
128. C DETERMINE IF ACTIVELY YIELDING
     
129. C
     
130.        IF (SMISES.GT.(ONE+TOLER)*SYIEL0) THEN
     
131. C
     
132. C ACTIVELY YIELDING
     
133. C SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS
     
134. C CALCULATE THE FLOW DIRECTION
     
135. C
     
136.        SHYDRO=(STRESSNEW(K,1)+STRESSNEW(K,2)+STRESSNEW(K,3))/THREE
     
137.        DO K1=1,NDIR
     
138.        FLOW(K1)=(STRESSNEW(K,K1)-SHYDRO)/SMISES
     
139.        END DO
     
140.        DO K1=NDIR+1,NDIR+NSHR
     
141.        FLOW(K1)=STRESSNEW(K,K1)/SMISES
     
142.        END DO
     
143. C
     
144. C SOLVE FOR EQUIVALENT VON MISES STRESS
     
145. C AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATION
     
146. C
     
147.        SYIELD=SYIEL0
     
148.        DEQPL=ZERO
     
149.        DO KEWTON=1, NEWTON
     
150.        RHS=SMISES-EG3*DEQPL-SYIELD
     
151.        DEQPL=DEQPL+RHS/(EG3+HARD(1))
     
152.        CALL VUHARD(SYIELD,HARD,EQPLAS+DEQPL,EQPLASRT,TOTALTIME,DT,TEMP,
     
153.      1            DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NUMFIELDV,
     
154.      2            NVALUE,PROPS(3))
     
155.        IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10
     
156.        END DO
     
157. C
     
158. C WRITE WARNING MESSAGE TO THE .MSG FILE
     
159. C
     
160.        WRITE(7,2) NEWTON
     
161.     2  FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT',
     
162.      1 'CONVERGE AFTER ',I3,' ITERATIONS')
     
163.    10  CONTINUE
     
164. C
     
165. C UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND
     
166. C EQUIVALENT PLASTIC STRAIN
     
167. C
     
168.        DO K1=1,NDIR
     
169.        STRESSNEW(K,K1)=FLOW(K1)*SYIELD+SHYDRO
     
170.        EPLAS(K1)=EPLAS(K1)+THREE/TWO*FLOW(K1)*DEQPL
     
171.        EELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPL
     
172.        END DO
     
173.        DO K1=NDIR+1,NDIR+NSHR
     
174.        STRESSNEW(K,K1)=FLOW(K1)*SYIELD
     
175.        EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL
     
176.        EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL
     
177.        END DO
     
178.        EQPLAS=EQPLAS+DEQPL
     
179. C
     
180. C CALCULATE PLASTIC DISSIPATION
     
181. C
     
182.        SPD=DEQPL*(SYIEL0+SYIELD)/TWO
     
183.          ENDIF
     
184. C
     
185. C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS
     
186. C IN STATE VARIABLE ARRAY
     
187. C
     
188.        DO K1=1, NDIR+NSHR
     
189.        STATENEW(K,K1)=EELAS(K1)
     
190.        STATENEW(K,K1+NDIR+NSHR)=EPLAS(K1)
     
191.        END DO
     
192.        STATENEW(K,1+2*(NDIR+NSHR)) = STATEOLD(K,1+2*(NDIR+NSHR)) + DEQPL
     
193. C
     
194. C UPDATE THE SPECIFIC INTERNAL ENERGY -
     
195. C
     
196.        IF ( NSHR .EQ. 1 ) THEN
     
197.        STRESSPOWER = HALF * (
     
198.      * ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
     
199.      * ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
     
200.      * ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3) ) +
     
201.      * ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4)
     
202.        ELSE
     
203.        STRESSPOWER = HALF * (
     
204.      * ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
     
205.      * ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
     
206.      * ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3) ) +
     
207.      * ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4) +
     
208.      * ( STRESSOLD(K,5) + STRESSNEW(K,5) ) * STRAININC(K,5) +
     
209.      * ( STRESSOLD(K,6) + STRESSNEW(K,6) ) * STRAININC(K,6)
     
210.        END IF
     
211.        ENERINTERNNEW(K) = ENERINTERNOLD(K) + STRESSPOWER / DENSITY(K)
     
212. C
     
213. C UPDATE THE DISSIPATED INELASTIC SPECIFIC ENERGY -
     
214. C
     
215.        PLASTICWORKINC = HALF * ( SYIEL0 + SYIELD ) * DEQPL
     
216.        ENERINELASNEW(K) = ENERINELASOLD(K) + PLASTICWORKINC / DENSITY(K)
     
217.        ENDDO   
     
218.        END IF
     
219. C
     
220.        RETURN
     
221.        END
     
222. 
     
223.        SUBROUTINE VUHARD(SYIELD,HARD,EQPLAS,EQPLASRT,TIME,DTIME,TEMP,
     
224.      1                   DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,
     
225.      2                   NUMFIELDV, NVALUE,TABLE)
     
226.        INCLUDE 'VABA_PARAM.INC'
     
227.        CHARACTER*80 CMNAME
     
228.        DIMENSION HARD(3)
     
229. C
     
230.        DIMENSION TABLE(2, NVALUE)
     
231. C
     
232.        PARAMETER(ZERO=0.D0)
     
233. 
     
234. C
     
235. C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
     
236. C
     
237.        SYIELD=TABLE(1, NVALUE)
     
238.        HARD(1)=ZERO
     
239. C
     
240. C IF MORE THAN ONE ENTRY, SEARCH TABLE
     
241. C
     
242.        IF(NVALUE.GT.1) THEN
     
243.        DO K1=1, NVALUE-1
     
244.        EQPL1=TABLE(2,K1+1)
     
245.        IF(EQPLAS.LT.EQPL1) THEN
     
246.        EQPL0=TABLE(2, K1)
     
247.        IF(EQPL1.LE.EQPL0) THEN
     
248.        WRITE(7, 1)
     
249.     1  FORMAT(//, 30X, '***ERROR - PLASTIC STRAIN MUST BE ',
     
250.      1 'ENTERED IN ASCENDING ORDER')
     
251.        CALL EXIT
     
252.        ENDIF
     
253. C
     
254. C CURRENT YIELD STRESS AND HARDENING
     
255. C
     
256.        DEQPL=EQPL1-EQPL0
     
257.        SYIEL0=TABLE(1, K1)
     
258.        SYIEL1=TABLE(1, K1+1)
     
259.        DSYIEL=SYIEL1-SYIEL0
     
260.        HARD(1)=DSYIEL/DEQPL
     
261.        SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD(1)
     
262.        GOTO 10
     
263.        ENDIF
     
264.        END DO
     
265.    10  CONTINUE
     
266.        ENDIF
     
267.        RETURN
     
268.        END

复制代码

小侯

RathLusia 发表于 2013-11-20 14:25
here is the VUMAT for Isotropic Hardening:

这个你验证过吗？

787613786

收了，走过路过

chijiuziyou

RathLusia 发表于 2013-11-20 14:25
here is the VUMAT for Isotropic Hardening:

这个你验证过吗

mecookie1990

走过路过顺手牵过

鸿雁河上

感谢楼主

mengkp

RathLusia 发表于 2013-11-20 14:25
here is the VUMAT for Isotropic Hardening:

同问，验证过吗？

mengkp

RathLusia 发表于 2013-11-20 14:25
here is the VUMAT for Isotropic Hardening:

已验证，感谢大神无私的教诲

邂逅

顺手拿走

alis

FSDSFSFSSFSFSFSF

枯藤老树昏鸦

66666666666666666666