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 extrapolatefluid(nk,ipofano,ifano,inum,vfa,v,ielfa, 20 & ithermal,imach,ikappa,xmach,xkappa,shcon,nshcon,ntmat_,ielmat, 21 & physcon,mi,iturb,xturb,gradtfa,gradvfa,gradpfa,gradkfa,gradofa, 22 & co,cofa,ifabou) 23! 24! extrapolates the field values at the center of the faces to 25! the nodes 26! 27 implicit none 28! 29 logical gradient 30! 31 integer nk,ipofano(*),ifano(2,*),inum(*),ielfa(4,*),i,l,indexf, 32 & iface,ithermal(*),imach,ikappa,imat,nshcon(*),ntmat_,mi(*), 33 & ielmat(mi(3),*),iturb,ifabou(*) 34! 35 real*8 vfa(0:7,*),v(0:mi(2),*),cp,r,xk,xmach(*),xkappa(*),t1l, 36 & shcon(0:3,ntmat_,*),physcon(*),xturb(2,*),gradtfa(3,*), 37 & gradvfa(3,3,*),gradpfa(3,*),gradkfa(3,*),gradofa(3,*), 38 & vfal(0:7),co(3,*),cofa(3,*) 39! 40c sum=0.d0 41! 42 do i=1,nk 43! 44! athermal calculations 45! 46 if(ithermal(1).eq.0) then 47 do l=1,4 48 v(l,i)=0.d0 49 enddo 50 inum(i)=0 51 indexf=ipofano(i) 52 do 53 if(indexf.eq.0) exit 54 iface=ifano(1,indexf) 55 if(ielfa(2,iface).gt.0) then 56 do l=1,3 57 vfal(l)=vfa(l,iface)+ 58 & gradvfa(l,1,iface)*(co(1,i)-cofa(1,iface))+ 59 & gradvfa(l,2,iface)*(co(2,i)-cofa(2,iface))+ 60 & gradvfa(l,3,iface)*(co(3,i)-cofa(3,iface)) 61 enddo 62 vfal(4)=vfa(4,iface)+ 63 & gradpfa(1,iface)*(co(1,i)-cofa(1,iface))+ 64 & gradpfa(2,iface)*(co(2,i)-cofa(2,iface))+ 65 & gradpfa(3,iface)*(co(3,i)-cofa(3,iface)) 66 else 67 do l=1,3 68 vfal(l)=vfa(l,iface) 69 enddo 70 vfal(4)=vfa(4,iface) 71 endif 72 do l=1,4 73 v(l,i)=v(l,i)+vfal(l) 74 enddo 75 inum(i)=inum(i)+1 76 indexf=ifano(2,indexf) 77 enddo 78 if(inum(i).gt.0) then 79 do l=1,4 80 v(l,i)=v(l,i)/inum(i) 81 enddo 82 endif 83 else 84! 85! 1) thermal incompressible calculations 86! 2) compressible calculations 87! 88 do l=0,4 89 v(l,i)=0.d0 90 enddo 91 inum(i)=0 92 indexf=ipofano(i) 93 do 94 if(indexf.eq.0) exit 95 iface=ifano(1,indexf) 96 if(ielfa(2,iface).gt.0) then 97! 98! internal node 99! 100 gradient=.true. 101 elseif(ielfa(2,iface).eq.0) then 102! 103! external node without boundary condition 104! 105 gradient=.false. 106 elseif(ielfa(2,iface).lt.0) then 107! 108! external node with boundary condition 109! 110 if(ifabou(-ielfa(2,iface)+5).lt.0) then 111! 112! sliding conditions 113! 114 gradient=.true. 115 else 116! 117! other conditions 118! 119 gradient=.false. 120 endif 121 endif 122 if(gradient) then 123 vfal(0)=vfa(0,iface)+ 124 & gradtfa(1,iface)*(co(1,i)-cofa(1,iface))+ 125 & gradtfa(2,iface)*(co(2,i)-cofa(2,iface))+ 126 & gradtfa(3,iface)*(co(3,i)-cofa(3,iface)) 127 do l=1,3 128 vfal(l)=vfa(l,iface)+ 129 & gradvfa(l,1,iface)*(co(1,i)-cofa(1,iface))+ 130 & gradvfa(l,2,iface)*(co(2,i)-cofa(2,iface))+ 131 & gradvfa(l,3,iface)*(co(3,i)-cofa(3,iface)) 132 enddo 133 vfal(4)=vfa(4,iface)+ 134 & gradpfa(1,iface)*(co(1,i)-cofa(1,iface))+ 135 & gradpfa(2,iface)*(co(2,i)-cofa(2,iface))+ 136 & gradpfa(3,iface)*(co(3,i)-cofa(3,iface)) 137 else 138 vfal(0)=vfa(0,iface) 139 do l=1,3 140 vfal(l)=vfa(l,iface) 141 enddo 142 vfal(4)=vfa(4,iface) 143 endif 144 do l=0,4 145 v(l,i)=v(l,i)+vfal(l) 146 enddo 147 if(imach.eq.1) then 148 t1l=vfal(0) 149 imat=ielmat(1,ielfa(1,iface)) 150 r=shcon(3,1,imat) 151 call materialdata_cp_sec(imat,ntmat_,t1l, 152 & shcon,nshcon,cp,physcon) 153 xk=cp/(cp-r) 154 xmach(i)=xmach(i)+dsqrt((vfal(1)**2+ 155 & vfal(2)**2+vfal(3)**2)/(xk*r*t1l)) 156 endif 157 if(ikappa.eq.1) then 158 xkappa(i)=xkappa(i)+xk 159 endif 160! 161 inum(i)=inum(i)+1 162 indexf=ifano(2,indexf) 163 enddo 164 if(inum(i).gt.0) then 165 do l=0,4 166 v(l,i)=v(l,i)/inum(i) 167 enddo 168 if(imach.eq.1) xmach(i)=xmach(i)/inum(i) 169 if(ikappa.eq.1) xkappa(i)=xkappa(i)/inum(i) 170 endif 171 endif 172! 173! turbulence output 174! 175 if(iturb.ne.0) then 176 do l=1,2 177 xturb(l,i)=0.d0 178 enddo 179 indexf=ipofano(i) 180 do 181 if(indexf.eq.0) exit 182 iface=ifano(1,indexf) 183 if(ielfa(2,iface).gt.0) then 184 vfal(6)=vfa(6,iface)+ 185 & gradkfa(1,iface)*(co(1,i)-cofa(1,iface))+ 186 & gradkfa(2,iface)*(co(2,i)-cofa(2,iface))+ 187 & gradkfa(3,iface)*(co(3,i)-cofa(3,iface)) 188c write(*,*) 'extrapolatefluid i1',i,vfal(6) 189c if(sum.lt.vfal(6)) sum=vfal(6) 190 vfal(7)=vfa(7,iface)+ 191 & gradofa(1,iface)*(co(1,i)-cofa(1,iface))+ 192 & gradofa(2,iface)*(co(2,i)-cofa(2,iface))+ 193 & gradofa(3,iface)*(co(3,i)-cofa(3,iface)) 194 else 195 vfal(6)=vfa(6,iface) 196c write(*,*) 'extrapolatefluid i2',i,vfal(6) 197c if(sum.lt.vfal(6)) sum=vfal(6) 198 vfal(7)=vfa(7,iface) 199 endif 200 do l=1,2 201 xturb(l,i)=xturb(l,i)+vfal(l+5) 202 enddo 203c inum(i)=inum(i)+1 204 indexf=ifano(2,indexf) 205 enddo 206 if(inum(i).gt.0) then 207 do l=1,2 208 xturb(l,i)=xturb(l,i)/inum(i) 209c write(*,*) 'extrapolatesum ',i,xturb(1,i) 210c if(sum.lt.abs(xturb(1,i))) sum=abs(xturb(1,i)) 211 enddo 212! 213! check for values with an exponent consisting of 3 digits 214! (are not correctly displayed in cgx) 215! 216 if(xturb(1,i).lt.1.e-30) xturb(1,i)=0.d0 217 endif 218 endif 219 enddo 220c write(*,*) 'extrapolatefluid sum',sum 221! 222 return 223 end 224