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 springforc_n2f_th(xl,vl,imat,elcon,nelcon, 20 & tnl,ncmat_,ntmat_,nope,kode,elconloc, 21 & plkcon,nplkcon,npmat_,mi,springarea,timeend,matname, 22 & node,noel,istep,iinc,iperturb) 23! 24! calculates the heat flux across a contact area 25! 26 implicit none 27! 28 character*80 matname(*),slname,msname 29! 30 integer i,j,imat,ncmat_,ntmat_,nope,nterms,iflag,mi(*), 31 & kode,niso,id,nplkcon(0:ntmat_,*),npmat_,nelcon(2,*), 32 & node,noel,istep,iinc,npred,iperturb(*) 33! 34 real*8 xl(3,10),ratio(9),t0l,t1l,al(3),vl(0:mi(2),10), 35 & pl(3,10),xn(3),dm,alpha,beta,tnl(10),pressure,dtemp, 36 & dist,conductance,eps,pi,springarea,timeend(2),ak(5), 37 & elcon(0:ncmat_,ntmat_,*),pproj(3),xsj2(3),xs2(3,7),val, 38 & shp2(7,9),xi,et,elconloc(21),plconloc(802),xk,d(2),flowm(2), 39 & xiso(200),yiso(200),plkcon(0:2*npmat_,ntmat_,*),temp(2), 40 & predef(2),coords(3),tmean 41! 42 iflag=2 43! 44! actual positions of the nodes belonging to the contact spring 45! (for geometrically linear static calculations the undeformed position 46! is taken) 47! 48 if(iperturb(2).eq.0) then 49 do i=1,nope 50 do j=1,3 51 pl(j,i)=xl(j,i) 52 enddo 53 enddo 54 else 55 do i=1,nope 56 do j=1,3 57 pl(j,i)=xl(j,i)+vl(j,i) 58 enddo 59 enddo 60 endif 61! 62 nterms=nope-1 63! 64! vector vr connects the dependent node with its projection 65! on the independent face 66! 67 do i=1,3 68 pproj(i)=pl(i,nope) 69 enddo 70 call attach_2d(pl,pproj,nterms,ratio,dist,xi,et) 71 do i=1,3 72 al(i)=pl(i,nope)-pproj(i) 73 enddo 74! 75! determining the jacobian vector on the surface 76! 77c if(nterms.eq.9) then 78c call shape9q(xi,et,pl,xsj2,xs2,shp2,iflag) 79 if(nterms.eq.8) then 80 call shape8q(xi,et,pl,xsj2,xs2,shp2,iflag) 81 elseif(nterms.eq.4) then 82 call shape4q(xi,et,pl,xsj2,xs2,shp2,iflag) 83 elseif(nterms.eq.6) then 84 call shape6tri(xi,et,pl,xsj2,xs2,shp2,iflag) 85c elseif(nterms.eq.7) then 86c call shape7tri(xi,et,pl,xsj2,xs2,shp2,iflag) 87 else 88 call shape3tri(xi,et,pl,xsj2,xs2,shp2,iflag) 89 endif 90! 91! normal on the surface 92! 93 dm=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+xsj2(3)*xsj2(3)) 94 do i=1,3 95 xn(i)=xsj2(i)/dm 96 enddo 97! 98! distance from surface along normal 99! 100 val=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3) 101! 102! representative area: usually the slave surface stored in 103! springarea; however, if no area was assigned because the 104! node does not belong to any element, the master surface 105! is used 106! 107 if(springarea.le.0.d0) then 108 if(nterms.eq.3) then 109 springarea=dm/2.d0 110 else 111 springarea=dm*4.d0 112 endif 113 endif 114! 115 if(elcon(1,1,imat).gt.0.d0) then 116! 117! exponential overclosure 118! 119 if(dabs(elcon(2,1,imat)).lt.1.d-30) then 120 pressure=0.d0 121 beta=1.d0 122 else 123! 124 alpha=elcon(2,1,imat) 125 beta=elcon(1,1,imat) 126 if(-beta*val.gt.23.d0-dlog(alpha)) then 127 beta=(dlog(alpha)-23.d0)/val 128 endif 129 pressure=dexp(-beta*val+dlog(alpha)) 130 endif 131 else 132! 133! linear overclosure 134! 135 pi=4.d0*datan(1.d0) 136 eps=elcon(1,1,imat)*pi/elcon(2,1,imat) 137 pressure=-elcon(2,1,imat)*val* 138 & (0.5d0+datan(-val/eps)/pi) 139 endif 140! 141! calculating the temperature difference across the contact 142! area and the mean temperature for the determination of the 143! conductance 144! 145 t1l=0.d0 146 do j=1,nterms 147 t1l=t1l+ratio(j)*vl(0,j) 148 enddo 149 dtemp=t1l-vl(0,nope) 150 tmean=(vl(0,nope)+t1l)/2.d0 151! 152! interpolating the material data according to temperature 153! 154 call materialdata_sp(elcon,nelcon,imat,ntmat_,i,tmean, 155 & elconloc,kode,plkcon,nplkcon,npmat_,plconloc,ncmat_) 156! 157! interpolating the material data according to pressure 158! 159 niso=int(plconloc(801)) 160! 161 if(niso.eq.0) then 162 d(1)=val 163 d(2)=pressure 164 temp(1)=vl(0,nope) 165 temp(2)=t1l 166 do j=1,3 167 coords(j)=xl(j,nope) 168 enddo 169 call gapcon(ak,d,flowm,temp,predef,timeend,matname(imat), 170 & slname,msname,coords,noel,node,npred,istep,iinc, 171 & springarea) 172 conductance=ak(1) 173 else 174 do i=1,niso 175 xiso(i)=plconloc(2*i-1) 176 yiso(i)=plconloc(2*i) 177 enddo 178 call ident(xiso,pressure,niso,id) 179 if(id.eq.0) then 180 xk=0.d0 181 conductance=yiso(1) 182 elseif(id.eq.niso) then 183 xk=0.d0 184 conductance=yiso(niso) 185 else 186 xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id)) 187 conductance=yiso(id)+xk*(pressure-xiso(id)) 188 endif 189 endif 190! 191! calculating the concentrated heat flow 192! 193 tnl(nope)=-springarea*conductance*dtemp 194 do j=1,nterms 195 tnl(j)=-ratio(j)*tnl(nope) 196 enddo 197! 198 return 199 end 200 201