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

Go to the source code of this file.

Functions/Subroutines

subroutine incplas_lin (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_lin()

subroutine incplas_lin ( 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
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(21),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.,0.,0.,0.,1.,0.,0.,0.,1./),(/3,3/))
60 !
61 c pnewdt=-1.d0
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.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.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 equivalent viscoplastic strain in this increment
446 !
447  depvisc=max(depvisc,ep-epini)
448 !
449  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
Hosted by OpenAircraft.com, (Michigan UAV, LLC)