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