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 frictionheating(ne0,ne,ipkon,lakon,ielmat,mi,elcon,
20     &  ncmat_,ntmat_,kon,islavsurf,pmastsurf,springarea,co,vold,
21     &  veold,pslavsurf,xloadact,nload,nload_,nelemload,iamload,
22     &  idefload,sideload,stx,nam,time,ttime,matname,istep,iinc)
23!
24!     determines the effect of friction heating
25!
26      implicit none
27!
28      character*8 lakon(*),lakonl
29      character*20 label,sideload(*)
30      character*80 matname(*)
31!
32      integer i,j,k,ne0,ne,indexe,ipkon(*),imat,mi(*),ielmat(mi(3),*),
33     &  ncmat_,ntmat_,kon(*),nope,igauss,jfaces,ifaces,nelems,ifacem,
34     &  nelemm,jfacem,islavsurf(2,*),nopes,nopem,konl(20),iflag,
35     &  mint2d,iamplitude,isector,nload,nload_,nelemload(2,*),nam,
36     &  iamload(2,*),idefload(*),istep,iinc
37!
38      real*8 elcon(0:ncmat_,ntmat_,*),pressure,stx(6,mi(1),*),
39     &  pmastsurf(6,*),area,springarea(2,*),pl(3,20),co(3,*),
40     &  vold(0:mi(2),*),areaslav,xi,et,vels(3),veold(0:mi(2),*),
41     &  xsj2m(3),xs2m(3,7),shp2m(7,9),xsj2s(3),xs2s(3,7),shp2s(7,9),
42     &  areamast,pslavsurf(3,*),value,velm(3),um,xloadact(2,*),weight,
43     &  shear,vnorm,f,eta,timeend(2),time,ttime,coords(3),xl(3,20)
44!
45!
46!
47      include "gauss.f"
48!
49!     mortar=1 is assumed (face-to-face penalty contact)
50!     ithermal=3 is assumed
51!
52      iamplitude=0
53      isector=0
54!
55      do i=ne0+1,ne
56         imat=ielmat(1,i)
57!
58!        heat conversion factor
59!
60         eta=elcon(9,1,imat)
61!
62!        surface weighting factor
63!
64         f=elcon(10,1,imat)
65!
66!        velocity
67!
68         vnorm=elcon(11,1,imat)
69!
70!        friction coefficient
71!
72         um=elcon(6,1,imat)
73!
74         pressure=stx(4,1,i)
75         if(pressure.lt.0.d0) cycle
76!
77         shear=dsqrt(stx(5,1,i)**2+stx(6,1,i)**2)
78         if(vnorm.lt.-0.5d0) then
79!
80!           if ||v||<0 => take differential velocity from the results
81!           no heat generation if no slip
82!
83            if(shear.lt.um*pressure*0.95d0) cycle
84         endif
85!
86         indexe=ipkon(i)
87         lakonl=lakon(i)
88!
89         nope=kon(ipkon(i))
90         nopem=ichar(lakonl(8:8))-48
91         nopes=nope-nopem
92!
93         igauss=kon(indexe+nope+1)
94         jfaces=kon(indexe+nope+2)
95!
96!        slave face
97!
98         ifaces=islavsurf(1,jfaces)
99         nelems=int(ifaces/10.d0)
100         jfaces=ifaces-10*nelems
101!
102!        master face
103!
104         ifacem=int(pmastsurf(3,igauss))
105         nelemm=int(ifacem/10.d0)
106         jfacem=ifacem-10*nelemm
107!
108!        contact area
109!
110         area=springarea(1,igauss)
111!
112!        slave and master nodes
113!
114         do j=1,nope
115            konl(j)=kon(indexe+j)
116            do k=1,3
117               pl(k,j)=co(k,konl(j))+vold(k,konl(j))
118            enddo
119         enddo
120!
121!        user subroutine called if vnorm=-0.01d0
122!
123         if((vnorm.lt.0.d0).and.(vnorm.gt.-0.5d0)) then
124!
125            xi=pslavsurf(1,igauss)
126            et=pslavsurf(2,igauss)
127!
128            do j=nopem+1,nopem+nopes
129               konl(j)=kon(indexe+j)
130               do k=1,3
131                  xl(k,j)=co(k,konl(j))
132               enddo
133            enddo
134!
135!           determining the jacobian vector on the surface
136!
137            iflag=1
138            if(nopes.eq.8) then
139               call shape8q(xi,et,xl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
140            elseif(nopes.eq.4) then
141               call shape4q(xi,et,xl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
142            elseif(nopes.eq.6) then
143               call shape6tri(xi,et,xl(1,nopem+1),xsj2s,xs2s,shp2s,
144     &                        iflag)
145            else
146               call shape3tri(xi,et,xl(1,nopem+1),xsj2s,xs2s,shp2s,
147     &                        iflag)
148            endif
149!
150!           position of the slave integration point
151!
152            do j=1,3
153               coords(j)=0.d0
154               do k=1,nopes
155                  coords(j)=coords(j)+shp2s(4,k)*xl(j,nopem+k)
156               enddo
157            enddo
158!
159            timeend(1)=time
160            timeend(2)=ttime+time
161            call fricheat(eta,f,vnorm,timeend,matname(imat),i,
162     &                    nelems,jfaces,nelemm,jfacem,um,
163     &                    istep,iinc,area,pressure,coords)
164         endif
165!
166!        heat flux into the slave face
167!
168         if(nopes.eq.8) then
169            mint2d=9
170         elseif(nopes.eq.6) then
171            mint2d=3
172         elseif(nopes.eq.4) then
173            mint2d=4
174         else
175            mint2d=1
176         endif
177!
178!        calculating the area of the slave face
179!
180         areaslav=0.d0
181!
182         do j=1,mint2d
183            if(nopes.eq.8) then
184               xi=gauss2d3(1,j)
185               et=gauss2d3(2,j)
186               weight=weight2d3(j)
187            elseif(nopes.eq.6) then
188               xi=gauss2d5(1,j)
189               et=gauss2d5(2,j)
190               weight=weight2d5(j)
191            elseif(nopes.eq.4) then
192               xi=gauss2d2(1,j)
193               et=gauss2d2(2,j)
194               weight=weight2d2(j)
195            else
196               xi=gauss2d4(1,j)
197               et=gauss2d4(2,j)
198               weight=weight2d4(j)
199            endif
200!
201            iflag=2
202            if(nopes.eq.8) then
203               call shape8q(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
204            elseif(nopes.eq.4) then
205               call shape4q(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
206            elseif(nopes.eq.6) then
207               call shape6tri(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,
208     &iflag)
209            else
210               call shape3tri(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,
211     &iflag)
212            endif
213!
214            areaslav=areaslav+dsqrt(xsj2s(1)**2+
215     &                              xsj2s(2)**2+
216     &                              xsj2s(3)**2)
217         enddo
218!
219         label(1:20)='S                   '
220         write(label(2:2),'(i1)') jfaces
221!
222         if(vnorm.gt.0.d0) then
223            value=um*pressure*vnorm*eta*f*area/areaslav
224            call loadadd(nelems,label,value,nelemload,sideload,xloadact,
225     &                nload,nload_,iamload,iamplitude,nam,isector,
226     &                idefload)
227         elseif(vnorm.lt.-0.5d0) then
228!
229!           calculate the differential velocity
230!
231!           determining the slave velocity
232!
233            xi=pslavsurf(1,igauss)
234            et=pslavsurf(2,igauss)
235            iflag=1
236            if(nopes.eq.8) then
237               call shape8q(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
238            elseif(nopes.eq.4) then
239               call shape4q(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
240            elseif(nopes.eq.6) then
241               call shape6tri(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,
242     &                        iflag)
243            else
244               call shape3tri(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,
245     &                        iflag)
246            endif
247!
248            do k=1,3
249               vels(k)=0.d0
250               do j=1,nopes
251                  vels(k)=vels(k)+shp2s(4,j)*
252     &                            veold(k,konl(nopem+j))
253               enddo
254            enddo
255!
256!           determining the master velocity
257!
258            xi=pmastsurf(1,igauss)
259            et=pmastsurf(2,igauss)
260            iflag=1
261            if(nopem.eq.8) then
262               call shape8q(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
263            elseif(nopem.eq.4) then
264               call shape4q(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
265            elseif(nopem.eq.6) then
266               call shape6tri(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
267            else
268               call shape3tri(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
269            endif
270!
271            do k=1,3
272               velm(k)=0.d0
273               do j=1,nopem
274                  velm(k)=velm(k)+shp2m(4,j)*
275     &                            veold(k,konl(j))
276               enddo
277            enddo
278!
279            vnorm=dsqrt((vels(1)-velm(1))**2+
280     &                  (vels(2)-velm(2))**2+
281     &                  (vels(3)-velm(3))**2)
282            value=um*pressure*vnorm*eta*f*area/areaslav
283            call loadadd(nelems,label,value,nelemload,sideload,xloadact,
284     &                nload,nload_,iamload,iamplitude,nam,isector,
285     &                idefload)
286         endif
287!
288!        heat flux into the master face
289!
290         if(nopem.eq.8) then
291            mint2d=9
292         elseif(nopem.eq.6) then
293            mint2d=3
294         elseif(nopem.eq.4) then
295            mint2d=4
296         else
297            mint2d=1
298         endif
299!
300!        calculating the area of the slave face
301!
302         areamast=0.d0
303!
304         do j=1,mint2d
305            if(nopem.eq.8) then
306               xi=gauss2d3(1,j)
307               et=gauss2d3(2,j)
308               weight=weight2d3(j)
309            elseif(nopem.eq.6) then
310               xi=gauss2d5(1,j)
311               et=gauss2d5(2,j)
312               weight=weight2d5(j)
313            elseif(nopem.eq.4) then
314               xi=gauss2d2(1,j)
315               et=gauss2d2(2,j)
316               weight=weight2d2(j)
317            else
318               xi=gauss2d4(1,j)
319               et=gauss2d4(2,j)
320               weight=weight2d4(j)
321            endif
322!
323            iflag=2
324            if(nopem.eq.8) then
325               call shape8q(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
326            elseif(nopem.eq.4) then
327               call shape4q(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
328            elseif(nopem.eq.6) then
329               call shape6tri(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
330            else
331               call shape3tri(xi,et,pl,xsj2m,xs2m,shp2m,iflag)
332            endif
333!
334            areamast=areamast+dsqrt(xsj2m(1)**2+
335     &                              xsj2m(2)**2+
336     &                              xsj2m(3)**2)
337         enddo
338!
339         label(1:20)='S                   '
340         write(label(2:2),'(i1)') jfacem
341!
342!        at this point vnorm was either given by the user or
343!        calculated for the slave surface (differential velocity)
344!
345         value=um*pressure*vnorm*eta*(1.d0-f)*area/areamast
346         call loadadd(nelemm,label,value,nelemload,sideload,xloadact,
347     &                nload,nload_,iamload,iamplitude,nam,isector,
348     &                idefload)
349      enddo
350!
351      return
352      end
353
354