xref: /phasta/phSolver/incompressible/BCprofile.f (revision 32eed1412d9a3e748fe578866226c6021c4817b2)
1c-----------------------------------------------------------------------
2        subroutine BCprofileScale(vbc_prof,BC,yold)
3        use pvsQbi
4        include "common.h"
5        real*8 vbc_prof(nshg,3)
6        real*8 BC(nshg,ndofBC),yold(nshg,ndof)
7        real*8 offphase
8        integer factor
9
10        r_amp =rampmdot(1,1)
11        r_freq=rampmdot(2,1)
12!  Usual sinusoidal in time syn jet is given in the next line.  It assumes that you are NOT changing the time step from previous runs since it computes the
13!  current time in the sin function to be (lstep+1)* Dt where lstep is a running total of all of the runs up to now.  This will be "ok" under the following
14!  assumption:  lstep*dt_current = m * T where m is an integer and T is the period of the jet.  That is because this will evaluate to 2*m*Pi and the sin
15!  that is zero
16
17        if(abs(rampmdot(1,2)) .lt. 1.0e-5) then
18           r_time_factor = r_amp*sin(two*pi*r_freq*(lstep+1)*Delt(1)) ! BC set for next time step - all phases are computed
19        else
20           r_time_factor = rampmdot(1,2)*r_amp   ! read parameter scales Vmax
21        endif
22
23        icount = 0
24        do kk=1,nshg
25          if(ndsurf(kk).ge.1 .and. ndsurf(kk).le.12) then ! this means all the SJ BC (1->12) for the Beta (Boeing) model
26
27            factor = idnint(rampmdot(2,2))
28
29            if(factor == 0) then
30              offphase = 0.d0
31
32            elseif(factor == 1) then
33              offphase = 1.d0
34
35            elseif(factor == 2) then
36              !Start to count from tip. If factor == 2 -> 1, 3, 5, etc
37              if(mod(ndsurf(kk)-1,factor) == 0) then
38                offphase = 1.d0
39              else
40                offphase = 0.d0
41              endif
42
43            elseif(factor == 3) then
44              !Start to count from tip. If factor == 3 -> 2, 5, 8, etc
45              if(mod(ndsurf(kk)+1,factor) == 0) then
46                offphase = 1.d0
47              else
48                offphase = 0.d0
49              endif
50
51            elseif(factor < 0) then
52              !Only one jet blowing. If factor = -5, then jet 5 only
53              !blows
54              if (ndsurf(kk) == -factor) then
55                offphase = 1.d0
56              else
57                offphase = 0.d0
58              endif
59            endif
60
61            BC(kk,3)=r_time_factor*vbc_prof(kk,1)*offphase
62            BC(kk,4)=r_time_factor*vbc_prof(kk,2)*offphase
63            BC(kk,5)=r_time_factor*vbc_prof(kk,3)*offphase
64
65            icount = icount + 1
66
67          endif
68        enddo
69
70!        if(istep.eq.0 .and. icount.ne.0)
71!     &     write(*,*) 'BCprofile count',myrank,icount
72
73        return
74        end
75
76c--------------------------------------------------------------
77        subroutine BCprofileInit(vbc_prof,x)
78
79        use pvsQbi
80        include "common.h"
81        real*8 vbc_prof(nshg,3), x(numnp,nsd)
82        real*8 rcenter(3),rnorml(3)
83
84!        open(unit=789, file='bcprofile.dat',status='unknown')
85        do kk=1,nshg
86          if(ndsurf(kk).eq.1) then     ! SJ1
87            rnorml(1)= -0.03805041312991605
88            rnorml(2)= 0.99925652732035097
89            rnorml(3)= 0.00620956265097128
90            rcenter(1)= 2.31206845913690984
91            rcenter(2)= 0.39594943486752049
92            rcenter(3)= 0.51243427055630897
93
94          elseif(ndsurf(kk).eq.2) then     ! SJ2
95            rnorml(1)= -0.03805041312991605
96            rnorml(2)= 0.99925652732035097
97            rnorml(3)= 0.00620956265097128
98            rcenter(1)= 2.28612828683167457
99            rcenter(2)= 0.39519558316178066
100            rcenter(3)= 0.47479183928994600
101
102          elseif(ndsurf(kk).eq.3) then     ! SJ3
103            rnorml(1)= -0.03805041312991605
104            rnorml(2)= 0.99925652732035097
105            rnorml(3)= 0.00620956265097128
106            rcenter(1)= 2.26019111197956502
107            rcenter(2)= 0.39444184406245603
108            rcenter(3)= 0.43714965468096900
109
110          elseif(ndsurf(kk).eq.4) then     ! SJ4
111            rnorml(1)= -0.03806860790482020
112            rnorml(2)= 0.99925599805332099
113            rnorml(3)= 0.00618315830704482
114            rcenter(1)= 2.22406898326471003
115            rcenter(2)= 0.39338941696418994
116            rcenter(3)= 0.38472970195368150
117
118          elseif(ndsurf(kk).eq.5) then     ! SJ5
119            rnorml(1)= -0.03807026719147060
120            rnorml(2)= 0.99925594973518905
121            rnorml(3)= 0.00618075034236115
122            rcenter(1)= 2.19813166796875015
123            rcenter(2)= 0.39263411631338135
124            rcenter(3)= 0.34709106843832799
125
126          elseif(ndsurf(kk).eq.6) then     ! SJ6
127            rnorml(1)= -0.03807076874492700
128            rnorml(2)= 0.99925593512835698
129            rnorml(3)= 0.00618002248559851
130            rcenter(1)= 2.17219205684703009
131            rcenter(2)= 0.39187870083238074
132            rcenter(3)= 0.30944277157095951
133
134          elseif(ndsurf(kk).eq.7) then     ! SJ7
135            rnorml(1)= -0.03786596844535725
136            rnorml(2)= 0.99926183452064998
137            rnorml(3)= 0.00647722966396847
138            rcenter(1)= 2.13621068971991468
139            rcenter(2)= 0.39081967107127247
140            rcenter(3)= 0.25712280227877149
141
142          elseif(ndsurf(kk).eq.8) then     ! SJ8
143            rnorml(1)= -0.03805101631143114
144            rnorml(2)= 0.99925650979092495
145            rnorml(3)= 0.00620868731108326
146            rcenter(1)= 2.11020507860295004
147            rcenter(2)= 0.39007468380778049
148            rcenter(3)= 0.21948491547550850
149
150          elseif(ndsurf(kk).eq.9) then     ! SJ9
151            rnorml(1)= -0.03805101631143114
152            rnorml(2)= 0.99925650979092495
153            rnorml(3)= 0.00620868731108326
154            rcenter(1)= 2.08426231456637501
155            rcenter(2)= 0.38931079611337499
156            rcenter(3)= 0.18184066095664247
157
158          elseif(ndsurf(kk).eq.10) then    ! SJ10
159            rnorml(1)= -0.03427503641589490
160            rnorml(2)= 0.99934408637497807
161            rnorml(3)= 0.01168840904698280
162            rcenter(1)= 2.04819760342475998
163            rcenter(2)= 0.38841663756682199
164            rcenter(3)= 0.12950419445094899
165
166          elseif(ndsurf(kk).eq.11) then    ! SJ11
167            rnorml(1)= -0.03208677268269435
168            rnorml(2)= 0.99937455421864496
169            rnorml(3)= 0.01486403037855546
170            rcenter(1)= 2.02226290188708502
171            rcenter(2)= 0.38805554104084300
172            rcenter(3)= 0.09186760148390730
173
174          elseif(ndsurf(kk).eq.12) then     ! SJ12
175            rnorml(1)= -0.02932216953382100
176            rnorml(2)= 0.99939176776025096
177            rnorml(3)= 0.01887604055064800
178            rcenter(1)= 1.99632633336121512
179            rcenter(2)= 0.38789394943524053
180            rcenter(3)= 0.05422828360356270
181
182          endif
183
184           rdistfromcsq=(x(kk,1)-rcenter(1))*(x(kk,1)-rcenter(1))+
185     &                  (x(kk,2)-rcenter(2))*(x(kk,2)-rcenter(2))+
186     &                  (x(kk,3)-rcenter(3))*(x(kk,3)-rcenter(3))
187
188           radius = 0.0183896  ! 0.724 inches in meters
189
190c.............Factors below are negative for desired blowing direction
191           if(ndsurf(kk).ge.1 .and. ndsurf(kk).le.12) then  ! All jet blowing profiles
192!           if(ndsurf(kk).eq.9) then  ! Only jet 9  blowing for now
193
194              vbc_prof(kk,1)=-rnorml(1)*(1-rdistfromcsq/(radius*radius))
195              vbc_prof(kk,2)=-rnorml(2)*(1-rdistfromcsq/(radius*radius))
196              vbc_prof(kk,3)=-rnorml(3)*(1-rdistfromcsq/(radius*radius))
197
198!              write(789,987) kk,vbc_prof(kk,1),vbc_prof(kk,2),
199!     &                          vbc_prof(kk,3)
200
201           else
202              vbc_prof(kk,:)=zero
203           endif
204
205        enddo
206!        close(789)
207!987     format(i6,3(2x,e14.7))
208
209        return
210        end
211