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