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 applybounfem(nodeboun,ndirboun,nboun,xbounact, 20 & ithermal,nk,vold,isolidsurf, 21 & nsolidsurf,xsolidsurf,nfreestream,ifreestream,iturbulent, 22 & vcon,shcon,nshcon,rhcon,nrhcon,ntmat_,physcon,v, 23 & compressible,nmpc,nodempc,ipompc,coefmpc,inomat, 24 & mi,ilboun,ilmpc,labmpc,coefmodmpc, 25 & ifreesurface,ierr,dgravity,depth,iexplicit) 26! 27! 1) applies temperature and velocity SPC's for 28! incompressible fluids (liquids) 29! materials (the pressure BC's have been applied in applybounp.f) 30! 2) applies temperature, velocity and pressure SPC's for 31! compressible fluids (gases) 32! 3) applies MPC's for the conservative variables 33! 34 implicit none 35! 36 character*20 labmpc(*) 37! 38 integer iturbulent,compressible,ifreesurface,ierr, 39 & nrhcon(*),mi(*),ntmat_,nodeboun(*), 40 & isolidsurf(*),j,ilboun(*),ilmpc(*), 41 & ndirboun(*),nshcon(*),nk,i,nboun,node,imat,ithermal(*), 42 & iexplicit, 43 & nsolidsurf,ifreestream(*),nfreestream, 44 & index,nodei,nmpc,nodempc(3,*),ipompc(*), 45 & ist,ndir,ndiri,inomat(*),nref 46! 47 real*8 rhcon(0:1,ntmat_,*),rho,vold(0:mi(2),*),xbounact(*), 48 & shcon(0:3,ntmat_,*),xnorm,coefmodmpc(*),dgravity, 49 & temp,xsolidsurf(*),sum,depth(*), 50 & vcon(nk,0:mi(2)),physcon(*), 51 & coefmpc(*),residu,correction,xkin,xtu,sumk,sumt, 52 & dvi,v(nk,0:mi(2)) 53! 54 nref=0 55! 56! SPC's: temperature, velocity and pressure (latter only for 57! compressible fluids) 58! 59 do j=1,nboun 60! 61! monotonically increasing DOF-order 62! 63 i=ilboun(j) 64! 65! pressure boundary conditions for incompressible materials 66! have already been treated 67! 68 ndir=ndirboun(i) 69 if((iexplicit.eq.0).and.(ndir.gt.3)) cycle 70! 71! check whether fluid node 72! 73 node=nodeboun(i) 74! 75! calculating the conservative variables for the previously 76! treated node, if any 77! 78 if(nref.le.nk) then 79 if(node.ne.nref) then 80 if(nref.ne.0) then 81 call phys2con(inomat,nref,vold,ntmat_,shcon,nshcon, 82 & physcon,compressible,vcon,rhcon,nrhcon,ithermal,mi, 83 & ifreesurface,ierr,dgravity,depth,nk) 84 endif 85 nref=node 86 endif 87 endif 88! 89! calculating the physical variables for the node at stake 90! 91 vold(ndir,node)=xbounact(i) 92 if(node.le.nk) v(node,ndir)=0.d0 93 enddo 94! 95! treating the remaining node 96! 97 if((nref.ne.0).and.(nref.le.nk)) then 98 call phys2con(inomat,nref,vold,ntmat_,shcon,nshcon, 99 & physcon,compressible,vcon,rhcon,nrhcon,ithermal,mi, 100 & ifreesurface,ierr,dgravity,depth,nk) 101 endif 102! 103! MPC's: temperature, velocity and pressure 104! 105 if(compressible.eq.0) then 106! 107! incompressible fluids 108! 109 nref=0 110! 111 do j=1,nmpc 112! 113! monotonically increasing DOF-order 114! 115 i=ilmpc(j) 116 ist=ipompc(i) 117! 118! pressure multiple point constraints for incompressible materials 119! have already been treated 120! 121 ndir=nodempc(2,ist) 122! 123! check whether fluid node 124! 125 node=nodempc(1,ist) 126! 127! calculating the value of the dependent DOF of the MPC 128! 129 index=nodempc(3,ist) 130 if(index.eq.0) cycle 131 sum=0.d0 132 do 133 nodei=nodempc(1,index) 134 ndiri=nodempc(2,index) 135 sum=sum+coefmpc(index)*vold(ndiri,nodei) 136 index=nodempc(3,index) 137 if(index.eq.0) exit 138 enddo 139! 140! calculating the conservative variables for the previously 141! treated node, if any 142! 143 if(nref.le.nk) then 144 if(node.ne.nref) then 145 if(nref.ne.0) then 146 call phys2con(inomat,nref,vold,ntmat_,shcon,nshcon, 147 & physcon,compressible,vcon,rhcon,nrhcon,ithermal,mi, 148 & ifreesurface,ierr,dgravity,depth,nk) 149 endif 150 nref=node 151 endif 152 endif 153! 154! calculating the physical variables for the node at stake 155! 156 vold(ndir,node)=-sum/coefmpc(ist) 157 if(node.le.nk) v(node,ndir)=0.d0 158 enddo 159! 160! treating the remaining node 161! 162 if((nref.ne.0).and.(nref.le.nk)) then 163 call phys2con(inomat,nref,vold,ntmat_,shcon,nshcon, 164 & physcon,compressible,vcon,rhcon,nrhcon,ithermal,mi, 165 & ifreesurface,ierr,dgravity,depth,nk) 166 endif 167 else 168! 169! compressible fluids 170! 171! MPC's are treated by distributing the residual proportional to 172! the coefficients 173! 174! Right now it is assumed that the MPC's are independent of each 175! other, i.e. degrees of freedom used in one MPC are not used in 176! any other MPC 177! 178 nref=0 179! 180 do j=1,nmpc 181 i=ilmpc(j) 182 index=ipompc(i) 183! 184! pressure multiple point constraints for incompressible materials 185! have already been treated 186! 187 ndir=nodempc(2,index) 188! 189! check whether fluid node 190! 191 node=nodempc(1,index) 192! 193! calculating the value of the dependent DOF of the MPC 194! 195 residu=coefmpc(index)*vold(ndir,node) 196 xnorm=1.d0 197 if(index.eq.0) cycle 198 do 199 index=nodempc(3,index) 200 if(index.eq.0) exit 201 nodei=nodempc(1,index) 202 ndiri=nodempc(2,index) 203 residu=residu+coefmpc(index)*vold(ndiri,nodei) 204 enddo 205! 206! correcting all terms of the MPC 207! 208 index=ipompc(i) 209 do 210 nodei=nodempc(1,index) 211 ndiri=nodempc(2,index) 212! 213 correction=-residu*coefmodmpc(index) 214 vold(ndiri,nodei)=vold(ndiri,nodei)+correction 215 if(nodei.le.nk) then 216 call phys2con(inomat,nodei,vold,ntmat_,shcon,nshcon, 217 & physcon,compressible,vcon,rhcon,nrhcon,ithermal,mi, 218 & ifreesurface,ierr,dgravity,depth,nk) 219 v(nodei,ndiri)=0.d0 220 endif 221 index=nodempc(3,index) 222 if(index.eq.0) exit 223 enddo 224 enddo 225! 226 endif 227! 228 if(iturbulent.ne.0) then 229! 230! freestream conditions for the iturbulent variables 231! 232 xtu=10.d0*physcon(5)/physcon(8) 233 xkin=10.d0**(-3.5d0)*xtu 234 do j=1,nfreestream 235 node=ifreestream(j) 236 imat=inomat(node) 237 if(imat.eq.0) cycle 238 temp=vold(0,node) 239 call materialdata_dvifem(imat,ntmat_,temp,shcon,nshcon,dvi) 240! 241! density 242! 243 rho=vcon(node,4) 244! 245 vcon(node,5)=xkin*dvi 246 vcon(node,6)=xtu*rho 247! 248 vold(5,node)=vcon(node,5)/rho 249 vold(6,node)=xtu 250! 251 v(node,5)=0.d0 252 v(node,6)=0.d0 253 enddo 254! 255! solid boundary conditions for the turbulent variables 256! 257 do j=1,nsolidsurf 258! 259! turbulent kinetic energy is applied at the wall 260! 261 node=isolidsurf(j) 262 imat=inomat(node) 263 if(imat.eq.0) cycle 264 temp=vold(0,node) 265 call materialdata_dvifem(imat,ntmat_,temp,shcon,nshcon,dvi) 266 rho=vcon(node,4) 267 vcon(node,5)=0.d0 268 vcon(node,6)=800.d0*dvi/(xsolidsurf(j)**2) 269! 270 vold(5,node)=0.d0 271 vold(6,node)=vcon(node,6)/rho 272! 273 v(node,5)=0.d0 274 v(node,6)=0.d0 275 enddo 276! 277! taking fluid pressure MPC's into account: it is assumed 278! that cyclic fluid pressure MPC's also apply to the iturbulent 279! conservative variables 280! 281 do i=1,nmpc 282 if(labmpc(i)(1:6).ne.'CYCLIC') cycle 283 ist=ipompc(i) 284 ndir=nodempc(2,ist) 285 if(ndir.ne.4) cycle 286 node=nodempc(1,ist) 287! 288! check whether fluid MPC 289! 290 imat=inomat(node) 291 if(imat.eq.0) cycle 292! 293 index=nodempc(3,ist) 294 if(index.eq.0) cycle 295! 296 sumk=0.d0 297 sumt=0.d0 298! 299 do 300 nodei=nodempc(1,index) 301 sumk=sumk+coefmpc(index)*vcon(nodei,5) 302 sumt=sumt+coefmpc(index)*vcon(nodei,6) 303 index=nodempc(3,index) 304 if(index.eq.0) exit 305 enddo 306 vcon(node,5)=-sumk/coefmpc(ist) 307 vcon(node,6)=-sumt/coefmpc(ist) 308 enddo 309 endif 310! 311 return 312 end 313 314