File: N:\mfix\model\dqmom\source_population_eq.f
1
2 SUBROUTINE Source_population_eq(x,y,dydx)
3
4 use param1, only: zero, small_number
5 USE constant
6 USE fldvar
7 USE physprop
8 USE rdf
9 USE scalars
10
11 IMPLICIT NONE
12
13 double precision x,y(*),dydx(*)
14 double precision K_v
15 double precision m1,m2, dav,theta,c11
16 double precision epstotal
17 integer L,n,np,i,j,k,IJK
18
19 k_v=PI/6
20 epstotal=0.0
21 theta=0.0
22
23
24 DO I=1,Nscalar
25 A(I)= 0.0
26 omega(I)=y(I)/(k_v*(y(Nscalar+I)**3))
27 ENDDO
28
29 n=2*Nscalar
30 np=2*Nscalar
31
32 IJK=IJK_INDEX
33
34
35 DO i=1,n
36 DO j=1,n
37 matrix_a(i,j)=ZERO
38 matrix_b(i,j)=ZERO
39 matrix_c(i,j)=ZERO
40 inv_a(i,j)=ZERO
41 ENDDO
42 ENDDO
43
44 DO i=1,n
45 dydx(I)=ZERO
46 ENDDO
47
48 DO i=0,2*Nscalar-1
49 S_bar(I)=ZERO
50 ENDDO
51
52
53
54 DO I=1, Nscalar
55 epstotal=ROP_s(IJK,I)/RO_S(IJK,I)+epstotal
56 ENDDO
57
58 If(epstotal>=small_number) THEN
59 Do J=1,Nscalar
60 theta= theta+(ROP_s(IJk,J)/RO_S(IJK,J))*Theta_m(IJK,J)
61 ENDDO
62 theta= theta/epstotal
63 else
64 theta=0.0
65 endif
66
67
68 DO I=1,Nscalar
69 DO J=1,Nscalar
70
71 m1=PI*(Scalar(IJK,I)**3)*RO_S(IJK,I)/6.0
72 m2=PI*(Scalar(IJK,J)**3)*RO_S(IJK,J)/6.0
73 dav=(Scalar(IJk,I)+Scalar(IJk,J))/2.0
74 c11=((theta*(m1+m2)**2)/(4*PI*m1*m2))**(0.5)*(4.0/dav)
75 BETA_A(I,J)=aggregation_eff*PI*(dav**3)*G_0(IJK,I,J)*c11
76 A(I)=A(I)+PI*(dav**3)*G_0(IJK,I,J)*c11*omega(J)*breakage_eff
77
78 ENDDO
79 ENDDO
80
81
82
83
84 DO L=0,2*Nscalar-1
85 DO I=1,Nscalar
86 DO J=1,Nscalar
87 If(y(I)<=1.0E-6) THEN
88 IF(y(J)<=1.0E-6) THEN
89 S_bar(L)=S_bar(L)+0.0
90 ELSE
91 S_bar(L)=S_bar(L)+&
92 BETA_a(I,J)*(1.0/2.0)*omega(I)*omega(J)*&
93 ((y(Nscalar+J)**3)**(L/3.0))
94 ENDIF
95 ELSE
96 IF(y(J)<=1.0E-6)THEN
97 S_bar(L)=S_bar(L)+&
98 BETA_a(I,J)*(1.0/2.0)*omega(I)*omega(J)*&
99 ((y(Nscalar+I)**3)**(L/3.0))-&
100 BETA_a(I,J)*(y(Nscalar+I)**L)*omega(I)*omega(J)
101 ELSE
102 S_bar(L)=S_bar(L)+&
103 BETA_a(I,J)*(1.0/2.0)*omega(I)*omega(J)*&
104 (((y(Nscalar+I)**3)+(y(Nscalar+J)**3))**(L/3.0))-&
105 BETA_a(I,J)*(y(Nscalar+I)**L)*omega(I)*omega(J)
106 ENDIF
107 ENDIF
108 ENDDO
109 ENDDO
110
111 DO I=1,Nscalar
112 IF(y(I)<1.0E-6) THEN
113 S_bar(L)=S_bar(L)+0.0
114 ELSE
115 S_bar(L)=S_bar(L)+&
116 A(I)*(2**((3.0-L)/3.0))*(y(Nscalar+I)**L)*omega(I)-&
117 A(I)*(y(Nscalar+I)**L)*omega(I)
118 ENDIF
119 ENDDO
120 ENDDO
121
122
123
124
125 IF(NSCALAR==2) THEN
126 IF(abs(y(3)-y(4))<=1.0E-4) THEN
127 y(3)=y(3)*0.999
128 y(4)=y(4)*1.001
129 ENDIF
130 ENDIF
131
132 IF(NSCALAR==3) THEN
133 IF(abs(y(4)-y(5))<=1.0E-4) THEN
134 y(4)=y(4)*0.999
135 y(5)=y(4)*1.001
136 ENDIF
137 IF(abs(y(5)-y(6))<=1.0E-4) THEN
138 y(5)=y(5)*0.999
139 y(6)=y(5)*1.001
140 ENDIF
141 IF(abs(y(4)-y(6))<=1.0E-4) THEN
142 y(4)=y(4)*0.999
143 y(6)=y(4)*1.001
144 ENDIF
145 ENDIF
146
147 DO j=1,Nscalar
148 matrix_a(1,j)=1.0
149 matrix_a(2,j)=0.0
150 ENDDO
151
152 DO j=Nscalar+1,n
153 matrix_a(1,j)=0.0
154 matrix_a(2,j)=1.0
155 ENDDO
156
157 DO i=3,n
158 DO j=1,Nscalar
159 matrix_a(i,j)=(1.0-i+1.0)*(y(j+Nscalar)**((i-1)*1.0))
160 ENDDO
161
162 DO j=Nscalar+1,n
163 matrix_a(i,j)=(i-1.0)*(y(j)**((i-2)*1.0))
164 ENDDO
165 ENDDO
166
167
168
169
170 DO i=1,n
171 inv_a(i,i)=1.0
172 ENDDO
173
174 call gaussj(matrix_a,n,np,inv_a,n,n)
175
176 DO i=1,n
177 IF(i<=Nscalar) THEN
178 IF(y(i)>1.0E-3) matrix_c(i,i)=-2*(y(i+Nscalar)**3)
179 ELSE
180 IF(y(i-Nscalar)>1.0E-3) matrix_c(i,i)= 4*(y(i)**3)
181 ENDIF
182 ENDDO
183
184 DO i=1,Nscalar
185 IF(y(i)>1.0E-3) THEN
186 matrix_c(i,i+Nscalar)=3*(y(i+Nscalar)**2)
187 matrix_c(i+Nscalar,i)=-3*(y(i+Nscalar)**4)
188 ENDIF
189 ENDDO
190
191
192 DO i = 1, n
193 DO j = 1, n
194 DO k = 1, n
195 matrix_b(i,j) = matrix_b(i,j) + &
196 matrix_c(i,k)*inv_a(k,j)
197 ENDDO
198 ENDDO
199 ENDDO
200
201 DO i=1,Nscalar
202 DO j=1,n
203 dydx(i)=dydx(i)+k_v*matrix_b(i,j)*s_bar(j-1)
204 ENDDO
205 ENDDO
206
207 DO i=Nscalar+1,n
208 IF(y(i-Nscalar)>1.0E-6) THEN
209 dydx(i)=-dydx(i-Nscalar)*y(i)/y(i-Nscalar)
210 ELSE
211 dydx(i)=dydx(i)+0.0
212 ENDIF
213 DO j=1,n
214 IF(y(i-Nscalar)>1.0E-6) THEN
215 dydx(i)=dydx(i)+(k_v/y(i-Nscalar))*(matrix_b(i,j)*s_bar(j-1))
216 ELSE
217 dydx(i)=dydx(i)+0.0
218 ENDIF
219 ENDDO
220 ENDDO
221
222
223 return
224 END
225
226
227
228
229