CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
incplas.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine incplas (elconloc, plconloc, xstate, xstateini, elas, emec, ithermal, icmd, beta, stre, vj, kode, ielas, amat, t1l, dtime, time, ttime, iel, iint, nstate_, mi, eloc, pgauss, nmethod, pnewdt, depvisc)
 

Function/Subroutine Documentation

◆ incplas()

subroutine incplas ( real*8, dimension(21)  elconloc,
real*8, dimension(802)  plconloc,
real*8, dimension(nstate_,mi(1),*)  xstate,
real*8, dimension(nstate_,mi(1),*)  xstateini,
real*8, dimension(21)  elas,
real*8, dimension(6)  emec,
integer  ithermal,
integer  icmd,
real*8, dimension(6)  beta,
real*8, dimension(6)  stre,
real*8  vj,
integer  kode,
integer  ielas,
character*80  amat,
real*8  t1l,
real*8  dtime,
real*8  time,
real*8  ttime,
integer  iel,
integer  iint,
integer  nstate_,
integer, dimension(*)  mi,
real*8, dimension(6)  eloc,
real*8, dimension(3)  pgauss,
integer  nmethod,
real*8  pnewdt,
real*8  depvisc 
)
23 !
24 ! calculates stiffness and stresses for the incremental plasticity
25 ! material law (Ref: J.C. Simo, A framework for finite strain
26 ! elastoplasticity, Comp. Meth. Appl. Mech. Engng., 66(1988)199-219
27 ! and 68(1988)1-31)
28 !
29 ! icmd=3: calculates stress at mechanical strain
30 ! else: calculates stress at mechanical strain and the stiffness
31 ! matrix
32 !
33 ! the stresses in the routine proposed by Simo are Kirchhoff
34 ! stresses. Since the stress in the hardening laws are Chauchy
35 ! stresses, they are converted into Kirchhoff stress by
36 ! multiplication with the Jacobian determinant
37 !
38  implicit none
39 !
40  character*80 amat
41 !
42  integer ithermal,icmd,i,j,kode,ivisco,ielastic,
43  & niso,nkin,ielas,iel,iint,nstate_,mi(*),id,leximp,lend,layer,
44  & kspt,kstep,kinc,iloop,nmethod,user_hardening,user_creep
45 !
46  real*8 elconloc(21),elas(21),emec(6),beta(6),stre(6),
47  & vj,plconloc(802),stbl(6),stril(6),xitril(6),
48  & ee,un,um,al,xk,cop,umb,umbb,dxitril,f0,d0,f1,d1,d2,xg(3,3),
49  & xs(3,3),xx(3,3),xn(3,3),xd(3,3),cpl(6),c(6),ci(6),
50  & c1,c2,c3,c4,c5,c6,c7,c8,c9,cplb(6),stblb(6),
51  & ftrial,xiso(200),yiso(200),xkin(200),ykin(200),
52  & fiso,dfiso,fkin,dfkin,fiso0,fkin0,ep,t1l,dtime,
53  & epini,a1,dsvm,xxa,xxn,vj2,vj23,sc(3,3),depvisc,
54  & cop1,cop2,fu1,fu2,fu,dcop,time,ttime,eloc(6),
55  & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
56  & g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11,g12,g13,g14,g15,g16,
57  & g17,g18,g28,g29,g30,g31,g32,g33,decra(5),deswa(5),serd,
58  & esw(2),ec(2),p,qtild,predef(1),dpred(1),timeabq(2),pgauss(3),
59  & dtemp,pnewdt
60 !
61 c data kk /1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3,
62 c & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,3,3,1,3,
63 c & 1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,1,2,2,3,1,3,2,3,
64 c & 2,3,2,3/
65 !
66 c pnewdt=-1.d0
67  leximp=1
68  lend=2
69  user_creep=0
70 !
71 ! localizing the plastic fields
72 !
73  do i=1,6
74  cpl(i)=-2.d0*xstateini(1+i,iint,iel)
75  stbl(i)=xstateini(7+i,iint,iel)
76  enddo
77  do i=1,3
78  cpl(i)=cpl(i)+1.d0
79  enddo
80  epini=xstateini(1,iint,iel)
81 !
82  ee=elconloc(1)
83  un=elconloc(2)
84  um=ee/(1.d0+un)
85  al=um*un/(1.d0-2.d0*un)
86  xk=al+um/3.d0
87  um=um/2.d0
88 !
89  ep=epini
90 !
91 ! right Cauchy-Green tensor (eloc contains the Lagrange strain,
92 ! including thermal strain)
93 !
94  c(1)=2.d0*emec(1)+1.d0
95  c(2)=2.d0*emec(2)+1.d0
96  c(3)=2.d0*emec(3)+1.d0
97  c(4)=2.d0*emec(4)
98  c(5)=2.d0*emec(5)
99  c(6)=2.d0*emec(6)
100 !
101 ! calculating the Jacobian
102 !
103  vj=c(1)*(c(2)*c(3)-c(6)*c(6))
104  & -c(4)*(c(4)*c(3)-c(6)*c(5))
105  & +c(5)*(c(4)*c(6)-c(2)*c(5))
106  if(vj.gt.1.d-30) then
107  vj=dsqrt(vj)
108  else
109 ccc write(*,*) '*WARNING in incplas: deformation inside-out'
110 !
111 ! deformation is reset to zero in order to continue the
112 ! calculation. Alternatively, a flag could be set forcing
113 ! a reiteration of the increment with a smaller size (to
114 ! be done)
115 !
116  c(1)=1.d0
117  c(2)=1.d0
118  c(3)=1.d0
119  c(4)=0.d0
120  c(5)=0.d0
121  c(6)=0.d0
122  vj=1.d0
123  endif
124 !
125 ! check whether the user activated a viscous calculation
126 ! (*VISCO instead of *STATIC)
127 !
128  if((nmethod.ne.1).or.(ithermal.eq.3)) then
129  ivisco=1
130  else
131  ivisco=0
132  endif
133 !
134 ! check for user subroutines
135 !
136  if((plconloc(801).lt.0.8d0).and.(plconloc(802).lt.0.8d0)) then
137  user_hardening=1
138  else
139  user_hardening=0
140  endif
141 c if(kode.eq.-52) then
142  if((kode.eq.-52).and.(ivisco.eq.1)) then
143  if(elconloc(3).lt.0.d0) then
144  user_creep=1
145  else
146  user_creep=0
147  xxa=elconloc(3)*(ttime+time)**elconloc(5)
148  if(xxa.lt.1.d-20) xxa=1.d-20
149  xxn=elconloc(4)
150  a1=xxa*dtime
151  endif
152  endif
153 !
154 ! inversion of the right Cauchy-Green tensor
155 !
156  vj2=vj*vj
157  ci(1)=(c(2)*c(3)-c(6)*c(6))/vj2
158  ci(2)=(c(1)*c(3)-c(5)*c(5))/vj2
159  ci(3)=(c(1)*c(2)-c(4)*c(4))/vj2
160  ci(4)=(c(5)*c(6)-c(4)*c(3))/vj2
161  ci(5)=(c(4)*c(6)-c(2)*c(5))/vj2
162  ci(6)=(c(4)*c(5)-c(1)*c(6))/vj2
163 !
164 ! reducing the plastic right Cauchy-Green tensor and
165 ! the back stress to "isochoric" quantities (b stands
166 ! for bar)
167 !
168  vj23=vj**(2.d0/3.d0)
169  do i=1,6
170  cplb(i)=cpl(i)/vj23
171  stblb(i)=stbl(i)/vj23
172  enddo
173 !
174 ! calculating the (n+1) trace and the (n+1) deviation of
175 ! the (n) "isochoric" plastic right Cauchy-Green tensor
176 !
177  umb=(c(1)*cplb(1)+c(2)*cplb(2)+c(3)*cplb(3)+
178  & 2.d0*(c(4)*cplb(4)+c(5)*cplb(5)+c(6)*cplb(6)))/3.d0
179  do i=1,6
180  cplb(i)=cplb(i)-umb*ci(i)
181  enddo
182 !
183 ! calculating the (n+1) trace and the (n+1) deviation of
184 ! the (n) "isochoric" back stress tensor
185 !
186  umbb=(c(1)*stblb(1)+c(2)*stblb(2)+c(3)*stblb(3)+
187  & 2.d0*(c(4)*stblb(4)+c(5)*stblb(5)+c(6)*stblb(6)))/3.d0
188  do i=1,6
189  stblb(i)=stblb(i)-umbb*ci(i)
190  enddo
191 !
192 ! calculating the trial stress
193 !
194  do i=1,6
195  stril(i)=um*cplb(i)-beta(i)
196  enddo
197 !
198 ! calculating the trial radius vector of the yield surface
199 !
200  do i=1,6
201  xitril(i)=stril(i)-stblb(i)
202  enddo
203 !
204  sc(1,1)=xitril(1)*c(1)+xitril(4)*c(4)+xitril(5)*c(5)
205  sc(1,2)=xitril(1)*c(4)+xitril(4)*c(2)+xitril(5)*c(6)
206  sc(1,3)=xitril(1)*c(5)+xitril(4)*c(6)+xitril(5)*c(3)
207  sc(2,1)=xitril(4)*c(1)+xitril(2)*c(4)+xitril(6)*c(5)
208  sc(2,2)=xitril(4)*c(4)+xitril(2)*c(2)+xitril(6)*c(6)
209  sc(2,3)=xitril(4)*c(5)+xitril(2)*c(6)+xitril(6)*c(3)
210  sc(3,1)=xitril(5)*c(1)+xitril(6)*c(4)+xitril(3)*c(5)
211  sc(3,2)=xitril(5)*c(4)+xitril(6)*c(2)+xitril(3)*c(6)
212  sc(3,3)=xitril(5)*c(5)+xitril(6)*c(6)+xitril(3)*c(3)
213 !
214  dxitril=sc(1,1)*sc(1,1)+sc(2,2)*sc(2,2)+sc(3,3)*sc(3,3)+
215  & 2.d0*(sc(1,2)*sc(2,1)+sc(1,3)*sc(3,1)+sc(2,3)*sc(3,2))
216 c g1=c(6)
217 c g2=xitril(6)
218 c g3=xitril(3)
219 c g4=xitril(2)
220 c g5=c(5)
221 c g6=xitril(5)
222 c g7=xitril(4)
223 c g8=c(4)
224 c g9=c(3)
225 c g10=c(2)
226 c g11=c(1)
227 c g12=xitril(1)
228 c g13=g12*g11
229 c g14=g10*g4
230 c g15=g9*g3
231 c g16=g8*g7
232 c g17=g6*g5
233 c g18=g2*g1
234 c g28=4*(g16 + g15)
235 c g29=4*g13
236 c g30=4*g14
237 c g31=4*g6*g1
238 c g32=4*g8*g2
239 c g33=4*g7*g5
240 c dxitril=(g13*g13 + g14*g14 + g15*g15 + g16*(g30 + g29 + 2*
241 c & g16) + g17*(g29 + g28 + 2*g17) + g18*(g30 + g28 + 2*
242 c & g18 + 4*g17) + g11*g7*(g31 + 2*g10*g7) + g9*g6*(g32 +
243 c & 2*g11*g6) + g10*g2*(g33 + 2*g9*g2) + g8*g4*(g31 + 2*
244 c & g12*g8) + g12*g5*(g32 + 2*g5*g3) + g3*g1*(g33 + 2*g4*
245 c & g1))
246 !
247 ! in principal dxitril is a norm and cannot be less than zero
248 ! however, due to round-off it can happen anyway.
249 ! In order to avoid this, one would have to use the deformation
250 ! gradient and calculate tau_ij*tau_kl*g_ik*g_jl instead of
251 ! S_IJ*S_KL*C_IK*C_JL, however, this would increase the
252 ! calculational effort substantially.
253 !
254  if(dxitril.lt.0.d0) then
255 ccc write(*,*) '*WARNING in incplas: dxitril < 0'
256  dxitril=0.d0
257  else
258  dxitril=dsqrt(dxitril)
259  endif
260 !
261 ! restoring the hardening curves for the actual temperature
262 ! plconloc contains the true stresses. By multiplying by
263 ! the Jacobian, yiso and ykin are Kirchhoff stresses, as
264 ! required by the hyperelastic theory (cf. Simo, 1988).
265 !
266  niso=int(plconloc(801))
267  nkin=int(plconloc(802))
268  if(niso.ne.0) then
269  do i=1,niso
270  xiso(i)=plconloc(2*i-1)
271  yiso(i)=vj*plconloc(2*i)
272  enddo
273  endif
274  if(nkin.ne.0) then
275  do i=1,nkin
276  xkin(i)=plconloc(399+2*i)
277  ykin(i)=vj*plconloc(400+2*i)
278  enddo
279  endif
280 !
281 ! if no viscous calculation is performed a pure creep calculatino
282 ! (without plasticity) is reduced to an elastic calculation
283 !
284  ielastic=0
285  if(ivisco.eq.0) then
286  if(niso.eq.2) then
287  if((dabs(yiso(1)).lt.1.d-10).and.
288  & (dabs(yiso(2)).lt.1.d-10)) then
289  ielastic=1
290  endif
291  endif
292  endif
293 !
294 ! check for yielding
295 !
296  if(user_hardening.eq.1) then
297  call uhardening(amat,iel,iint,t1l,epini,ep,dtime,
298  & fiso,dfiso,fkin,dfkin)
299  fiso=fiso*vj
300  else
301  if(niso.ne.0) then
302  call ident(xiso,ep,niso,id)
303  if(id.eq.0) then
304  fiso=yiso(1)
305  elseif(id.eq.niso) then
306  fiso=yiso(niso)
307  else
308  dfiso=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
309  fiso=yiso(id)+dfiso*(ep-xiso(id))
310  endif
311  elseif(nkin.ne.0) then
312  fiso=ykin(1)
313  else
314  fiso=0.d0
315  endif
316  endif
317 !
318  ftrial=dxitril-dsqrt(2.d0/3.d0)*fiso
319  if((ftrial.le.1.d-10).or.(ielas.eq.1).or.(ielastic.eq.1)
320  & .or.(dtime.lt.1.d-30)) then
321 !
322 ! no plastic deformation
323 ! beta contains the Cauchy residual stresses
324 !
325 c write(*,*) 'no plastic deformation'
326  c8=xk*(vj2-1.d0)/2.d0
327 !
328 ! residual stresses are de facto PK2 stresses
329 ! (Piola-Kirchhoff of the second kind)
330 !
331  stre(1)=c8*ci(1)+stril(1)-beta(1)
332  stre(2)=c8*ci(2)+stril(2)-beta(2)
333  stre(3)=c8*ci(3)+stril(3)-beta(3)
334  stre(4)=c8*ci(4)+stril(4)-beta(4)
335  stre(5)=c8*ci(5)+stril(5)-beta(5)
336  stre(6)=c8*ci(6)+stril(6)-beta(6)
337 !
338 ! updating the plastic fields
339 !
340  do i=1,3
341  cpl(i)=cpl(i)-1.d0
342  enddo
343  do i=1,6
344  xstate(1+i,iint,iel)=-cpl(i)/2.d0
345  xstate(7+i,iint,iel)=stbl(i)
346  enddo
347  xstate(1,iint,iel)=ep
348 !
349  if(icmd.ne.3) then
350 !
351  umb=um*umb
352 !
353 ! calculating the local stiffness matrix
354 !
355  xg(1,1)=(c(2)*c(3)-c(6)*c(6))/vj2
356  xg(2,2)=(c(1)*c(3)-c(5)*c(5))/vj2
357  xg(3,3)=(c(1)*c(2)-c(4)*c(4))/vj2
358  xg(1,2)=(c(5)*c(6)-c(4)*c(3))/vj2
359  xg(1,3)=(c(4)*c(6)-c(2)*c(5))/vj2
360  xg(2,3)=(c(4)*c(5)-c(1)*c(6))/vj2
361  xg(2,1)=xg(1,2)
362  xg(3,1)=xg(1,3)
363  xg(3,2)=xg(2,3)
364 !
365  xs(1,1)=stril(1)
366  xs(2,2)=stril(2)
367  xs(3,3)=stril(3)
368  xs(1,2)=stril(4)
369  xs(2,1)=stril(4)
370  xs(1,3)=stril(5)
371  xs(3,1)=stril(5)
372  xs(2,3)=stril(6)
373  xs(3,2)=stril(6)
374 !
375  elas(1)=umb*(xg(1,1)*xg(1,1)+xg(1,1)*xg(1,1)-
376  & 2.d0*xg(1,1)*xg(1,1)/3.d0)
377  & -2.d0*(xs(1,1)*xg(1,1)+xg(1,1)*xs(1,1))/3.d0
378  & +xk*vj2*xg(1,1)*xg(1,1)
379  & -xk*(vj2-1.d0)*(xg(1,1)*xg(1,1)
380  & +xg(1,1)*xg(1,1))/2.d0
381  elas(2)=umb*(xg(1,2)*xg(1,2)+xg(1,2)*xg(1,2)-
382  & 2.d0*xg(1,1)*xg(2,2)/3.d0)
383  & -2.d0*(xs(1,1)*xg(2,2)+xg(1,1)*xs(2,2))/3.d0
384  & +xk*vj2*xg(1,1)*xg(2,2)
385  & -xk*(vj2-1.d0)*(xg(1,2)*xg(1,2)
386  & +xg(1,2)*xg(1,2))/2.d0
387  elas(3)=umb*(xg(2,2)*xg(2,2)+xg(2,2)*xg(2,2)-
388  & 2.d0*xg(2,2)*xg(2,2)/3.d0)
389  & -2.d0*(xs(2,2)*xg(2,2)+xg(2,2)*xs(2,2))/3.d0
390  & +xk*vj2*xg(2,2)*xg(2,2)
391  & -xk*(vj2-1.d0)*(xg(2,2)*xg(2,2)
392  & +xg(2,2)*xg(2,2))/2.d0
393  elas(4)=umb*(xg(1,3)*xg(1,3)+xg(1,3)*xg(1,3)-
394  & 2.d0*xg(1,1)*xg(3,3)/3.d0)
395  & -2.d0*(xs(1,1)*xg(3,3)+xg(1,1)*xs(3,3))/3.d0
396  & +xk*vj2*xg(1,1)*xg(3,3)
397  & -xk*(vj2-1.d0)*(xg(1,3)*xg(1,3)
398  & +xg(1,3)*xg(1,3))/2.d0
399  elas(5)=umb*(xg(2,3)*xg(2,3)+xg(2,3)*xg(2,3)-
400  & 2.d0*xg(2,2)*xg(3,3)/3.d0)
401  & -2.d0*(xs(2,2)*xg(3,3)+xg(2,2)*xs(3,3))/3.d0
402  & +xk*vj2*xg(2,2)*xg(3,3)
403  & -xk*(vj2-1.d0)*(xg(2,3)*xg(2,3)
404  & +xg(2,3)*xg(2,3))/2.d0
405  elas(6)=umb*(xg(3,3)*xg(3,3)+xg(3,3)*xg(3,3)-
406  & 2.d0*xg(3,3)*xg(3,3)/3.d0)
407  & -2.d0*(xs(3,3)*xg(3,3)+xg(3,3)*xs(3,3))/3.d0
408  & +xk*vj2*xg(3,3)*xg(3,3)
409  & -xk*(vj2-1.d0)*(xg(3,3)*xg(3,3)
410  & +xg(3,3)*xg(3,3))/2.d0
411  elas(7)=umb*(xg(1,1)*xg(1,2)+xg(1,2)*xg(1,1)-
412  & 2.d0*xg(1,1)*xg(1,2)/3.d0)
413  & -2.d0*(xs(1,1)*xg(1,2)+xg(1,1)*xs(1,2))/3.d0
414  & +xk*vj2*xg(1,1)*xg(1,2)
415  & -xk*(vj2-1.d0)*(xg(1,1)*xg(1,2)
416  & +xg(1,2)*xg(1,1))/2.d0
417  elas(8)=umb*(xg(2,1)*xg(2,2)+xg(2,2)*xg(2,1)-
418  & 2.d0*xg(2,2)*xg(1,2)/3.d0)
419  & -2.d0*(xs(2,2)*xg(1,2)+xg(2,2)*xs(1,2))/3.d0
420  & +xk*vj2*xg(2,2)*xg(1,2)
421  & -xk*(vj2-1.d0)*(xg(2,1)*xg(2,2)
422  & +xg(2,2)*xg(2,1))/2.d0
423  elas(9)=umb*(xg(3,1)*xg(3,2)+xg(3,2)*xg(3,1)-
424  & 2.d0*xg(3,3)*xg(1,2)/3.d0)
425  & -2.d0*(xs(3,3)*xg(1,2)+xg(3,3)*xs(1,2))/3.d0
426  & +xk*vj2*xg(3,3)*xg(1,2)
427  & -xk*(vj2-1.d0)*(xg(3,1)*xg(3,2)
428  & +xg(3,2)*xg(3,1))/2.d0
429  elas(10)=umb*(xg(1,1)*xg(2,2)+xg(1,2)*xg(2,1)-
430  & 2.d0*xg(1,2)*xg(1,2)/3.d0)
431  & -2.d0*(xs(1,2)*xg(1,2)+xg(1,2)*xs(1,2))/3.d0
432  & +xk*vj2*xg(1,2)*xg(1,2)
433  & -xk*(vj2-1.d0)*(xg(1,1)*xg(2,2)
434  & +xg(1,2)*xg(2,1))/2.d0
435  elas(11)=umb*(xg(1,1)*xg(1,3)+xg(1,3)*xg(1,1)-
436  & 2.d0*xg(1,1)*xg(1,3)/3.d0)
437  & -2.d0*(xs(1,1)*xg(1,3)+xg(1,1)*xs(1,3))/3.d0
438  & +xk*vj2*xg(1,1)*xg(1,3)
439  & -xk*(vj2-1.d0)*(xg(1,1)*xg(1,3)
440  & +xg(1,3)*xg(1,1))/2.d0
441  elas(12)=umb*(xg(2,1)*xg(2,3)+xg(2,3)*xg(2,1)-
442  & 2.d0*xg(2,2)*xg(1,3)/3.d0)
443  & -2.d0*(xs(2,2)*xg(1,3)+xg(2,2)*xs(1,3))/3.d0
444  & +xk*vj2*xg(2,2)*xg(1,3)
445  & -xk*(vj2-1.d0)*(xg(2,1)*xg(2,3)
446  & +xg(2,3)*xg(2,1))/2.d0
447  elas(13)=umb*(xg(3,1)*xg(3,3)+xg(3,3)*xg(3,1)-
448  & 2.d0*xg(3,3)*xg(1,3)/3.d0)
449  & -2.d0*(xs(3,3)*xg(1,3)+xg(3,3)*xs(1,3))/3.d0
450  & +xk*vj2*xg(3,3)*xg(1,3)
451  & -xk*(vj2-1.d0)*(xg(3,1)*xg(3,3)
452  & +xg(3,3)*xg(3,1))/2.d0
453  elas(14)=umb*(xg(1,1)*xg(2,3)+xg(1,3)*xg(2,1)-
454  & 2.d0*xg(1,2)*xg(1,3)/3.d0)
455  & -2.d0*(xs(1,2)*xg(1,3)+xg(1,2)*xs(1,3))/3.d0
456  & +xk*vj2*xg(1,2)*xg(1,3)
457  & -xk*(vj2-1.d0)*(xg(1,1)*xg(2,3)
458  & +xg(1,3)*xg(2,1))/2.d0
459  elas(15)=umb*(xg(1,1)*xg(3,3)+xg(1,3)*xg(3,1)-
460  & 2.d0*xg(1,3)*xg(1,3)/3.d0)
461  & -2.d0*(xs(1,3)*xg(1,3)+xg(1,3)*xs(1,3))/3.d0
462  & +xk*vj2*xg(1,3)*xg(1,3)
463  & -xk*(vj2-1.d0)*(xg(1,1)*xg(3,3)
464  & +xg(1,3)*xg(3,1))/2.d0
465  elas(16)=umb*(xg(1,2)*xg(1,3)+xg(1,3)*xg(1,2)-
466  & 2.d0*xg(1,1)*xg(2,3)/3.d0)
467  & -2.d0*(xs(1,1)*xg(2,3)+xg(1,1)*xs(2,3))/3.d0
468  & +xk*vj2*xg(1,1)*xg(2,3)
469  & -xk*(vj2-1.d0)*(xg(1,2)*xg(1,3)
470  & +xg(1,3)*xg(1,2))/2.d0
471  elas(17)=umb*(xg(2,2)*xg(2,3)+xg(2,3)*xg(2,2)-
472  & 2.d0*xg(2,2)*xg(2,3)/3.d0)
473  & -2.d0*(xs(2,2)*xg(2,3)+xg(2,2)*xs(2,3))/3.d0
474  & +xk*vj2*xg(2,2)*xg(2,3)
475  & -xk*(vj2-1.d0)*(xg(2,2)*xg(2,3)
476  & +xg(2,3)*xg(2,2))/2.d0
477  elas(18)=umb*(xg(3,2)*xg(3,3)+xg(3,3)*xg(3,2)-
478  & 2.d0*xg(3,3)*xg(2,3)/3.d0)
479  & -2.d0*(xs(3,3)*xg(2,3)+xg(3,3)*xs(2,3))/3.d0
480  & +xk*vj2*xg(3,3)*xg(2,3)
481  & -xk*(vj2-1.d0)*(xg(3,2)*xg(3,3)
482  & +xg(3,3)*xg(3,2))/2.d0
483  elas(19)=umb*(xg(1,2)*xg(2,3)+xg(1,3)*xg(2,2)-
484  & 2.d0*xg(1,2)*xg(2,3)/3.d0)
485  & -2.d0*(xs(1,2)*xg(2,3)+xg(1,2)*xs(2,3))/3.d0
486  & +xk*vj2*xg(1,2)*xg(2,3)
487  & -xk*(vj2-1.d0)*(xg(1,2)*xg(2,3)
488  & +xg(1,3)*xg(2,2))/2.d0
489  elas(20)=umb*(xg(1,2)*xg(3,3)+xg(1,3)*xg(3,2)-
490  & 2.d0*xg(1,3)*xg(2,3)/3.d0)
491  & -2.d0*(xs(1,3)*xg(2,3)+xg(1,3)*xs(2,3))/3.d0
492  & +xk*vj2*xg(1,3)*xg(2,3)
493  & -xk*(vj2-1.d0)*(xg(1,2)*xg(3,3)
494  & +xg(1,3)*xg(3,2))/2.d0
495  elas(21)=umb*(xg(2,2)*xg(3,3)+xg(2,3)*xg(3,2)-
496  & 2.d0*xg(2,3)*xg(2,3)/3.d0)
497  & -2.d0*(xs(2,3)*xg(2,3)+xg(2,3)*xs(2,3))/3.d0
498  & +xk*vj2*xg(2,3)*xg(2,3)
499  & -xk*(vj2-1.d0)*(xg(2,2)*xg(3,3)
500  & +xg(2,3)*xg(3,2))/2.d0
501 !
502  endif
503 !
504  return
505  endif
506 !
507 ! plastic deformation
508 !
509  umb=um*umb
510  umbb=umb-umbb
511 !
512 ! calculating the consistency parameter
513 !
514  c1=2.d0/3.d0
515  c2=dsqrt(c1)
516  c3=c1/um
517  c4=c2/um
518 !
519  iloop=0
520  cop=0.d0
521 !
522  loop: do
523  iloop=iloop+1
524  ep=epini+c2*cop
525 !
526  if(user_hardening.eq.1) then
527  call uhardening(amat,iel,iint,t1l,epini,ep,dtime,
528  & fiso,dfiso,fkin,dfkin)
529  fiso=fiso*vj
530  dfiso=dfiso*vj
531  fkin=fkin*vj
532  dfkin=dfkin*vj
533  else
534  if(niso.ne.0) then
535  call ident(xiso,ep,niso,id)
536  if(id.eq.0) then
537  fiso=yiso(1)
538  dfiso=0.d0
539  elseif(id.eq.niso) then
540  fiso=yiso(niso)
541  dfiso=0.d0
542  else
543  dfiso=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
544  fiso=yiso(id)+dfiso*(ep-xiso(id))
545  endif
546  elseif(nkin.ne.0) then
547  fiso=ykin(1)
548  dfiso=0.d0
549  else
550  fiso=0.d0
551  dfiso=0.d0
552  endif
553 !
554  if(nkin.ne.0) then
555  call ident(xkin,ep,nkin,id)
556  if(id.eq.0) then
557  fkin=ykin(1)
558  dfkin=0.d0
559  elseif(id.eq.nkin) then
560  fkin=ykin(nkin)
561  dfkin=0.d0
562  else
563  dfkin=(ykin(id+1)-ykin(id))/(xkin(id+1)-xkin(id))
564  fkin=ykin(id)+dfkin*(ep-xkin(id))
565  endif
566  elseif(niso.ne.0) then
567  fkin=yiso(1)
568  dfkin=0.d0
569  else
570  fkin=0.d0
571  dfkin=0.d0
572  endif
573  endif
574 !
575  if(dabs(cop).lt.1.d-10) then
576  fiso0=fiso
577  fkin0=fkin
578  endif
579 !
580  if((kode.eq.-51).or.(ivisco.eq.0)) then
581  dcop=(ftrial-c2*(fiso-fiso0)
582  & -umbb*(2.d0*cop+c4*(fkin-fkin0)))/
583  & (-c1*dfiso-umbb*(2.d0+c3*dfkin))
584  else
585  if(user_creep.eq.1) then
586  if(ithermal.eq.0) then
587 ccc write(*,*) '*ERROR in incplas: no temperature defined'
588  call exit(201)
589  endif
590  timeabq(1)=time
591  timeabq(2)=ttime+time
592  qtild=(ftrial-c2*(fiso-fiso0)
593  & -umbb*(2.d0*cop+c4*(fkin-fkin0)))/(c2*vj)
594 !
595 ! the Von Mises stress must be positive
596 !
597  if(qtild.lt.1.d-10) qtild=1.d-10
598  ec(1)=epini
599  call creep(decra,deswa,xstateini(1,iint,iel),serd,ec,
600  & esw,p,qtild,t1l,dtemp,predef,dpred,timeabq,dtime,
601  & amat,leximp,lend,pgauss,nstate_,iel,iint,layer,kspt,
602  & kstep,kinc)
603  dsvm=1.d0/decra(5)
604  dcop=-(decra(1)-c2*cop)/
605  & (c2*(decra(5)*(dfiso+umbb*(3.d0+dfkin/um))+1.d0))
606  else
607  qtild=(ftrial-c2*(fiso-fiso0)
608  & -umbb*(2.d0*cop+c4*(fkin-fkin0)))/(c2*vj)
609 !
610 ! the Von Mises stress must be positive
611 !
612  if(qtild.lt.1.d-10) qtild=1.d-10
613  decra(1)=a1*qtild**xxn
614  decra(5)=xxn*decra(1)/qtild
615  dsvm=1.d0/decra(5)
616  dcop=-(decra(1)-c2*cop)/
617  & (c2*(decra(5)*(dfiso+umbb*(3.d0+dfkin/um))+1.d0))
618  endif
619  endif
620  cop=cop-dcop
621 !
622  if((dabs(dcop).lt.cop*1.d-4).or.
623  & (dabs(dcop).lt.1.d-10)) exit
624 !
625 ! check for endless loops or a negative consistency
626 ! parameter
627 !
628  if((iloop.gt.15).or.(cop.le.0.d0)) then
629  iloop=1
630  cop=0.d0
631  do
632  if(iloop.gt.100) then
633 c NOTE: write statements cause problems for
634 c parallellized execution
635 c write(*,*) '*WARNING in incplas: material loop'
636 c write(*,*) ' did not converge in integration'
637 c write(*,*) ' point',iint,'in element',iel,';'
638 c write(*,*) ' the increment size is reduced'
639 c write(*,*)
640  pnewdt=0.25d0
641  return
642  endif
643  ep=epini+c2*cop
644 !
645  if(user_hardening.eq.1) then
646  call uhardening(amat,iel,iint,t1l,epini,ep,dtime,
647  & fiso,dfiso,fkin,dfkin)
648  fiso=fiso*vj
649  fkin=fkin*vj
650  else
651  if(niso.ne.0) then
652  call ident(xiso,ep,niso,id)
653  if(id.eq.0) then
654  fiso=yiso(1)
655  elseif(id.eq.niso) then
656  fiso=yiso(niso)
657  else
658  dfiso=(yiso(id+1)-yiso(id))/
659  & (xiso(id+1)-xiso(id))
660  fiso=yiso(id)+dfiso*(ep-xiso(id))
661  endif
662  elseif(nkin.ne.0) then
663  fiso=ykin(1)
664  else
665  fiso=0.d0
666  endif
667 !
668  if(nkin.ne.0) then
669  call ident(xkin,ep,nkin,id)
670  if(id.eq.0) then
671  fkin=ykin(1)
672  elseif(id.eq.nkin) then
673  fkin=ykin(nkin)
674  else
675  dfkin=(ykin(id+1)-ykin(id))/
676  & (xkin(id+1)-xkin(id))
677  fkin=ykin(id)+dfkin*(ep-xkin(id))
678  endif
679  elseif(niso.ne.0) then
680  fkin=yiso(1)
681  else
682  fkin=0.d0
683  endif
684  endif
685 !
686  if(dabs(cop).lt.1.d-10) then
687  fiso0=fiso
688  fkin0=fkin
689  endif
690 !
691  if((kode.eq.-51).or.(ivisco.eq.0)) then
692  fu=(ftrial-c2*(fiso-fiso0)
693  & -umbb*(2.d0*cop+c4*(fkin-fkin0)))
694  else
695  if(user_creep.eq.1) then
696  timeabq(1)=time
697  timeabq(2)=ttime+time
698  qtild=(ftrial-c2*(fiso-fiso0)
699  & -umbb*(2.d0*cop+c4*(fkin-fkin0)))/(c2*vj)
700 !
701 ! the Von Mises stress must be positive
702 !
703  if(qtild.lt.1.d-10) qtild=1.d-10
704  ec(1)=epini
705  call creep(decra,deswa,xstateini(1,iint,iel),serd,
706  & ec,esw,p,qtild,t1l,dtemp,predef,dpred,timeabq,
707  & dtime,amat,leximp,lend,pgauss,nstate_,iel,
708  & iint,layer,kspt,kstep,kinc)
709  dsvm=1.d0/decra(5)
710  fu=decra(1)-c2*cop
711  else
712  qtild=(ftrial-c2*(fiso-fiso0)
713  & -umbb*(2.d0*cop+c4*(fkin-fkin0)))/(c2*vj)
714 !
715 ! the Von Mises stress must be positive
716 !
717  if(qtild.lt.1.d-10) qtild=1.d-10
718  decra(1)=a1*qtild**xxn
719  decra(5)=xxn*decra(1)/qtild
720  dsvm=1.d0/decra(5)
721  fu=decra(1)-c2*cop
722  endif
723  endif
724 !
725  if(iloop.eq.1) then
726 c write(*,*) 'cop,fu ',cop,fu
727  cop1=0.d0
728  fu1=fu
729  iloop=2
730  cop=1.d-10
731  elseif(iloop.eq.2) then
732  if(fu*fu1.le.0.d0) then
733 c write(*,*) cop,fu
734  iloop=3
735  fu2=fu
736  cop2=cop
737  cop=(cop1+cop2)/2.d0
738  dcop=(cop2-cop1)/2.d0
739  else
740 c write(*,*) cop,fu
741  cop=cop*10.d0
742  if(cop.gt.100.d0) then
743  write(*,*) '*ERROR: no convergence in incplas'
744  pnewdt=0.25d0
745  return
746 c call exit(201)
747  endif
748  endif
749  else
750 c write(*,*) cop,fu
751  if(fu*fu1.ge.0.d0) then
752  cop1=cop
753  fu1=fu
754  else
755  cop2=cop
756  fu2=fu
757  endif
758  cop=(cop1+cop2)/2.d0
759  dcop=(cop2-cop1)/2.d0
760  if((dabs(dcop).lt.cop*1.d-4).or.
761  & (dabs(dcop).lt.1.d-10)) exit loop
762  endif
763  enddo
764  endif
765 !
766  enddo loop
767 !
768 ! updating the equivalent plastic strain
769 !
770  ep=epini+c2*cop
771 !
772 ! updating the back stress
773 !
774  c5=2.d0*umbb*cop/dxitril
775  c6=c5/(3.d0*um)
776  c7=c6*dfkin*vj23
777  do i=1,6
778  stbl(i)=stbl(i)+c7*xitril(i)
779  enddo
780 !
781 ! updating the stress
782 ! vj: Jacobian of the total deformation gradient
783 !
784  c8=xk*(vj2-1.d0)/2.d0
785 !
786  do i=1,6
787  stre(i)=c8*ci(i)-beta(i)+stril(i)-c5*xitril(i)
788  enddo
789 !
790 ! updating the plastic right Cauchy-Green tensor
791 !
792  c9=c6*3.d0*vj23
793  do i=1,6
794  cpl(i)=cpl(i)-c9*xitril(i)
795  enddo
796 !
797  if(icmd.ne.3) then
798 !
799 ! calculating the local stiffness matrix
800 !
801  xg(1,1)=(c(2)*c(3)-c(6)*c(6))/vj2
802  xg(2,2)=(c(1)*c(3)-c(5)*c(5))/vj2
803  xg(3,3)=(c(1)*c(2)-c(4)*c(4))/vj2
804  xg(1,2)=(c(5)*c(6)-c(4)*c(3))/vj2
805  xg(1,3)=(c(4)*c(6)-c(2)*c(5))/vj2
806  xg(2,3)=(c(4)*c(5)-c(1)*c(6))/vj2
807  xg(2,1)=xg(1,2)
808  xg(3,1)=xg(1,3)
809  xg(3,2)=xg(2,3)
810 !
811  xs(1,1)=stril(1)
812  xs(2,2)=stril(2)
813  xs(3,3)=stril(3)
814  xs(1,2)=stril(4)
815  xs(2,1)=stril(4)
816  xs(1,3)=stril(5)
817  xs(3,1)=stril(5)
818  xs(2,3)=stril(6)
819  xs(3,2)=stril(6)
820 !
821  f0=2.d0*umbb*cop/dxitril
822  d0=1.d0+(dfkin/um+dfiso/umbb)/3.d0
823 !
824 ! creep contribution
825 !
826  if((kode.eq.-52).and.(ivisco.eq.1)) then
827  d0=d0+dsvm/(3.d0*umbb)
828  endif
829 !
830  f1=1.d0/d0-f0
831  d1=2.d0*f1*umbb-((1.d0+dfkin/(3.d0*um))/d0-1.d0)*
832  & 4.d0*cop*dxitril/3.d0
833  d2=2d0*dxitril*f1
834 !
835  xx(1,1)=xitril(1)
836  xx(2,2)=xitril(2)
837  xx(3,3)=xitril(3)
838  xx(1,2)=xitril(4)
839  xx(2,1)=xitril(4)
840  xx(1,3)=xitril(5)
841  xx(3,1)=xitril(5)
842  xx(2,3)=xitril(6)
843  xx(3,2)=xitril(6)
844 !
845  xn(1,1)=xitril(1)/dxitril
846  xn(2,2)=xitril(2)/dxitril
847  xn(3,3)=xitril(3)/dxitril
848  xn(1,2)=xitril(4)/dxitril
849  xn(2,1)=xitril(4)/dxitril
850  xn(1,3)=xitril(5)/dxitril
851  xn(3,1)=xitril(5)/dxitril
852  xn(2,3)=xitril(6)/dxitril
853  xn(3,2)=xitril(6)/dxitril
854 !
855  do i=1,3
856  do j=i,3
857  xd(i,j)=xn(i,1)*xn(1,j)*c(1)+xn(i,1)*xn(2,j)*c(4)+
858  & xn(i,1)*xn(3,j)*c(5)+xn(i,2)*xn(1,j)*c(4)+
859  & xn(i,2)*xn(2,j)*c(2)+xn(i,2)*xn(3,j)*c(6)+
860  & xn(i,3)*xn(1,j)*c(5)+xn(i,3)*xn(2,j)*c(6)+
861  & xn(i,3)*xn(3,j)*c(3)
862  enddo
863  enddo
864  xd(2,1)=xd(1,2)
865  xd(3,1)=xd(1,3)
866  xd(3,2)=xd(2,3)
867 !
868 ! deviatoric part
869 !
870  c1=(xd(1,1)*c(1)+xd(2,2)*c(2)+xd(3,3)*c(3)+
871  & 2.d0*(xd(1,2)*c(4)+xd(1,3)*c(5)+xd(2,3)*c(6)))/3.d0
872  do i=1,3
873  do j=i,3
874  xd(i,j)=xd(i,j)-c1*xg(i,j)
875  enddo
876  enddo
877  xd(2,1)=xd(1,2)
878  xd(3,1)=xd(1,3)
879  xd(3,2)=xd(2,3)
880 !
881  elas(1)=(umb-f0*umbb)*(xg(1,1)*xg(1,1)+xg(1,1)*xg(1,1)-
882  & 2.d0*xg(1,1)*xg(1,1)/3.d0)
883  & -2.d0*(xs(1,1)*xg(1,1)+xg(1,1)*xs(1,1))/3.d0
884  & +f0*2.d0*(xx(1,1)*xg(1,1)+xg(1,1)*xx(1,1))/3.d0
885  & -d1*xn(1,1)*xn(1,1)-d2*(xn(1,1)*xd(1,1)+
886  & xd(1,1)*xn(1,1))/2.d0+xk*vj2*xg(1,1)*xg(1,1)
887  & -xk*(vj2-1.d0)*(xg(1,1)*xg(1,1)+xg(1,1)*xg(1,1))/2.d0
888  elas(2)=(umb-f0*umbb)*(xg(1,2)*xg(1,2)+xg(1,2)*xg(1,2)-
889  & 2.d0*xg(1,1)*xg(2,2)/3.d0)
890  & -2.d0*(xs(1,1)*xg(2,2)+xg(1,1)*xs(2,2))/3.d0
891  & +f0*2.d0*(xx(1,1)*xg(2,2)+xg(1,1)*xx(2,2))/3.d0
892  & -d1*xn(1,1)*xn(2,2)-d2*(xn(1,1)*xd(2,2)+
893  & xd(1,1)*xn(2,2))/2.d0+xk*vj2*xg(1,1)*xg(2,2)
894  & -xk*(vj2-1.d0)*(xg(1,2)*xg(1,2)+xg(1,2)*xg(1,2))/2.d0
895  elas(3)=(umb-f0*umbb)*(xg(2,2)*xg(2,2)+xg(2,2)*xg(2,2)-
896  & 2.d0*xg(2,2)*xg(2,2)/3.d0)
897  & -2.d0*(xs(2,2)*xg(2,2)+xg(2,2)*xs(2,2))/3.d0
898  & +f0*2.d0*(xx(2,2)*xg(2,2)+xg(2,2)*xx(2,2))/3.d0
899  & -d1*xn(2,2)*xn(2,2)-d2*(xn(2,2)*xd(2,2)+
900  & xd(2,2)*xn(2,2))/2.d0+xk*vj2*xg(2,2)*xg(2,2)
901  & -xk*(vj2-1.d0)*(xg(2,2)*xg(2,2)+xg(2,2)*xg(2,2))/2.d0
902  elas(4)=(umb-f0*umbb)*(xg(1,3)*xg(1,3)+xg(1,3)*xg(1,3)-
903  & 2.d0*xg(1,1)*xg(3,3)/3.d0)
904  & -2.d0*(xs(1,1)*xg(3,3)+xg(1,1)*xs(3,3))/3.d0
905  & +f0*2.d0*(xx(1,1)*xg(3,3)+xg(1,1)*xx(3,3))/3.d0
906  & -d1*xn(1,1)*xn(3,3)-d2*(xn(1,1)*xd(3,3)+
907  & xd(1,1)*xn(3,3))/2.d0+xk*vj2*xg(1,1)*xg(3,3)
908  & -xk*(vj2-1.d0)*(xg(1,3)*xg(1,3)+xg(1,3)*xg(1,3))/2.d0
909  elas(5)=(umb-f0*umbb)*(xg(2,3)*xg(2,3)+xg(2,3)*xg(2,3)-
910  & 2.d0*xg(2,2)*xg(3,3)/3.d0)
911  & -2.d0*(xs(2,2)*xg(3,3)+xg(2,2)*xs(3,3))/3.d0
912  & +f0*2.d0*(xx(2,2)*xg(3,3)+xg(2,2)*xx(3,3))/3.d0
913  & -d1*xn(2,2)*xn(3,3)-d2*(xn(2,2)*xd(3,3)+
914  & xd(2,2)*xn(3,3))/2.d0+xk*vj2*xg(2,2)*xg(3,3)
915  & -xk*(vj2-1.d0)*(xg(2,3)*xg(2,3)+xg(2,3)*xg(2,3))/2.d0
916  elas(6)=(umb-f0*umbb)*(xg(3,3)*xg(3,3)+xg(3,3)*xg(3,3)-
917  & 2.d0*xg(3,3)*xg(3,3)/3.d0)
918  & -2.d0*(xs(3,3)*xg(3,3)+xg(3,3)*xs(3,3))/3.d0
919  & +f0*2.d0*(xx(3,3)*xg(3,3)+xg(3,3)*xx(3,3))/3.d0
920  & -d1*xn(3,3)*xn(3,3)-d2*(xn(3,3)*xd(3,3)+
921  & xd(3,3)*xn(3,3))/2.d0+xk*vj2*xg(3,3)*xg(3,3)
922  & -xk*(vj2-1.d0)*(xg(3,3)*xg(3,3)+xg(3,3)*xg(3,3))/2.d0
923  elas(7)=(umb-f0*umbb)*(xg(1,1)*xg(1,2)+xg(1,2)*xg(1,1)-
924  & 2.d0*xg(1,1)*xg(1,2)/3.d0)
925  & -2.d0*(xs(1,1)*xg(1,2)+xg(1,1)*xs(1,2))/3.d0
926  & +f0*2.d0*(xx(1,1)*xg(1,2)+xg(1,1)*xx(1,2))/3.d0
927  & -d1*xn(1,1)*xn(1,2)-d2*(xn(1,1)*xd(1,2)+
928  & xd(1,1)*xn(1,2))/2.d0+xk*vj2*xg(1,1)*xg(1,2)
929  & -xk*(vj2-1.d0)*(xg(1,1)*xg(1,2)+xg(1,2)*xg(1,1))/2.d0
930  elas(8)=(umb-f0*umbb)*(xg(2,1)*xg(2,2)+xg(2,2)*xg(2,1)-
931  & 2.d0*xg(2,2)*xg(1,2)/3.d0)
932  & -2.d0*(xs(2,2)*xg(1,2)+xg(2,2)*xs(1,2))/3.d0
933  & +f0*2.d0*(xx(2,2)*xg(1,2)+xg(2,2)*xx(1,2))/3.d0
934  & -d1*xn(2,2)*xn(1,2)-d2*(xn(2,2)*xd(1,2)+
935  & xd(2,2)*xn(1,2))/2.d0+xk*vj2*xg(2,2)*xg(1,2)
936  & -xk*(vj2-1.d0)*(xg(2,1)*xg(2,2)+xg(2,2)*xg(2,1))/2.d0
937  elas(9)=(umb-f0*umbb)*(xg(3,1)*xg(3,2)+xg(3,2)*xg(3,1)-
938  & 2.d0*xg(3,3)*xg(1,2)/3.d0)
939  & -2.d0*(xs(3,3)*xg(1,2)+xg(3,3)*xs(1,2))/3.d0
940  & +f0*2.d0*(xx(3,3)*xg(1,2)+xg(3,3)*xx(1,2))/3.d0
941  & -d1*xn(3,3)*xn(1,2)-d2*(xn(3,3)*xd(1,2)+
942  & xd(3,3)*xn(1,2))/2.d0+xk*vj2*xg(3,3)*xg(1,2)
943  & -xk*(vj2-1.d0)*(xg(3,1)*xg(3,2)+xg(3,2)*xg(3,1))/2.d0
944  elas(10)=(umb-f0*umbb)*(xg(1,1)*xg(2,2)+xg(1,2)*xg(2,1)-
945  & 2.d0*xg(1,2)*xg(1,2)/3.d0)
946  & -2.d0*(xs(1,2)*xg(1,2)+xg(1,2)*xs(1,2))/3.d0
947  & +f0*2.d0*(xx(1,2)*xg(1,2)+xg(1,2)*xx(1,2))/3.d0
948  & -d1*xn(1,2)*xn(1,2)-d2*(xn(1,2)*xd(1,2)+
949  & xd(1,2)*xn(1,2))/2.d0+xk*vj2*xg(1,2)*xg(1,2)
950  & -xk*(vj2-1.d0)*(xg(1,1)*xg(2,2)+xg(1,2)*xg(2,1))/2.d0
951  elas(11)=(umb-f0*umbb)*(xg(1,1)*xg(1,3)+xg(1,3)*xg(1,1)-
952  & 2.d0*xg(1,1)*xg(1,3)/3.d0)
953  & -2.d0*(xs(1,1)*xg(1,3)+xg(1,1)*xs(1,3))/3.d0
954  & +f0*2.d0*(xx(1,1)*xg(1,3)+xg(1,1)*xx(1,3))/3.d0
955  & -d1*xn(1,1)*xn(1,3)-d2*(xn(1,1)*xd(1,3)+
956  & xd(1,1)*xn(1,3))/2.d0+xk*vj2*xg(1,1)*xg(1,3)
957  & -xk*(vj2-1.d0)*(xg(1,1)*xg(1,3)+xg(1,3)*xg(1,1))/2.d0
958  elas(12)=(umb-f0*umbb)*(xg(2,1)*xg(2,3)+xg(2,3)*xg(2,1)-
959  & 2.d0*xg(2,2)*xg(1,3)/3.d0)
960  & -2.d0*(xs(2,2)*xg(1,3)+xg(2,2)*xs(1,3))/3.d0
961  & +f0*2.d0*(xx(2,2)*xg(1,3)+xg(2,2)*xx(1,3))/3.d0
962  & -d1*xn(2,2)*xn(1,3)-d2*(xn(2,2)*xd(1,3)+
963  & xd(2,2)*xn(1,3))/2.d0+xk*vj2*xg(2,2)*xg(1,3)
964  & -xk*(vj2-1.d0)*(xg(2,1)*xg(2,3)+xg(2,3)*xg(2,1))/2.d0
965  elas(13)=(umb-f0*umbb)*(xg(3,1)*xg(3,3)+xg(3,3)*xg(3,1)-
966  & 2.d0*xg(3,3)*xg(1,3)/3.d0)
967  & -2.d0*(xs(3,3)*xg(1,3)+xg(3,3)*xs(1,3))/3.d0
968  & +f0*2.d0*(xx(3,3)*xg(1,3)+xg(3,3)*xx(1,3))/3.d0
969  & -d1*xn(3,3)*xn(1,3)-d2*(xn(3,3)*xd(1,3)+
970  & xd(3,3)*xn(1,3))/2.d0+xk*vj2*xg(3,3)*xg(1,3)
971  & -xk*(vj2-1.d0)*(xg(3,1)*xg(3,3)+xg(3,3)*xg(3,1))/2.d0
972  elas(14)=(umb-f0*umbb)*(xg(1,1)*xg(2,3)+xg(1,3)*xg(2,1)-
973  & 2.d0*xg(1,2)*xg(1,3)/3.d0)
974  & -2.d0*(xs(1,2)*xg(1,3)+xg(1,2)*xs(1,3))/3.d0
975  & +f0*2.d0*(xx(1,2)*xg(1,3)+xg(1,2)*xx(1,3))/3.d0
976  & -d1*xn(1,2)*xn(1,3)-d2*(xn(1,2)*xd(1,3)+
977  & xd(1,2)*xn(1,3))/2.d0+xk*vj2*xg(1,2)*xg(1,3)
978  & -xk*(vj2-1.d0)*(xg(1,1)*xg(2,3)+xg(1,3)*xg(2,1))/2.d0
979  elas(15)=(umb-f0*umbb)*(xg(1,1)*xg(3,3)+xg(1,3)*xg(3,1)-
980  & 2.d0*xg(1,3)*xg(1,3)/3.d0)
981  & -2.d0*(xs(1,3)*xg(1,3)+xg(1,3)*xs(1,3))/3.d0
982  & +f0*2.d0*(xx(1,3)*xg(1,3)+xg(1,3)*xx(1,3))/3.d0
983  & -d1*xn(1,3)*xn(1,3)-d2*(xn(1,3)*xd(1,3)+
984  & xd(1,3)*xn(1,3))/2.d0+xk*vj2*xg(1,3)*xg(1,3)
985  & -xk*(vj2-1.d0)*(xg(1,1)*xg(3,3)+xg(1,3)*xg(3,1))/2.d0
986  elas(16)=(umb-f0*umbb)*(xg(1,2)*xg(1,3)+xg(1,3)*xg(1,2)-
987  & 2.d0*xg(1,1)*xg(2,3)/3.d0)
988  & -2.d0*(xs(1,1)*xg(2,3)+xg(1,1)*xs(2,3))/3.d0
989  & +f0*2.d0*(xx(1,1)*xg(2,3)+xg(1,1)*xx(2,3))/3.d0
990  & -d1*xn(1,1)*xn(2,3)-d2*(xn(1,1)*xd(2,3)+
991  & xd(1,1)*xn(2,3))/2.d0+xk*vj2*xg(1,1)*xg(2,3)
992  & -xk*(vj2-1.d0)*(xg(1,2)*xg(1,3)+xg(1,3)*xg(1,2))/2.d0
993  elas(17)=(umb-f0*umbb)*(xg(2,2)*xg(2,3)+xg(2,3)*xg(2,2)-
994  & 2.d0*xg(2,2)*xg(2,3)/3.d0)
995  & -2.d0*(xs(2,2)*xg(2,3)+xg(2,2)*xs(2,3))/3.d0
996  & +f0*2.d0*(xx(2,2)*xg(2,3)+xg(2,2)*xx(2,3))/3.d0
997  & -d1*xn(2,2)*xn(2,3)-d2*(xn(2,2)*xd(2,3)+
998  & xd(2,2)*xn(2,3))/2.d0+xk*vj2*xg(2,2)*xg(2,3)
999  & -xk*(vj2-1.d0)*(xg(2,2)*xg(2,3)+xg(2,3)*xg(2,2))/2.d0
1000  elas(18)=(umb-f0*umbb)*(xg(3,2)*xg(3,3)+xg(3,3)*xg(3,2)-
1001  & 2.d0*xg(3,3)*xg(2,3)/3.d0)
1002  & -2.d0*(xs(3,3)*xg(2,3)+xg(3,3)*xs(2,3))/3.d0
1003  & +f0*2.d0*(xx(3,3)*xg(2,3)+xg(3,3)*xx(2,3))/3.d0
1004  & -d1*xn(3,3)*xn(2,3)-d2*(xn(3,3)*xd(2,3)+
1005  & xd(3,3)*xn(2,3))/2.d0+xk*vj2*xg(3,3)*xg(2,3)
1006  & -xk*(vj2-1.d0)*(xg(3,2)*xg(3,3)+xg(3,3)*xg(3,2))/2.d0
1007  elas(19)=(umb-f0*umbb)*(xg(1,2)*xg(2,3)+xg(1,3)*xg(2,2)-
1008  & 2.d0*xg(1,2)*xg(2,3)/3.d0)
1009  & -2.d0*(xs(1,2)*xg(2,3)+xg(1,2)*xs(2,3))/3.d0
1010  & +f0*2.d0*(xx(1,2)*xg(2,3)+xg(1,2)*xx(2,3))/3.d0
1011  & -d1*xn(1,2)*xn(2,3)-d2*(xn(1,2)*xd(2,3)+
1012  & xd(1,2)*xn(2,3))/2.d0+xk*vj2*xg(1,2)*xg(2,3)
1013  & -xk*(vj2-1.d0)*(xg(1,2)*xg(2,3)+xg(1,3)*xg(2,2))/2.d0
1014  elas(20)=(umb-f0*umbb)*(xg(1,2)*xg(3,3)+xg(1,3)*xg(3,2)-
1015  & 2.d0*xg(1,3)*xg(2,3)/3.d0)
1016  & -2.d0*(xs(1,3)*xg(2,3)+xg(1,3)*xs(2,3))/3.d0
1017  & +f0*2.d0*(xx(1,3)*xg(2,3)+xg(1,3)*xx(2,3))/3.d0
1018  & -d1*xn(1,3)*xn(2,3)-d2*(xn(1,3)*xd(2,3)+
1019  & xd(1,3)*xn(2,3))/2.d0+xk*vj2*xg(1,3)*xg(2,3)
1020  & -xk*(vj2-1.d0)*(xg(1,2)*xg(3,3)+xg(1,3)*xg(3,2))/2.d0
1021  elas(21)=(umb-f0*umbb)*(xg(2,2)*xg(3,3)+xg(2,3)*xg(3,2)-
1022  & 2.d0*xg(2,3)*xg(2,3)/3.d0)
1023  & -2.d0*(xs(2,3)*xg(2,3)+xg(2,3)*xs(2,3))/3.d0
1024  & +f0*2.d0*(xx(2,3)*xg(2,3)+xg(2,3)*xx(2,3))/3.d0
1025  & -d1*xn(2,3)*xn(2,3)-d2*(xn(2,3)*xd(2,3)+
1026  & xd(2,3)*xn(2,3))/2.d0+xk*vj2*xg(2,3)*xg(2,3)
1027  & -xk*(vj2-1.d0)*(xg(2,2)*xg(3,3)+xg(2,3)*xg(3,2))/2.d0
1028 !
1029  endif
1030 !
1031 ! updating the plastic fields
1032 !
1033  do i=1,3
1034  cpl(i)=cpl(i)-1.d0
1035  enddo
1036  do i=1,6
1037  xstate(1+i,iint,iel)=-cpl(i)/2.d0
1038  xstate(7+i,iint,iel)=stbl(i)
1039  enddo
1040  xstate(1,iint,iel)=ep
1041 !
1042 ! maximum equivalent viscoplastic strain in this increment
1043 !
1044  depvisc=max(depvisc,ep-epini)
1045 !
1046  return
ident
subroutine ident(x, px, n, id)
Definition: ident.f:26
max
#define max(a, b)
Definition: cascade.c:32
c1
static double * c1
Definition: mafillvcompmain.c:30
creep
subroutine creep(decra, deswa, statev, serd, ec, esw, p, qtild, temp, dtemp, predef, dpred, time, dtime, cmname, leximp, lend, coords, nstatv, noel, npt, layer, kspt, kstep, kinc)
Definition: creep.f:22
uhardening
subroutine uhardening(amat, iel, iint, t1l, epini, ep, dtime, fiso, dfiso, fkin, dfkin)
Definition: uhardening.f:21
f1
static double * f1
Definition: objectivemain_se.c:47
Hosted by OpenAircraft.com, (Michigan UAV, LLC)