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 advecstiff(nope,voldl,ithermal,xl,nelemload,nelemadvec, 20 & nload,lakon,xload,istep,time,ttime,dtime,sideload,vold,mi, 21 & xloadold,reltime,nmethod,s,iinc,iponoel,inoel,ielprop,prop, 22 & ielmat,shcon,nshcon,rhcon,nrhcon,ntmat_,ipkon,kon,cocon,ncocon, 23 & ipobody,xbody,ibody) 24! 25! calculates the stiffness of an advective element. 26! An advective element consists of a face with a forced convection 27! film condition and a network node 28! 29 implicit none 30! 31 character*8 lakonl,lakon(*) 32 character*20 sideload(*),sideloadl 33! 34 integer nope,i,ithermal(*),j,nelemload(2,*),nelemadvec,nload,id, 35 & nelem,ig,mint2d,iflag,istep,jltyp,nfield,mi(*),nmethod,k,iinc, 36 & node,nopes,iponoel(*),inoel(2,*),ielprop(*),ielmat(mi(3),*), 37 & ipkon(*), 38 & nshcon(*),nrhcon(*),ntmat_,kon(*),ncocon(2,*),ipobody(2,*), 39 & ibody(3,*) 40! 41 real*8 tl2(9),voldl(0:mi(2),9),xl(3,9),sinktemp,xi,et,weight, 42 & xl2(3,8),xsj2(3),shp2(7,9),coords(3),xs2(3,7),dxsj2,areaj, 43 & temp,xload(2,*),timeend(2),time,ttime,dtime,field,reltime, 44 & vold(0:mi(2),*),xloadold(2,*),s(60,60),sref,sref2,prop(*), 45 & shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*),cocon(0:6,ntmat_,*), 46 & xbody(7,*),heatnod,heatfac 47! 48! 49! 50 include "gauss.f" 51! 52 iflag=2 53! 54 timeend(1)=time 55 timeend(2)=ttime+time 56! 57! number of nodes in the advective face 58! 59 nopes=nope-1 60! 61! temperature and displacements in the element's nodes 62! 63 do i=1,nope 64 tl2(i)=voldl(0,i) 65 enddo 66 do i=1,nopes 67 if(ithermal(2).eq.2) then 68 do j=1,3 69 xl2(j,i)=xl(j,i) 70 enddo 71 else 72 do j=1,3 73 xl2(j,i)=xl(j,i)+voldl(j,i) 74 enddo 75 endif 76 enddo 77! 78 call nident2(nelemload,nelemadvec,nload,id) 79! 80! the second entry in nelemload points to the original 81! film loading 82! 83 id=nelemload(2,id) 84! 85! number of the original element 86! 87 nelem=nelemload(1,id) 88 lakonl=lakon(nelem) 89 read(sideload(id)(2:2),'(i1)') ig 90! 91! number of integration points 92! 93 if(lakonl(4:5).eq.'8R') then 94 mint2d=1 95 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then 96 if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'E')) then 97 mint2d=2 98 else 99 mint2d=4 100 endif 101 elseif(lakonl(4:4).eq.'2') then 102 mint2d=9 103 elseif(lakonl(4:5).eq.'10') then 104 mint2d=3 105 elseif(lakonl(4:4).eq.'4') then 106 mint2d=1 107 elseif(lakonl(4:5).eq.'15') then 108 if(ig.le.2) then 109 mint2d=3 110 else 111 mint2d=4 112 endif 113 elseif(lakonl(4:4).eq.'6') then 114 mint2d=1 115 endif 116! 117 do i=1,mint2d 118! 119! copying the sink temperature to ensure the same 120! value in each integration point (sinktemp can be 121! changed in subroutine film: requirement from the 122! thermal people) 123! 124 sinktemp=tl2(nope) 125! 126 if((lakonl(4:5).eq.'8R').or. 127 & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then 128 xi=gauss2d1(1,i) 129 et=gauss2d1(2,i) 130 weight=weight2d1(i) 131 elseif((lakonl(4:4).eq.'8').or. 132 & (lakonl(4:6).eq.'20R').or. 133 & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then 134 xi=gauss2d2(1,i) 135 et=gauss2d2(2,i) 136 weight=weight2d2(i) 137 elseif(lakonl(4:4).eq.'2') then 138 xi=gauss2d3(1,i) 139 et=gauss2d3(2,i) 140 weight=weight2d3(i) 141 elseif((lakonl(4:5).eq.'10').or. 142 & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then 143 xi=gauss2d5(1,i) 144 et=gauss2d5(2,i) 145 weight=weight2d5(i) 146 elseif((lakonl(4:4).eq.'4').or. 147 & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then 148 xi=gauss2d4(1,i) 149 et=gauss2d4(2,i) 150 weight=weight2d4(i) 151 endif 152! 153 if(nopes.eq.8) then 154 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 155 elseif(nopes.eq.4) then 156 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 157 elseif(nopes.eq.6) then 158 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 159 else 160 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 161 endif 162! 163 dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 164 & xsj2(3)*xsj2(3)) 165 areaj=dxsj2*weight 166! 167 temp=0.d0 168 do j=1,nopes 169 temp=temp+tl2(j)*shp2(4,j) 170 enddo 171! 172! for nonuniform load: determine the coordinates of the 173! point (transferred into the user subroutine) 174! 175 if((sideload(id)(3:4).eq.'NU').or. 176 & (sideload(id)(5:6).eq.'NU')) then 177 do k=1,3 178 coords(k)=0.d0 179 do j=1,nopes 180 coords(k)=coords(k)+xl2(k,j)*shp2(4,j) 181 enddo 182 enddo 183 read(sideload(id)(2:2),'(i1)') jltyp 184 jltyp=jltyp+10 185 node=nelemload(2,id) 186 sideloadl=sideload(id) 187 sideloadl(1:1)='F' 188 call film(xload(1,id),sinktemp,temp,istep, 189 & iinc,timeend,nelem,i,coords,jltyp,field,nfield, 190 & sideloadl,node,areaj,vold,mi, 191 & ipkon,kon,lakon,iponoel,inoel,ielprop,prop,ielmat, 192 & shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon, 193 & ipobody,xbody,ibody,heatnod,heatfac) 194c if(nmethod.eq.1) xload(1,id)=xloadold(1,id)+ 195c & (xload(1,id)-xloadold(1,id))*reltime 196 endif 197! 198 sref=xload(1,id)*areaj 199 shp2(4,nope)=-1.d0 200! 201 do j=1,nope 202 sref2=sref*shp2(4,j) 203 s(nope,j)=s(nope,j)-sref2 204 do k=1,nopes 205 s(k,j)=s(k,j)+sref2*shp2(4,k) 206 enddo 207 enddo 208 enddo 209! 210! for some axisymmetric and plane strain elements only half the 211! number of integration points was used 212! 213 if(((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')).and. 214 & ((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'E'))) then 215 do i=1,nope 216 do j=1,nope 217 s(i,j)=2.d0*s(i,j) 218 enddo 219 enddo 220 endif 221! 222 return 223 end 224 225