-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathm_bubbles.f90
executable file
·254 lines (197 loc) · 7.51 KB
/
m_bubbles.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Module computing bubble dynamics
!
! Last update: May 19, 2009
! Author: Keita Ando
! Department of Mechanical Engineering
! Division of Engineering and Applied Science
! California Institute of Technology, Pasadena CA 91125
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
MODULE m_bubbles
USE m_globalvar
IMPLICIT NONE
! bubble wall properties
REAL(KIND(0.D0)) :: chi_vw
REAL(KIND(0.D0)) :: k_bw
REAL(KIND(0.D0)) :: rho_bw
CONTAINS
!========================================================================
SUBROUTINE s_bwproperty( pb )
REAL(KIND(0.D0)), INTENT(IN) :: pb
REAL(KIND(0.D0)) :: x_vw
IF ( vapor=='n' ) THEN
! mass fraction of vapor
chi_vw = 0.D0
! mole fraction of vapor & thermal conductivity of gas mixture
x_vw = 0.D0
k_bw = k_n(iR0)
! gas mixture density
rho_bw = pb/( R_n*Tw )
ELSE
! mass fraction of vapor
chi_vw = 1.D0/( 1.D0+R_v/R_n*(pb/pv-1.D0) )
! mole fraction of vapor & thermal conductivity of gas mixture
x_vw = M_n*chi_vw/( M_v+(M_n-M_v)*chi_vw )
k_bw = x_vw*k_v(iR0)/( x_vw+(1.D0-x_vw)*phi_vn ) &
+ ( 1.D0-x_vw )*k_n(iR0)/( x_vw*phi_nv+1.D0-x_vw )
! gas mixture density
rho_bw = pv/( chi_vw*R_v*Tw )
END IF
END SUBROUTINE s_bwproperty
!========================================================================
FUNCTION f_vflux( R,V,mass_v )
REAL(KIND(0.D0)), INTENT(IN) :: R
REAL(KIND(0.D0)), INTENT(IN) :: V
REAL(KIND(0.D0)), INTENT(IN) :: mass_v
REAL(KIND(0.D0)) :: chi_bar
REAL(KIND(0.D0)) :: grad_chi
REAL(KIND(0.D0)) :: f_vflux
IF ( vapor=='y' ) THEN
IF ( thermal=='transfer' ) THEN
! constant transfer model
chi_bar = mass_v/( mass_v+mass_n0(iR0) )
grad_chi = -Re_trans_c(iR0)*( chi_bar-chi_vw )
f_vflux = rho_bw*grad_chi/Pe_c/( 1.D0-chi_vw )/R
ELSE
! irrerevant procedure (polytropic)
f_vflux = pv*V/( R_v*Tw )
END IF
ELSE IF ( vapor=='n' ) THEN
f_vflux = 0.D0
END IF
END FUNCTION f_vflux
!========================================================================
FUNCTION f_bpres_dot( vflux,R,V,pb,mass_v )
REAL(KIND(0.D0)), INTENT(IN) :: vflux
REAL(KIND(0.D0)), INTENT(IN) :: R
REAL(KIND(0.D0)), INTENT(IN) :: V
REAL(KIND(0.D0)), INTENT(IN) :: pb
REAL(KIND(0.D0)), INTENT(IN) :: mass_v
REAL(KIND(0.D0)) :: T_bar
REAL(KIND(0.D0)) :: grad_T
REAL(KIND(0.D0)) :: tmp1, tmp2
REAL(KIND(0.D0)) :: f_bpres_dot
IF ( thermal=='transfer'.AND.model=='Preston' ) THEN
! constant transfer model
T_bar = Tw*( pb/pb0(iR0) )*( R/R0(iR0) )**3 &
* ( mass_n0(iR0)+mass_v0(iR0) )/( mass_n0(iR0)+mass_v )
grad_T = -Re_trans_T(iR0)*( T_bar-Tw )
f_bpres_dot = 3.D0*gamma_b*( -V*pb+vflux*R_v*Tw &
+ pb0(iR0)*k_bw*grad_T/Pe_T(iR0)/R )/R
ELSE IF ( thermal=='transfer'.AND.model=='Sugiyama' ) THEN
! only for gas bubble
T_bar = Tw*( pb/pb0(iR0) )*( R/R0(iR0) )**3
tmp1 = -V*pb - pb0(iR0)*k_bw/Pe_T(iR0)/R &
* ( Re_trans_T(iR0)*(T_bar-Tw) &
+ Im_trans_T(iR0)*T_bar/omegaN(iR0)*3.D0*V/R )
tmp2 = R/3.D0/gamma_b &
+ pb0(iR0)*k_bw*Im_trans_T(iR0)*T_bar &
/ Pe_T(iR0)/R/omegaN(iR0)/pb
f_bpres_dot = tmp1/tmp2
ELSE
f_bpres_dot = -3.D0*gamma_b*V/R*( pb-pv )
END IF
END FUNCTION f_bpres_dot
!========================================================================
FUNCTION f_bwpres_1( R,V,pb )
REAL(KIND(0.D0)), INTENT(IN) :: R
REAL(KIND(0.D0)), INTENT(IN) :: V
REAL(KIND(0.D0)), INTENT(IN) :: pb
REAL(KIND(0.D0)) :: tmp
REAL(KIND(0.D0)) :: f_bwpres_1
IF ( polytropic=='y' ) THEN
tmp = ( R0(iR0)/R )**( 3.D0*gamma_b )
tmp = ( Ca+2.D0/We/R0(iR0) )*tmp - Ca + 1.D0
ELSE IF ( polytropic=='n' ) THEN
tmp = pb
END IF
f_bwpres_1 = tmp - 4.D0*Re_inv*V/R - 2.D0/( We*R )
END FUNCTION f_bwpres_1
!========================================================================
FUNCTION f_bwpres( R,V,pb )
REAL(KIND(0.D0)), INTENT(IN) :: R
REAL(KIND(0.D0)), INTENT(IN) :: V
REAL(KIND(0.D0)), INTENT(IN) :: pb
REAL(KIND(0.D0)) :: tmp
REAL(KIND(0.D0)) :: f_bwpres
IF ( polytropic=='y' ) THEN
tmp = ( R0(iR0)/R )**( 3.D0*gamma_b )
tmp = ( Ca+2.D0/We/R0(iR0) )*tmp - Ca
ELSE IF ( polytropic=='n' ) THEN
tmp = pb - 1.D0
END IF
f_bwpres = tmp - 4.D0*Re_inv*V/R - 2.D0/( R*We )
END FUNCTION f_bwpres
!========================================================================
FUNCTION f_enthal( Cpinf,Cpbw )
REAL(KIND(0.D0)), INTENT(IN) :: Cpinf
REAL(KIND(0.D0)), INTENT(IN) :: Cpbw
REAL(KIND(0.D0)) :: tmp1, tmp2, tmp3
REAL(KIND(0.D0)) :: f_enthal
tmp1 = ( n_tait-1.D0 )/n_tait
tmp2 = ( Cpbw/(1.D0+B_tait)+1.D0 )**tmp1
tmp3 = ( Cpinf/(1.D0+B_tait)+1.D0 )**tmp1
f_enthal = ( tmp2-tmp3 )*n_tait*( 1.D0+B_tait )/( n_tait-1.D0 )
END FUNCTION f_enthal
!========================================================================
FUNCTION f_enthal_dot( R,V,Cpbw,Cpinf,pb_dot,Cpinf_dot )
REAL(KIND(0.D0)), INTENT(IN) :: R
REAL(KIND(0.D0)), INTENT(IN) :: V
REAL(KIND(0.D0)), INTENT(IN) :: Cpbw
REAL(KIND(0.D0)), INTENT(IN) :: Cpinf
REAL(KIND(0.D0)), INTENT(IN) :: pb_dot
REAL(KIND(0.D0)), INTENT(IN) :: Cpinf_dot
REAL(KIND(0.D0)) :: tmp1, tmp2
REAL(KIND(0.D0)) :: f_enthal_dot
IF ( polytropic=='y' ) THEN
tmp1 = ( R0(iR0)/R )**( 3.D0*gamma_b )
tmp1 = -3.D0*gamma_b*( Ca+2.D0/We/R0(iR0) )*tmp1*V/R
ELSE IF ( polytropic=='n' ) THEN
tmp1 = pb_dot
END IF
tmp2 = ( 2.D0/We+4.D0*Re_inv*V )*V/R**2
f_enthal_dot = &
( Cpbw/(1.D0+B_tait)+1.D0 )**( -1.D0/n_tait )*( tmp1+tmp2 ) &
- ( Cpinf/(1.D0+B_tait)+1.D0 )**( -1.D0/n_tait )*Cpinf_dot
END FUNCTION f_enthal_dot
!========================================================================
FUNCTION f_sound( Cpinf,H )
REAL(KIND(0.D0)), INTENT(IN) :: Cpinf
REAL(KIND(0.D0)), INTENT(IN) :: H
REAL(KIND(0.D0)) :: tmp
REAL(KIND(0.D0)) :: f_sound
tmp = ( Cpinf/(1.D0+B_tait)+1.D0 )**( (n_tait-1.D0)/n_tait )
tmp = n_tait*( 1.D0+B_tait )*tmp
f_sound = DSQRT( tmp+(n_tait-1.D0)*H )
END FUNCTION f_sound
!========================================================================
FUNCTION f_gilmore( R,V,pb,pb_dot,pinf,pinf_dot )
REAL(KIND(0.D0)), INTENT(IN) :: R
REAL(KIND(0.D0)), INTENT(IN) :: V
REAL(KIND(0.D0)), INTENT(IN) :: pb
REAL(KIND(0.D0)), INTENT(IN) :: pb_dot
REAL(KIND(0.D0)), INTENT(IN) :: pinf
REAL(KIND(0.D0)), INTENT(IN) :: pinf_dot
REAL(KIND(0.D0)) :: Cpinf
REAL(KIND(0.D0)) :: Cpbw
REAL(KIND(0.D0)) :: H
REAL(KIND(0.D0)) :: H_dot
REAL(KIND(0.D0)) :: A
REAL(KIND(0.D0)) :: tmp1, tmp2, tmp3
REAL(KIND(0.D0)) :: f_gilmore
Cpinf = pinf - pl0
Cpbw = f_bwpres( R,V,pb )
H = f_enthal( Cpinf,Cpbw )
H_dot = f_enthal_dot( R,V,Cpbw,Cpinf,pb_dot,pinf_dot )
A = f_sound( Cpinf,H )
tmp1 = V/A
tmp2 = 1.D0 + 4.D0*Re_inv/A/R*( Cpbw/(1.D0+B_tait)+1.D0 ) &
**( -1.D0/n_tait )
tmp3 = 1.5D0*V**2*( tmp1/3.D0-1.D0 ) + H*( 1.D0+tmp1 ) &
+ R*H_dot*( 1.D0-tmp1 )/A
f_gilmore = tmp3/( R*(1.D0-tmp1)*tmp2 )
END FUNCTION f_gilmore
!========================================================================
END MODULE m_bubbles