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