CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
resultstherm.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine resultstherm (co, kon, ipkon, lakon, v, elcon, nelcon, rhcon, nrhcon, ielmat, ielorien, norien, orab, ntmat_, t0, iperturb, fn, shcon, nshcon, iout, qa, vold, ipompc, nodempc, coefmpc, nmpc, dtime, time, ttime, plkcon, nplkcon, xstateini, xstiff, xstate, npmat_, matname, mi, ncmat_, nstate_, cocon, ncocon, qfx, ikmpc, ilmpc, istep, iinc, springarea, calcul_fn, calcul_qa, nal, nea, neb, ithermal, nelemload, nload, nmethod, reltime, sideload, xload, xloadold, pslavsurf, pmastsurf, mortar, clearini, plicon, nplicon, ielprop, prop, iponoel, inoel, network, ipobody, xbody, ibody)
 

Function/Subroutine Documentation

◆ resultstherm()

subroutine resultstherm ( real*8, dimension(3,*)  co,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
real*8, dimension(0:mi(2),*)  v,
real*8, dimension(0:ncmat_,ntmat_,*)  elcon,
integer, dimension(2,*)  nelcon,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
integer, dimension(mi(3),*)  ielmat,
integer, dimension(mi(3),*)  ielorien,
integer  norien,
real*8, dimension(7,*)  orab,
integer  ntmat_,
real*8, dimension(*)  t0,
integer, dimension(*)  iperturb,
real*8, dimension(0:mi(2),*)  fn,
real*8, dimension(0:3,ntmat_,*)  shcon,
integer, dimension(*)  nshcon,
integer  iout,
real*8, dimension(*)  qa,
real*8, dimension(0:mi(2),*)  vold,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
real*8  dtime,
real*8  time,
real*8  ttime,
real*8, dimension(0:2*npmat_,ntmat_,*)  plkcon,
integer, dimension(0:ntmat_,*)  nplkcon,
real*8, dimension(nstate_,mi(1),*)  xstateini,
real*8, dimension(27,mi(1),*)  xstiff,
real*8, dimension(nstate_,mi(1),*)  xstate,
integer  npmat_,
character*80, dimension(*)  matname,
integer, dimension(*)  mi,
integer  ncmat_,
integer  nstate_,
real*8, dimension(0:6,ntmat_,*)  cocon,
integer, dimension(2,*)  ncocon,
real*8, dimension(3,mi(1),*)  qfx,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
integer  istep,
integer  iinc,
real*8, dimension(2,*)  springarea,
integer  calcul_fn,
integer  calcul_qa,
integer  nal,
integer  nea,
integer  neb,
integer, dimension(2)  ithermal,
integer, dimension(2,*)  nelemload,
integer  nload,
integer  nmethod,
real*8  reltime,
character*20, dimension(*)  sideload,
real*8, dimension(2,*)  xload,
real*8, dimension(2,*)  xloadold,
real*8, dimension(3,*)  pslavsurf,
real*8, dimension(6,*)  pmastsurf,
integer  mortar,
real*8, dimension(3,9,*)  clearini,
real*8, dimension(0:2*npmat_,ntmat_,*)  plicon,
integer, dimension(0:ntmat_,*)  nplicon,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop,
integer, dimension(*)  iponoel,
integer, dimension(2,*)  inoel,
integer  network,
integer, dimension(2,*)  ipobody,
real*8, dimension(7,*)  xbody,
integer, dimension(3,*)  ibody 
)
30 !
31 ! calculates the heat flux and the material tangent at the integration
32 ! points and the internal concentrated flux at the nodes
33 !
34  implicit none
35 !
36  character*8 lakon(*),lakonl
37  character*20 sideload(*)
38  character*80 amat,matname(*)
39 !
40  integer kon(*),konl(26),iperm(20),ikmpc(*),ilmpc(*),mi(*),
41  & nelcon(2,*),nrhcon(*),ielmat(mi(3),*),ielorien(mi(3),*),
42  & ntmat_,ipkon(*),ipompc(*),nodempc(3,*),mortar,igauss,
43  & ncocon(2,*),iflag,nshcon(*),istep,iinc,mt,mattyp,
44  & i,j,k,m1,kk,i1,m3,indexe,nope,norien,iperturb(*),iout,
45  & nal,nmpc,kode,imat,mint3d,iorien,istiff,ncmat_,nstate_,
46  & nplkcon(0:ntmat_,*),npmat_,calcul_fn,calcul_qa,nea,neb,
47  & nelemload(2,*),nload,ithermal(2),nmethod,nopered,iloc,
48  & jfaces,node,nplicon(0:ntmat_,*),null,ielprop(*),
49  & iponoel(*),inoel(2,*),network,ipobody(2,*),ibody(3,*)
50 !
51  real*8 co(3,*),v(0:mi(2),*),shp(4,26),reltime,
52  & xl(3,26),vl(0:mi(2),26),elcon(0:ncmat_,ntmat_,*),
53  & rhcon(0:1,ntmat_,*),qfx(3,mi(1),*),orab(7,*),
54  & rho,fn(0:mi(2),*),tnl(19),timeend(2),q(0:mi(2),26),
55  & vkl(0:3,3),t0(*),vold(0:mi(2),*),coefmpc(*),
56  & springarea(2,*),elconloc(21),cocon(0:6,ntmat_,*),
57  & shcon(0:3,ntmat_,*),sph,c1,xi,et,ze,xsj,qa(*),t0l,t1l,dtime,
58  & weight,pgauss(3),coconloc(6),qflux(3),time,ttime,
59  & t1lold,plkcon(0:2*npmat_,ntmat_,*),xstiff(27,mi(1),*),
60  & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
61  & xload(2,*),xloadold(2,*),clearini(3,9,*),pslavsurf(3,*),
62  & pmastsurf(6,*),plicon(0:2*npmat_,ntmat_,*),prop(*),
63  & xbody(7,*)
64 !
65  include "gauss.f"
66 !
67  iflag=3
68  null=0
69  iperm=(/5,6,7,8,1,2,3,4,13,14,15,16,9,10,11,12,17,18,19,20/)
70 !
71  mt=mi(2)+1
72 !
73 ! calculation of temperatures and thermal flux
74 !
75  nal=0
76 !
77  do i=nea,neb
78 !
79  if(ipkon(i).lt.0) cycle
80 !
81 ! no 3D-fluid elements
82 !
83  if(lakon(i)(1:1).eq.'F') cycle
84  if(lakon(i)(1:7).eq.'DCOUP3D') cycle
85 !
86  imat=ielmat(1,i)
87  amat=matname(imat)
88  if(norien.gt.0) then
89  iorien=ielorien(1,i)
90  else
91  iorien=0
92  endif
93 !
94  indexe=ipkon(i)
95  if(lakon(i)(4:5).eq.'20') then
96  nope=20
97  elseif(lakon(i)(4:4).eq.'8') then
98  nope=8
99  elseif(lakon(i)(4:5).eq.'10') then
100  nope=10
101  elseif(lakon(i)(4:4).eq.'4') then
102  nope=4
103  elseif(lakon(i)(4:5).eq.'15') then
104  nope=15
105  elseif(lakon(i)(4:4).eq.'6') then
106  nope=6
107  elseif((lakon(i)(1:2).eq.'ES').and.(lakon(i)(7:7).ne.'A')) then
108 !
109 ! contact spring and advection elements (no dashpot elements
110 ! = ED... elements and no genuine spring elements)
111 !
112  if(lakon(i)(7:7).eq.'C') then
113 !
114 ! contact spring elements
115 !
116  if(mortar.eq.1) then
117 !
118 ! face-to-face penalty
119 !
120  nope=kon(ipkon(i))
121  elseif(mortar.eq.0) then
122 !
123 ! node-to-face penalty
124 !
125  nope=ichar(lakon(i)(8:8))-47
126  konl(nope+1)=kon(indexe+nope+1)
127  endif
128  else
129 !
130 ! advection elements
131 !
132  nope=ichar(lakon(i)(8:8))-47
133  endif
134 c nope=ichar(lakon(i)(8:8))-47
135 c!
136 c! local contact spring number
137 c!
138 c if(lakon(i)(7:7).eq.'C') konl(nope+1)=kon(indexe+nope+1)
139  elseif((lakon(i)(1:2).eq.'D ').or.
140  & ((lakon(i)(1:1).eq.'D').and.(network.eq.1))) then
141 !
142 ! no entry or exit elements
143 !
144  if((kon(indexe+1).eq.0).or.(kon(indexe+3).eq.0)) cycle
145  nope=3
146  else
147  cycle
148  endif
149 !
150  if(lakon(i)(4:5).eq.'8R') then
151  mint3d=1
152  elseif(lakon(i)(4:7).eq.'20RB') then
153  if((lakon(i)(8:8).eq.'R').or.(lakon(i)(8:8).eq.'C')) then
154  mint3d=50
155  else
156  call beamintscheme(lakon(i),mint3d,ielprop(i),prop,
157  & null,xi,et,ze,weight)
158  endif
159  elseif((lakon(i)(4:4).eq.'8').or.
160  & (lakon(i)(4:6).eq.'20R')) then
161  if(lakon(i)(6:7).eq.'RA') then
162  mint3d=4
163  else
164  mint3d=8
165  endif
166  elseif(lakon(i)(4:4).eq.'2') then
167  mint3d=27
168  elseif(lakon(i)(4:5).eq.'10') then
169  mint3d=4
170  elseif(lakon(i)(4:4).eq.'4') then
171  mint3d=1
172  elseif(lakon(i)(4:5).eq.'15') then
173  mint3d=9
174  elseif(lakon(i)(4:4).eq.'6') then
175  mint3d=2
176  else
177  mint3d=0
178  endif
179 !
180  do j=1,nope
181  konl(j)=kon(indexe+j)
182  do k=1,3
183  xl(k,j)=co(k,konl(j))
184  vl(k,j)=v(k,konl(j))
185  enddo
186  vl(0,j)=v(0,konl(j))
187  enddo
188 !
189 ! q contains the nodal forces per element; initialisation of q
190 !
191  if((iperturb(1).ge.2).or.((iperturb(1).le.0).and.(iout.lt.1)))
192  & then
193  do m1=1,nope
194  q(0,m1)=fn(0,konl(m1))
195  enddo
196  endif
197 !
198 ! calculating the concentrated flux for the contact elements
199 !
200  if(mint3d.eq.0) then
201 !
202  lakonl=lakon(i)
203 !
204 ! initialization of tnl
205 !
206  do j=1,nope
207  tnl(j)=0.d0
208  enddo
209 !
210 ! spring elements (including contact springs)
211 !
212  if(lakonl(2:2).eq.'S') then
213  if(lakonl(7:7).eq.'C') then
214 !
215 ! contact element
216 !
217  kode=nelcon(1,imat)
218  if(kode.eq.-51) then
219  timeend(1)=time
220  timeend(2)=ttime+time
221  if(mortar.eq.0) then
222  call springforc_n2f_th(xl,vl,imat,elcon,nelcon,
223  & tnl,ncmat_,ntmat_,nope,kode,elconloc,
224  & plkcon,nplkcon,npmat_,mi,
225  & springarea(1,konl(nope+1)),timeend,matname,
226  & konl(nope),i,istep,iinc,iperturb)
227  elseif(mortar.eq.1) then
228  iloc=kon(indexe+nope+1)
229  jfaces=kon(indexe+nope+2)
230  igauss=kon(indexe+nope+1)
231  node=0
232  call springforc_f2f_th(xl,vl,imat,elcon,nelcon,
233  & tnl,ncmat_,ntmat_,nope,lakonl,kode,elconloc,
234  & plicon,nplicon,npmat_,mi,springarea(1,iloc),
235  & nmethod,reltime,jfaces,igauss,
236  & pslavsurf,pmastsurf,clearini,timeend,istep,
237  & iinc,plkcon,nplkcon,node,i,matname)
238  endif
239  endif
240  elseif(lakonl(7:7).eq.'F') then
241 !
242 ! advective element
243 !
244  call advecforc(nope,vl,ithermal,xl,nelemload,
245  & i,nload,lakon,xload,istep,time,ttime,
246  & dtime,sideload,v,mi,xloadold,reltime,nmethod,
247  & tnl,iinc,iponoel,inoel,ielprop,prop,ielmat,shcon,
248  & nshcon,rhcon,nrhcon,ntmat_,ipkon,kon,cocon,
249  & ncocon,ipobody,xbody,ibody)
250  endif
251 !
252  elseif((lakonl(1:2).eq.'D ').or.
253  & ((lakonl(1:1).eq.'D').and.(network.eq.1))) then
254 !
255 ! generic networkelement
256 !
257  call networkforc(vl,tnl,imat,konl,mi,ntmat_,shcon,
258  & nshcon,rhcon,nrhcon)
259 
260  endif
261 !
262  do j=1,nope
263  fn(0,konl(j))=fn(0,konl(j))+tnl(j)
264  enddo
265  endif
266 !
267  do kk=1,mint3d
268  if(lakon(i)(4:5).eq.'8R') then
269  xi=gauss3d1(1,kk)
270  et=gauss3d1(2,kk)
271  ze=gauss3d1(3,kk)
272  weight=weight3d1(kk)
273  elseif(lakon(i)(4:7).eq.'20RB') then
274  if((lakon(i)(8:8).eq.'R').or.(lakon(i)(8:8).eq.'C')) then
275  xi=gauss3d13(1,kk)
276  et=gauss3d13(2,kk)
277  ze=gauss3d13(3,kk)
278  weight=weight3d13(kk)
279  else
280  call beamintscheme(lakon(i),mint3d,ielprop(i),
281  & prop,kk,xi,et,ze,weight)
282  endif
283  elseif((lakon(i)(4:4).eq.'8').or.
284  & (lakon(i)(4:6).eq.'20R'))
285  & then
286  xi=gauss3d2(1,kk)
287  et=gauss3d2(2,kk)
288  ze=gauss3d2(3,kk)
289  weight=weight3d2(kk)
290  elseif(lakon(i)(4:4).eq.'2') then
291  xi=gauss3d3(1,kk)
292  et=gauss3d3(2,kk)
293  ze=gauss3d3(3,kk)
294  weight=weight3d3(kk)
295  elseif(lakon(i)(4:5).eq.'10') then
296  xi=gauss3d5(1,kk)
297  et=gauss3d5(2,kk)
298  ze=gauss3d5(3,kk)
299  weight=weight3d5(kk)
300  elseif(lakon(i)(4:4).eq.'4') then
301  xi=gauss3d4(1,kk)
302  et=gauss3d4(2,kk)
303  ze=gauss3d4(3,kk)
304  weight=weight3d4(kk)
305  elseif(lakon(i)(4:5).eq.'15') then
306  xi=gauss3d8(1,kk)
307  et=gauss3d8(2,kk)
308  ze=gauss3d8(3,kk)
309  weight=weight3d8(kk)
310  elseif(lakon(i)(4:4).eq.'6') then
311  xi=gauss3d7(1,kk)
312  et=gauss3d7(2,kk)
313  ze=gauss3d7(3,kk)
314  weight=weight3d7(kk)
315  endif
316 !
317  if(nope.eq.20) then
318  if(lakon(i)(7:7).eq.'A') then
319  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
320  elseif((lakon(i)(7:7).eq.'E').or.
321  & (lakon(i)(7:7).eq.'S')) then
322  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
323  else
324  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
325  endif
326  elseif(nope.eq.8) then
327  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
328  elseif(nope.eq.10) then
329  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
330  elseif(nope.eq.4) then
331  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
332  elseif(nope.eq.15) then
333  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
334  else
335  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
336  endif
337  c1=xsj*weight
338 !
339 ! vkl(m2,m3) contains the derivative of the m2-
340 ! component of the displacement with respect to
341 ! direction m3
342 !
343  do m3=1,3
344  vkl(0,m3)=0.d0
345  enddo
346 !
347  do m1=1,nope
348  do m3=1,3
349  vkl(0,m3)=vkl(0,m3)+shp(m3,m1)*vl(0,m1)
350  enddo
351  enddo
352 !
353  kode=ncocon(1,imat)
354 !
355 ! calculating the temperature difference in
356 ! the integration point
357 !
358  t1lold=0.d0
359  t1l=0.d0
360  if((lakon(i)(4:5).eq.'8 ').or.
361  & (lakon(i)(4:5).eq.'8I')) then
362  do i1=1,8
363  t1lold=t1lold+vold(0,konl(i1))/8.d0
364  t1l=t1l+v(0,konl(i1))/8.d0
365  enddo
366  elseif(lakon(i)(4:6).eq.'20 ') then
367  nopered=20
368  call lintemp_th(t0,vold,konl,nopered,kk,t0l,t1lold,mi)
369  call lintemp_th(t0,v,konl,nopered,kk,t0l,t1l,mi)
370  elseif(lakon(i)(4:6).eq.'10T') then
371  call linscal10(vold,konl,t1lold,mi(2),shp)
372  call linscal10(v,konl,t1l,mi(2),shp)
373  else
374  do i1=1,nope
375  t1lold=t1lold+shp(4,i1)*vold(0,konl(i1))
376  t1l=t1l+shp(4,i1)*v(0,konl(i1))
377  enddo
378  endif
379 !
380 ! calculating the coordinates of the integration point
381 ! for material orientation purposes (for cylindrical
382 ! coordinate systems)
383 !
384  if((iorien.gt.0).or.(kode.le.-100)) then
385  do j=1,3
386  pgauss(j)=0.d0
387  do i1=1,nope
388  pgauss(j)=pgauss(j)+shp(4,i1)*co(j,konl(i1))
389  enddo
390  enddo
391  endif
392 !
393 ! material data; for linear elastic materials
394 ! this includes the calculation of the stiffness
395 ! matrix
396 !
397  istiff=0
398 !
399  call materialdata_th(cocon,ncocon,imat,iorien,pgauss,orab,
400  & ntmat_,coconloc,mattyp,t1l,rhcon,nrhcon,rho,shcon,
401  & nshcon,sph,xstiff,kk,i,istiff,mi(1))
402 !
403  call thermmodel(amat,i,kk,kode,coconloc,vkl,dtime,
404  & time,ttime,mi(1),nstate_,xstateini,xstate,qflux,xstiff,
405  & iorien,pgauss,orab,t1l,t1lold,vold,co,lakon(i),konl,
406  & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,nmethod,
407  & iperturb)
408 !
409  qfx(1,kk,i)=qflux(1)
410  qfx(2,kk,i)=qflux(2)
411  qfx(3,kk,i)=qflux(3)
412  if(lakon(i)(6:7).eq.'RA') then
413  qfx(1,kk+4,i)=qflux(1)
414  qfx(2,kk+4,i)=qflux(2)
415  qfx(3,kk+4,i)=qflux(3)
416  endif
417 !
418 ! calculation of the nodal flux
419 !
420  if(calcul_fn.eq.1)then
421 !
422 ! calculating fn using skl
423 !
424  if(lakon(i)(6:7).eq.'RA') then
425  do m1=1,nope
426  fn(0,konl(m1))=fn(0,konl(m1))
427  & -c1*(qflux(1)*(shp(1,m1)+shp(1,iperm(m1)))
428  & +qflux(2)*(shp(2,m1)+shp(2,iperm(m1)))
429  & +qflux(3)*(shp(3,m1)+shp(3,iperm(m1))))
430  enddo
431  else
432  do m1=1,nope
433  do m3=1,3
434  fn(0,konl(m1))=fn(0,konl(m1))-
435  & c1*qflux(m3)*shp(m3,m1)
436  enddo
437  enddo
438  endif
439  endif
440  enddo
441 !
442 ! q contains the contributions to the nodal force in the nodes
443 ! belonging to the element at stake from other elements (elements
444 ! already treated). These contributions have to be
445 ! subtracted to get the contributions attributable to the element
446 ! at stake only
447 !
448  if(calcul_qa.eq.1) then
449  do m1=1,nope
450  qa(2)=qa(2)+dabs(fn(0,konl(m1))-q(0,m1))
451  enddo
452  nal=nal+nope
453  endif
454  enddo
455 !
456  return
subroutine shape20h_pl(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h_pl.f:20
subroutine shape6w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape6w.f:20
subroutine linscal10(scal, konl, scall, idim, shp)
Definition: linscal10.f:20
subroutine springforc_n2f_th(xl, vl, imat, elcon, nelcon, tnl, ncmat_, ntmat_, nope, kode, elconloc, plkcon, nplkcon, npmat_, mi, springarea, timeend, matname, node, noel, istep, iinc, iperturb)
Definition: springforc_n2f_th.f:23
subroutine networkforc(vl, tnl, imat, konl, mi, ntmat_, shcon, nshcon, rhcon, nrhcon)
Definition: networkforc.f:21
subroutine shape10tet(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape10tet.f:20
subroutine shape8h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape8h.f:20
subroutine springforc_f2f_th(xl, vl, imat, elcon, nelcon, tnl, ncmat_, ntmat_, nope, lakonl, kode, elconloc, plicon, nplicon, npmat_, mi, springarea, nmethod, reltime, jfaces, igauss, pslavsurf, pmastsurf, clearini, timeend, istep, iinc, plkcon, nplkcon, node, noel, matname)
Definition: springforc_f2f_th.f:25
subroutine advecforc(nope, voldl, ithermal, xl, nelemload, nelemadvec, nload, lakon, xload, istep, time, ttime, dtime, sideload, vold, mi, xloadold, reltime, nmethod, tnl, iinc, iponoel, inoel, ielprop, prop, ielmat, shcon, nshcon, rhcon, nrhcon, ntmat_, ipkon, kon, cocon, ncocon, ipobody, ibody, xbody)
Definition: advecforc.f:24
subroutine shape15w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape15w.f:20
static double * c1
Definition: mafillvcompmain.c:30
subroutine shape20h_ax(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h_ax.f:20
subroutine materialdata_th(cocon, ncocon, imat, iorien, pgauss, orab, ntmat_, coconloc, mattyp, t1l, rhcon, nrhcon, rho, shcon, nshcon, sph, xstiff, iint, iel, istiff, mi)
Definition: materialdata_th.f:21
subroutine lintemp_th(t0, vold, konl, nope, jj, t0l, t1l, mi)
Definition: lintemp_th.f:20
subroutine shape20h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h.f:20
subroutine thermmodel(amat, iel, iint, kode, coconloc, vkl, dtime, time, ttime, mi, nstate_, xstateini, xstate, qflux, xstiff, iorien, pgauss, orab, t1l, t1lold, vold, co, lakonl, konl, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, nmethod, iperturb)
Definition: thermmodel.f:23
subroutine beamintscheme(lakonl, mint3d, npropstart, prop, kk, xi, et, ze, weight)
Definition: beamintscheme.f:21
subroutine shape4tet(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape4tet.f:20
static ITG * nal
Definition: results.c:31
Hosted by OpenAircraft.com, (Michigan UAV, LLC)