1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine incplas_lin(elconloc,plconloc,xstate,xstateini,
20     &  elas,emec,ithermal,icmd,beta,stre,vj,kode,
21     &  ielas,amat,t1l,dtime,time,ttime,iel,iint,nstate_,mi,
22     &  eloc,pgauss,nmethod,pnewdt,depvisc)
23!
24!     calculates stiffness and stresses for the incremental plasticity
25!     material law
26!
27!     icmd=3: calculates stress at mechanical strain
28!     else: calculates stress at mechanical strain and the stiffness
29!           matrix
30!
31!     This routine is meant for small strains. Reference:
32!     Dhondt, Guido; The Finite Element Method for Three-dimensional
33!     Thermomechanical Applications (Wiley, 2004), Section 5.3
34!
35      implicit none
36!
37      character*80 amat
38!
39      integer ithermal(*),icmd,i,k,l,m,n,kode,ivisco,ielastic,kel(4,21),
40     &  niso,nkin,ielas,iel,iint,nstate_,mi(*),id,leximp,lend,layer,
41     &  kspt,kstep,kinc,iloop,nmethod,user_hardening,user_creep
42!
43      real*8 elconloc(*),elas(21),emec(6),beta(6),stre(6),
44     &  vj,plconloc(802),stbl(6),stril(6),xitril(6),
45     &  ee,un,um,al,cop,dxitril,xn(3,3),epl(6),c1,c2,c3,c4,c7,
46     &  c8,ftrial,xiso(200),yiso(200),xkin(200),ykin(200),
47     &  fiso,dfiso,fkin,dfkin,fiso0,fkin0,ep,t1l,dtime,
48     &  epini,a1,dsvm,xxa,xxn,dkl(3,3),el(6),tracee,traces,
49     &  dcop,time,ttime,eloc(6),xstate(nstate_,mi(1),*),
50     &  xstateini(nstate_,mi(1),*),decra(5),deswa(5),serd,
51     &  esw(2),ec(2),p,qtild,predef(1),dpred(1),timeabq(2),pgauss(3),
52     &  dtemp,pnewdt,um2,depvisc
53!
54!
55      kel=reshape((/1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3,
56     &          1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,
57     &          3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,
58     &          1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/))
59      dkl=reshape((/1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,1.d0/),
60     &      (/3,3/))
61!
62      leximp=1
63      lend=2
64      user_creep=0
65!
66!     localizing the plastic fields
67!
68      do i=1,6
69         epl(i)=xstateini(1+i,iint,iel)
70         stbl(i)=xstateini(7+i,iint,iel)
71      enddo
72      epini=xstateini(1,iint,iel)
73!
74      ee=elconloc(1)
75      un=elconloc(2)
76      um2=ee/(1.d0+un)
77      al=um2*un/(1.d0-2.d0*un)
78      um=um2/2.d0
79!
80      ep=epini
81!
82!     check whether the user activated a viscous calculation
83!     (*VISCO instead of *STATIC)
84!
85      if((nmethod.ne.1).or.(ithermal(1).eq.3)) then
86         ivisco=1
87      else
88         ivisco=0
89      endif
90!
91!     check for user subroutines
92!
93      if((plconloc(801).lt.0.8d0).and.(plconloc(802).lt.0.8d0)) then
94         user_hardening=1
95      else
96         user_hardening=0
97      endif
98      if((kode.eq.-52).and.(ivisco.eq.1)) then
99         if(elconloc(3).lt.0.d0) then
100            user_creep=1
101         else
102            user_creep=0
103            xxa=elconloc(3)*(ttime+time)**elconloc(5)
104            if(xxa.lt.1.d-20) xxa=1.d-20
105            xxn=elconloc(4)
106            a1=xxa*dtime
107         endif
108      endif
109!
110!     trial elastic strain
111!
112      do i=1,6
113         el(i)=emec(i)-epl(i)
114      enddo
115!
116!     trial stress
117!
118      tracee=el(1)+el(2)+el(3)
119      do i=1,6
120         stre(i)=um2*el(i)-beta(i)
121      enddo
122      do i=1,3
123         stre(i)=stre(i)+al*tracee
124      enddo
125!
126!     trial deviatoric stress
127!
128      traces=(stre(1)+stre(2)+stre(3))/3.d0
129      do i=1,3
130         stril(i)=stre(i)-traces
131      enddo
132      do i=4,6
133         stril(i)=stre(i)
134      enddo
135!
136!     subtracting the back stress -> trial radius vector
137!
138      do i=1,6
139         xitril(i)=stril(i)-stbl(i)
140      enddo
141!
142!     size of trial radius vector
143!
144      dxitril=0.d0
145      do i=1,3
146         dxitril=dxitril+xitril(i)*xitril(i)
147      enddo
148      do i=4,6
149         dxitril=dxitril+2.d0*xitril(i)*xitril(i)
150      enddo
151      dxitril=dsqrt(dxitril)
152!
153!        restoring the hardening curves for the actual temperature
154!        plconloc contains the true stresses. By multiplying by
155!        the Jacobian, yiso and ykin are Kirchhoff stresses, as
156!        required by the hyperelastic theory (cf. Simo, 1988).
157!
158      niso=int(plconloc(801))
159      nkin=int(plconloc(802))
160      if(niso.ne.0) then
161         do i=1,niso
162            xiso(i)=plconloc(2*i-1)
163            yiso(i)=plconloc(2*i)
164         enddo
165      endif
166      if(nkin.ne.0) then
167         do i=1,nkin
168            xkin(i)=plconloc(399+2*i)
169            ykin(i)=plconloc(400+2*i)
170         enddo
171      endif
172!
173!     if no viscous calculation is performed a pure creep calculatino
174!     (without plasticity) is reduced to an elastic calculation
175!
176      ielastic=0
177      if(ivisco.eq.0) then
178         if(niso.eq.2) then
179            if((dabs(yiso(1)).lt.1.d-10).and.
180     &         (dabs(yiso(2)).lt.1.d-10)) then
181               ielastic=1
182            endif
183         endif
184      endif
185!
186!     check for yielding
187!
188      if(user_hardening.eq.1) then
189         call uhardening(amat,iel,iint,t1l,epini,ep,dtime,
190     &        fiso,dfiso,fkin,dfkin)
191      else
192         if(niso.ne.0) then
193            call ident(xiso,ep,niso,id)
194            if(id.eq.0) then
195               fiso=yiso(1)
196            elseif(id.eq.niso) then
197               fiso=yiso(niso)
198            else
199               dfiso=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
200               fiso=yiso(id)+dfiso*(ep-xiso(id))
201            endif
202         elseif(nkin.ne.0) then
203            fiso=ykin(1)
204         else
205            fiso=0.d0
206         endif
207      endif
208!
209      c1=2.d0/3.d0
210      c2=dsqrt(c1)
211!
212      ftrial=dxitril-c2*fiso
213      if((ftrial.le.1.d-10).or.(ielas.eq.1).or.(ielastic.eq.1)
214     &    .or.(dtime.lt.1.d-30)) then
215!
216!        updating the plastic fields
217!
218         do i=1,6
219            xstate(1+i,iint,iel)=epl(i)
220            xstate(7+i,iint,iel)=stbl(i)
221         enddo
222         xstate(1,iint,iel)=ep
223!
224         if(icmd.ne.3) then
225            elas(1)=al+um2
226            elas(2)=al
227            elas(3)=al+um2
228            elas(4)=al
229            elas(5)=al
230            elas(6)=al+um2
231            elas(7)=0.d0
232            elas(8)=0.d0
233            elas(9)=0.d0
234            elas(10)=um
235            elas(11)=0.d0
236            elas(12)=0.d0
237            elas(13)=0.d0
238            elas(14)=0.d0
239            elas(15)=um
240            elas(16)=0.d0
241            elas(17)=0.d0
242            elas(18)=0.d0
243            elas(19)=0.d0
244            elas(20)=0.d0
245            elas(21)=um
246         endif
247!
248         return
249      endif
250!
251!        plastic deformation
252!
253!        calculating the consistency parameter
254!
255      iloop=0
256      cop=0.d0
257!
258      loop: do
259         iloop=iloop+1
260         ep=epini+c2*cop
261!
262         if(user_hardening.eq.1) then
263            call uhardening(amat,iel,iint,t1l,epini,ep,dtime,
264     &           fiso,dfiso,fkin,dfkin)
265         else
266            if(niso.ne.0) then
267               call ident(xiso,ep,niso,id)
268               if(id.eq.0) then
269                  fiso=yiso(1)
270                  dfiso=0.d0
271               elseif(id.eq.niso) then
272                  fiso=yiso(niso)
273                  dfiso=0.d0
274               else
275                  dfiso=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
276                  fiso=yiso(id)+dfiso*(ep-xiso(id))
277               endif
278            elseif(nkin.ne.0) then
279               fiso=ykin(1)
280               dfiso=0.d0
281            else
282               fiso=0.d0
283               dfiso=0.d0
284            endif
285!
286            if(nkin.ne.0) then
287               call ident(xkin,ep,nkin,id)
288               if(id.eq.0) then
289                  fkin=ykin(1)
290                  dfkin=0.d0
291               elseif(id.eq.nkin) then
292                  fkin=ykin(nkin)
293                  dfkin=0.d0
294               else
295                  dfkin=(ykin(id+1)-ykin(id))/(xkin(id+1)-xkin(id))
296                  fkin=ykin(id)+dfkin*(ep-xkin(id))
297               endif
298            elseif(niso.ne.0) then
299               fkin=yiso(1)
300               dfkin=0.d0
301            else
302               fkin=0.d0
303               dfkin=0.d0
304            endif
305         endif
306!
307         if(dabs(cop).lt.1.d-10) then
308            fiso0=fiso
309            fkin0=fkin
310         endif
311!
312         if((kode.eq.-51).or.(ivisco.eq.0)) then
313            dcop=(dxitril-c2*(fiso+fkin-fkin0)-um2*cop)/
314     &           (um2+c1*(dfiso+dfkin))
315         else
316            if(user_creep.eq.1) then
317               if(ithermal(1).eq.0) then
318                  write(*,*) '*ERROR in incplas: no temperature defined'
319                  call exit(201)
320               endif
321               timeabq(1)=time
322               timeabq(2)=ttime+time
323!
324!              Eqn. (5.139)
325!
326               qtild=(dxitril-um2*cop)/c2-fiso-(fkin-fkin0)
327!
328!              the Von Mises stress must be positive
329!
330               if(qtild.lt.1.d-10) qtild=1.d-10
331               ec(1)=epini
332               call creep(decra,deswa,xstateini(1,iint,iel),serd,ec,
333     &             esw,p,qtild,t1l,dtemp,predef,dpred,timeabq,dtime,
334     &             amat,leximp,lend,pgauss,nstate_,iel,iint,layer,kspt,
335     &             kstep,kinc)
336               dsvm=1.d0/decra(5)
337               dcop=(decra(1)-c2*cop)/
338     &                 (c2*(decra(5)*(dfiso+um*(3.d0+dfkin/um))+1.d0))
339            else
340               qtild=(dxitril-um2*cop)/c2-fiso-(fkin-fkin0)
341!
342!              the Von Mises stress must be positive
343!
344               if(qtild.lt.1.d-10) qtild=1.d-10
345               decra(1)=a1*qtild**xxn
346               decra(5)=xxn*decra(1)/qtild
347               dsvm=1.d0/decra(5)
348               dcop=(decra(1)-c2*cop)/
349     &                 (c2*(decra(5)*(dfiso+um*(3.d0+dfkin/um))+1.d0))
350            endif
351         endif
352         cop=cop+dcop
353!
354         if((dabs(dcop).lt.cop*1.d-4).or.
355     &      (dabs(dcop).lt.1.d-10)) exit
356!
357!        check for endless loops or a negative consistency
358!        parameter
359!
360         if((iloop.gt.15).or.(cop.le.0.d0)) then
361            pnewdt=0.25d0
362            return
363         endif
364!
365      enddo loop
366!
367!        updating the equivalent plastic strain
368!
369      ep=epini+c2*cop
370!
371!     updating the back stress
372!
373      c7=c2*(fkin-fkin0)/dxitril
374      do i=1,6
375         stbl(i)=stbl(i)-c7*xitril(i)
376      enddo
377!
378!     updating the plastic strain tensor
379!
380      c7=cop/dxitril
381      do i=1,6
382         epl(i)=epl(i)+c7*xitril(i)
383      enddo
384!
385!     trial elastic strain
386!
387      do i=1,6
388         el(i)=emec(i)-epl(i)
389      enddo
390!
391!     trial stress
392!
393      tracee=el(1)+el(2)+el(3)
394      do i=1,6
395         stre(i)=um2*el(i)-beta(i)
396      enddo
397      do i=1,3
398         stre(i)=stre(i)+al*tracee
399      enddo
400!
401!     calculating the local stiffness matrix
402!
403      if(icmd.ne.3) then
404         c7=um2**2*cop/dxitril
405!
406         if((kode.eq.-51).or.(ivisco.eq.0)) then
407            c8=um2**2/(um2+c1*(dfiso+dfkin))
408         else
409            c8=um2**2/(um2+c1*(dfiso+dfkin+1.d0/decra(5)))
410         endif
411!
412         c1=um-c7/2.d0
413         c2=al+c7/3.d0
414         c3=(c7-c8)/(dxitril**2)
415!
416         xn(1,1)=xitril(1)
417         xn(2,2)=xitril(2)
418         xn(3,3)=xitril(3)
419         xn(1,2)=xitril(4)
420         xn(1,3)=xitril(5)
421         xn(2,3)=xitril(6)
422         xn(2,1)=xn(1,2)
423         xn(3,1)=xn(1,3)
424         xn(3,2)=xn(2,3)
425!
426         do i=1,21
427            k=kel(1,i)
428            l=kel(2,i)
429            m=kel(3,i)
430            n=kel(4,i)
431            elas(i)=c1*(dkl(k,m)*dkl(l,n)+dkl(k,n)*dkl(l,m))
432     &             +c2*dkl(k,l)*dkl(m,n)
433     &             +c3*xn(k,l)*xn(m,n)
434         enddo
435      endif
436!
437!        updating the plastic fields
438!
439      do i=1,6
440         xstate(1+i,iint,iel)=epl(i)
441         xstate(7+i,iint,iel)=stbl(i)
442      enddo
443      xstate(1,iint,iel)=ep
444!
445!     maximum difference in the equivalent viscoplastic strain
446!     in this increment based on the viscoplastic strain rate at
447!     the start and the end of the increment
448!
449      if(ivisco.eq.1) then
450         depvisc=max(depvisc,abs(decra(1)-dtime*xstateini(14,iint,iel)))
451         xstate(14,iint,iel)=decra(1)/dtime
452      endif
453!
454!     maximum difference in the equivalent viscoplastic strain
455!     in this increment based on the viscoplastic strain rate at
456!     the start and the end of the increment
457!
458      if(ivisco.eq.1) then
459         depvisc=max(depvisc,abs(decra(1)-dtime*xstateini(14,iint,iel)))
460         xstate(14,iint,iel)=decra(1)/dtime
461      endif
462c      depvisc=max(depvisc,ep-epini)
463!
464      return
465      end
466