xref: /phasta/phSolver/compressible/setBlowing_Duct2.f (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08)
1c================================================================
2c Set jet inlet BCs based on contraction mdot dynamically,
3c called by itrdrv.f
4c ===============================================================
5
6        subroutine setBlowing_Duct2(x,BC,y,iTurbWall,istp)
7
8        use blowingDuct ! njetinlet, jetinletf
9        use timedataC ! varts data
10        include "common.h"
11        include "mpif.h"
12        include "auxmpi.h"
13        integer istp
14        integer i, nn
15        real*8 xcoor,ycoor,zcoor
16        real*8 BC(nshg,ndofBC)
17        real*8 x(nshg,nsd)
18        real*8 y(nshg,ndof)
19        real*8 yVel, Temp, pres, mdot, rho, Area
20        real*8 xminm,xmaxm,zminm,zmaxm,xlm,zlm
21        real rx,rz
22        integer isfID
23        real*8 contraPres, contraTemp, contraXVel, contraRho,
24     &         contraArea, contraMdot
25        real*8 blowingPres,blowingTemp,blowingRho,blowingArea,
26     &         MdotRatio,blowingMdot,blowVel
27        integer contraProbeNo,blowingProbeNo
28        integer iTurbWall(nshg)
29
30c--- The following operations must be fulfilled on each proces because even if
31c--- there is no jet mouth on that process, it can be the host process that contains
32c--- all the varts data
33
34        if(myrank.eq.0)then ! only rank0 process can acess all varts data
35          contraProbeNo=1 ! in xyzts.dat the FIRST probe is at contraction inlet
36          blowingProbeNo=2
37          contraArea=1.17475**2 ! meter^2
38          blowingArea=0.00237064  !meter^2
39
40          contraPres=varts(contraProbeNo,1)
41          contraTemp=varts(contraProbeNo,5)
42          contraXVel=varts(contraProbeNo,2)
43          contraRho =contraPres/(287*contraTemp)
44          contraMdot=contraRho*contraArea*contraXVel !kg/s
45
46          blowingPres=varts(blowingProbeNo,1)
47          blowingTemp=varts(blowingProbeNo,5)
48          blowingRho = blowingPres/(287*blowingTemp)
49
50          MdotRatio =  (BlowingIniMdotDuct) +
51     &   real(istp-1)/real(nBlowingStepsDuct-1)
52     &   *(BlowingFnlMdotDuct-BlowingIniMdotDuct)
53
54          blowingMdot= (MdotRatio/100.0)*contraMdot   ! The absolute value of mdot
55          blowingVel = blowingMdot/(blowingArea*blowingRho)
56          write(*,*)'Blowing Velocity:', blowingVel
57        endif
58
59        call MPI_BARRIER(MPI_COMM_WORLD,ierr)
60        call MPI_BCAST(blowingVel,1,MPI_REAL8,master,
61     &                 MPI_COMM_WORLD,ierr)
62        call MPI_BARRIER(MPI_COMM_WORLD,ierr)
63
64
65c--- Only if current process contains jet inlet nodes
66        if(njetinlet .gt. 0)then
67          do i=1,njetinlet
68            nn=jetinletf(i)
69            if(iTurbWall(nn).eq.0)then
70c            xcoor=x(nn,1)
71c            ycoor=x(nn,2)
72c            zcoor=x(nn,3)
73c            rx=(xcoor-xminm)*(xmaxm-xcoor)/(xlm/2)**2
74c            rz=(zcoor-zminm)*(zmaxm-zcoor)/(zlm/2)**2
75c            rx=max(0.0,rx)
76c            rz=max(0.0,rz)
77c            pres = y(nn,4) ! the quantity updating each time step
78c           Temp = y(nn,5)
79c           rho = pres/(287.0*Temp)
80c              yVel = (blowingMdot/rho/Area)*rx*rz
81c............................
82              BC(nn,2) = 317          ! Temp
83              BC(nn,4) = blowingVel   ! set and scale y velocity
84              BC(nn,3) = 0
85              BC(nn,5) = 0
86              BC(nn,7) = 1.825e-5
87            endif
88          enddo
89
90        endif   ! only if current process contains jet inlet nodes
91
92        return
93        end
94
95