MFIX  2016-1
qmomk_collision_mod.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module: qmomk_collisions C
4 ! Purpose: Collision operators for the kinetic equation C
5 ! Contains the following subroutines and functions: C
6 ! collisions_istantaneous, collisions_bgk, C
7 ! compute_collision_time, collisions_boltzmann_one_specie, C
8 ! collisions_boltzmann_two_specie, C
9 ! solve_boltzmann_collisions_one_specie, C
10 ! solve_boltzmann_collisions_two_specie, radial_g0 C
11 ! C
12 ! Author: A. Passalacqua Date: C
13 ! Reviewer: Date: C
14 ! C
15 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
16 
17 #include "version.inc"
18 
19 MODULE qmomk_collision
20 
21 
22  USE qmomk_parameters
23  USE qmomk_quadrature
24 
25  IMPLICIT NONE
26 
27  PRIVATE
28 
29  PUBLIC :: collisions_istantaneous
30  PUBLIC :: collisions_bgk
31  PUBLIC :: compute_collision_time
32  !PUBLIC :: GLOBAL_COLLISION_TIME
33  PUBLIC :: solve_boltzmann_collisions_one_specie
34  PUBLIC :: solve_boltzmann_collisions_two_species
35 
36  PUBLIC :: radial_g0
37 
38 CONTAINS
39 
40 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
41 ! C
42 ! Instantaneous collisions C
43 ! C
44 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
45 
46  SUBROUTINE collisions_istantaneous(M, dp)
47 
48  IMPLICIT NONE
49 
50  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: M
51  DOUBLE PRECISION, INTENT(IN) :: dp
52 
53  DOUBLE PRECISION :: m0, m1, m2, m3, m11, m22, m33
54 
55  DOUBLE PRECISION :: up1, up2, up3
56  DOUBLE PRECISION :: sig1, sig2, sig3, seq
57 
58  DOUBLE PRECISION :: d11, d12, d13, d22, d23, d33
59  DOUBLE PRECISION :: d111, d112, d113, d122, d123, d133, d222, d223, d233, d333
60 
61  m0 = m(1)
62  m1 = m(2)
63  m2 = m(3)
64  m3 = m(4)
65  m11 = m(5)
66  m22 = m(8)
67  m33 = m(10)
68 
69  up1 = m1/m0
70  up2 = m2/m0
71  up3 = m3/m0
72  sig1 = max(0.d0, m11/m0 - up1**2)
73  sig2 = max(0.d0, m22/m0 - up2**2)
74  sig3 = max(0.d0, m33/m0 - up3**2)
75  seq = ( sig1 + sig2 + sig3 )/3.d0
76 
77  d11 = m0*(seq + up1**2)
78  d12 = m0*up1*up2
79  d13 = m0*up1*up3
80  d22 = m0*(seq + up2**2)
81  d23 = m0*up2*up3
82  d33 = m0*(seq + up3**2)
83  d111 = m1*(3.d0*seq + up1**2)
84  d112 = m2*(seq + up1**2)
85  d113 = m3*(seq + up1**2)
86  d122 = m1*(seq + up2**2)
87  d123 = m0*up1*up2*up3
88  d133 = m1*(seq + up3**2)
89  d222 = m2*(3.d0*seq + up2**2)
90  d223 = m3*(seq + up2**2)
91  d233 = m2*(seq + up3**2)
92  d333 = m3*(3.d0*seq + up3**2)
93 
94  m(5) = d11
95  m(6) = d12
96  m(7) = d13
97  m(8) = d22
98  m(9) = d23
99  m(10) = d33
100  m(11) = d111
101  m(12) = d112
102  m(13) = d113
103  m(14) = d122
104  m(15) = d123
105  m(16) = d133
106  m(17) = d222
107  m(18) = d223
108  m(19) = d233
109  m(20) = d333
110 
111  END SUBROUTINE collisions_istantaneous
112 
113 
114 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
115 ! C
116 ! Bathnagar-Gross-Krook (BGK) collision operator C
117 ! C
118 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
119 
120  SUBROUTINE collisions_bgk(M, Dt, tcol, dp, e)
121 
122  USE param1, only: small_number
123  USE constant, only: pi, ep_star
124 
125  IMPLICIT NONE
126 
127  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: M
128  DOUBLE PRECISION, INTENT(IN) :: Dt, tcol, dp, e
129 
130  DOUBLE PRECISION :: m0, m1, m2, m3, m11, m12, m13, m22, m23, m33
131  DOUBLE PRECISION :: m111, m112, m113, m122, m123, m133, m222
132  DOUBLE PRECISION :: m223, m233, m333
133 
134  DOUBLE PRECISION :: up1, up2, up3, b, a1, b1, gamma, om
135  DOUBLE PRECISION :: sig1, sig2, sig3, sig11, sig12, sig13, sig22, sig23
136  DOUBLE PRECISION :: sig33, seq, T
137 
138  DOUBLE PRECISION :: d11, d12, d13, d22, d23, d33
139  DOUBLE PRECISION :: d111, d112, d113, d122, d123, d133, d222, d223, d233, d333
140 
141  DOUBLE PRECISION :: kap, ALPHA_MAX
142 
143  m0 = m(1)
144  m1 = m(2)
145  m2 = m(3)
146  m3 = m(4)
147  m11 = m(5)
148  m12 = m(6)
149  m13 = m(7)
150  m22 = m(8)
151  m23 = m(9)
152  m33 = m(10)
153  m111 = m(11)
154  m112 = m(12)
155  m113 = m(13)
156  m122 = m(14)
157  m123 = m(15)
158  m133 = m(16)
159  m222 = m(17)
160  m223 = m(18)
161  m233 = m(19)
162  m333 = m(20)
163 
164  up1 = m1/m0
165  up2 = m2/m0
166  up3 = m3/m0
167  sig1 = max(0.d0, m11/m0 - up1**2)
168  sig2 = max(0.d0, m22/m0 - up2**2)
169  sig3 = max(0.d0, m33/m0 - up3**2)
170  t = ( sig1 + sig2 + sig3 )/3.d0
171 
172  !A.P. This parameter is hardcoded:
173  ! b = 0 for BGK, b = -1/2 for ES-BGK
174  b = 0.d0
175  gamma = 1.d0 - b
176 
177  om = (1 + e)/2.d0
178  a1 = gamma*(om**2)
179  b1 = a1 - 2.d0*gamma*om + 1.d0
180 
181  sig11 = a1*t + b1*sig1
182  sig22 = a1*t + b1*sig2
183  sig33 = a1*t + b1*sig3
184 
185  sig12 = b1*(m12/m0 - up1*up2)
186  sig13 = b1*(m13/m0 - up1*up3)
187  sig23 = b1*(m23/m0 - up2*up3)
188 
189  seq = ( sig11 + sig22 + sig33 )/3.d0
190 
191  d11 = m0*(sig11 + up1**2)
192  d12 = m0*(sig12 + up1*up2)
193  d13 = m0*(sig13 + up1*up3)
194  d22 = m0*(sig22 + up2**2)
195  d23 = m0*(sig23 + up2*up3)
196  d33 = m0*(sig33 + up3**2)
197 
198  d111 = -(-3.d0*d11 + 2.d0*m0*up1**2)*up1
199  d112 = -(-d11*up2 - (-2.d0*d12 + 2.d0*up1*up2*m0)*up1)
200 
201  d113 = -(-d11*up3+(-2.d0*d13+2.d0*up1*up3*m0)*up1)
202  d122 = -(-2.d0*d12*up2+(2.d0*(up2**2)*m0-d22)*up1)
203  d123 = -(-d12*up3-d13*up2+(2.d0*up2*up3*m0-d23)*up1)
204  d133 = -(-2.d0*d12*up3+(2.d0*(up3**2)*m0-d33)*up1)
205  d222 = -(-3.d0*d22+2.d0*(up2**2)*m0)*up2
206  d223 = -(-d22*up2+(-2.d0*d23+2.d0*up2*up3*m0)*up2)
207  d233 = -(-2.d0*d23*up3+(2.d0*(up3**2)*m0-d33)*up2)
208  d333 = -(-3.d0*d33+2.d0*(up3**2)*m0)*up3
209 
210 ! JEC: Probable BUG: alpha_max was not defined. include definition as found in
211 ! other places
212  alpha_max = 1. - ep_star
213 
214  IF (m0 > alpha_max - small_number) THEN
215  kap = 0.
216  ELSE
217  kap = exp(-12.0*radial_g0(m0)*m0*dt*sqrt(t/pi)/(dp*gamma))
218  END IF
219 
220  m(5) = d11 - kap*( d11 - m11 )
221  m(6) = d12 - kap*( d12 - m12 )
222  m(7) = d13 - kap*( d13 - m13 )
223  m(8) = d22 - kap*( d22 - m22 )
224  m(9) = d23 - kap*( d23 - m23 )
225  m(10) = d33 - kap*( d33 - m33 )
226  m(11) = d111 - kap*( d111 - m111 )
227  m(12) = d112 - kap*( d112 - m112 )
228  m(13) = d113 - kap*( d113 - m113 )
229  m(14) = d122 - kap*( d122 - m122 )
230  m(15) = d123 - kap*( d123 - m123 )
231  m(16) = d133 - kap*( d133 - m133 )
232  m(17) = d222 - kap*( d222 - m222 )
233  m(18) = d223 - kap*( d223 - m223 )
234  m(19) = d233 - kap*( d233 - m233 )
235  m(20) = d333 - kap*( d333 - m333 )
236  END SUBROUTINE collisions_bgk
237 
238 
239 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
240 ! C
241 ! Collision time C
242 ! C
243 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
244 
245  SUBROUTINE compute_collision_time (M, dp, theta, tcol, dt)
246 
247  USE param1, only: small_number
248  USE constant, only: pi, ep_star
249 ! JEC: should ep_star be replaced by ep_star_array?
250 
251  IMPLICIT NONE
252 
253  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NMOM) :: M
254  DOUBLE PRECISION, INTENT(IN) :: theta, dp, dt
255  DOUBLE PRECISION, INTENT(OUT) :: tcol
256 
257  DOUBLE PRECISION :: ALPHA_MAX
258 
259  !A.P.: Assuming m0 = M(1) = solids volume fraction due to the moments scaling
260  alpha_max = 1. - ep_star
261 
262  IF (m(1) > alpha_max - small_number) THEN
263  tcol = dt ! If this condition is satisfied, moments are set to
264  ! equilibrium. No need to integrate on more than one time step
265  ELSE
266  tcol = dp*sqrt(pi/(theta+small_number))/(12.0*m(1)*radial_g0(m(1))+small_number)
267  END IF
268  END SUBROUTINE compute_collision_time
269 
270  ! Global collision time based on the mean volume fraction (1 specie)
271  !SUBROUTINE GLOBAL_COLLISION_TIME (alpha, dp, theta, tcol, dt)
272  ! DOUBLE PRECISION, INTENT(IN) :: alpha
273  ! DOUBLE PRECISION, INTENT(IN) :: theta, dp, dt
274  ! DOUBLE PRECISION, INTENT(OUT) :: tcol
275 
276  ! DOUBLE PRECISION :: ALPHA_MAX
277 
278  ! ALPHA_MAX = 1. - EP_STAR
279 
280  ! IF (alpha > ALPHA_MAX - SMALL_NUMBER) THEN
281  ! tcol = dt ! If this condition is satisfied, moments are set to
282  ! equilibrium. No need to integrate on more than one time step
283  ! ELSE
284  ! tcol = dp*SQRT(pi/(theta+SMALL_NUMBER))/(12.0*alpha*RADIAL_G0(M(1))+SMALL_NUMBER)
285  ! END IF
286  !END SUBROUTINE
287 
288 
289 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
290 ! C
291 ! Boltzmann collision operator for the monodispersed case C
292 ! C
293 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
294 
295  SUBROUTINE collisions_boltzmann_one_specie (N, U, V, W, dp, e, Coll)
296 
297  USE constant, only: pi
298 
299  IMPLICIT NONE
300 
301  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: N, U, V, W
302  DOUBLE PRECISION, INTENT(IN) :: dp, e
303  DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NMOM) :: Coll
304 
305  DOUBLE PRECISION, DIMENSION(3) :: a1, a, vr, v1r, v2r, g1, g2
306 
307  DOUBLE PRECISION :: C200, C110, C101, C020, C011, C002, C300, C210, C201
308  DOUBLE PRECISION :: C120, C111, C102, C030, C021, C012, C003
309  DOUBLE PRECISION :: d, g, gsqr, o, ni, ui, vi, wi, nj, uj, vj, wj, cross
310  DOUBLE PRECISION :: denom, L11, L12, L13, L21, L22, L23, L31, L32, L33
311 
312  INTEGER :: i, j
313 
314  a1 = (/0.0462,0.0971,0.8235/) ! chosen at random (cannot be colinear with vr)
315  a = 1/((a1(1))**2 + (a1(2))**2 + (a1(3))**2) * a1;
316 
317  c200 = 0.
318  c110 = 0.
319  c101 = 0.
320  c020 = 0.
321  c011 = 0.
322  c002 = 0.
323  c300 = 0.
324  c210 = 0.
325  c201 = 0.
326  c120 = 0.
327  c111 = 0.
328  c102 = 0.
329  c030 = 0.
330  c021 = 0.
331  c012 = 0.
332  c003 = 0.
333 
334  o = (1+e)/2
335  DO i= 1, qmomk_nn
336  ni = n(i)
337  ui = u(i)
338  vi = v(i)
339  wi = w(i)
340  v1r(1) = ui
341  v1r(2) = vi
342  v1r(3) = wi
343  DO j = 1, 8
344  IF (i == j) THEN
345  cycle ! vi - vj = 0 always => nothing to add to source terms
346  END IF
347  nj = n(j)
348  uj = u(j)
349  vj = v(j)
350  wj = w(j)
351 
352  v2r(1) = uj
353  v2r(2) = vj
354  v2r(3) = wj
355 
356  vr(1) = v1r(1) - v2r(1)
357  vr(2) = v1r(2) - v2r(2)
358  vr(3) = v1r(3) - v2r(3)
359 
360  g = sqrt( vr(1)**2 + vr(2)**2 + vr(3)**2 )
361 
362  g1(1) = vr(1)/g
363  g1(2) = vr(2)/g
364  g1(3) = vr(3)/g
365 
366  cross = a(1)*g1(1) + a(2)*g1(2) + a(3)*g1(3)
367 
368  g2(1) = a(1) - cross*g1(1)
369  g2(2) = a(2) - cross*g1(2)
370  g2(3) = a(3) - cross*g1(3)
371 
372  denom = sqrt( g2(1)**2 + g2(2)**2 + g2(3)**2 )
373 
374  l11 = g2(1)/denom
375  l12 = g2(2)/denom
376  l13 = g2(3)/denom
377  l21 = ( g1(3)*a(2)-g1(2)*a(3) )/denom
378  l22 = ( g1(1)*a(3)-g1(3)*a(1) )/denom
379  l23 = ( g1(2)*a(1)-g1(1)*a(2) )/denom
380  l31 = g1(1)
381  l32 = g1(2)
382  l33 = g1(3)
383 
384  gsqr = g**2
385 
386  c200 = c200 + ni*nj*(pi*gsqr*o*(4*g*o*(l31**2)+6*uj*l31-6*ui*l31+g*o*(l21**2)+g*o*(l11**2)))/6
387 
388  c110 = c110 + ni*nj*(pi*gsqr*o*(4*g*o*l31*l32+3*uj*l32-3*ui*l32+3*vj*l31-3*vi*l31+g*o*l21*l22 + &
389  g*o*l11*l12))/6
390 
391  c101 = c101 + ni*nj*(pi*gsqr*o*(4*g*o*l31*l33+3*uj*l33-3*ui*l33+3*wj*l31-3*wi*l31+g*o*l21*l23 + &
392  g*o*l11*l13))/6
393 
394  c020 = c020 + ni*nj*(pi*gsqr*o*(4*g*o*(l32**2)+6*vj*l32-6*vi*l32+g*o*(l22**2)+g*o*(l12**2)))/6
395 
396  c011 = c011 + ni*nj*(pi*gsqr*o*(4*g*o*l32*l33+3*vj*l33-3*vi*l33+3*wj*l32-3*wi*l32+g*o*l22*l23 + &
397  g*o*l12*l13))/6
398 
399  c002 = c002 + ni*nj*(pi*gsqr*o*(4*g*o*(l33**2)+6*wj*l33-6*wi*l33+g*o*(l23**2)+g*o*(l13**2)))/6
400 
401  c300 = c300 + ni*nj*(pi*gsqr*o*(uj+ui)*(4*g*o*(l31**2)+6*uj*l31-6*ui*l31+g*o*(l21**2)+g*o*(l11**2)))/4
402 
403  c210 = c210 + ni*nj*(pi*gsqr*o*(8*g*o*uj*l31*l32+8*g*o*ui*l31*l32+6*(uj**2)*l32-6*(ui**2)*l32 + &
404  4*g*o*vj*(l31**2)+4*g*o*vi*(l31**2)+12*uj*vj*l31-12*ui*vi*l31+2*g*o*uj*l21*l22 + &
405  2*g*o*ui*l21*l22+g*o*vj*(l21**2)+g*o*vi*(l21**2)+2*g*o*uj*l11*l12+2*g*o*ui*l11*l12 + &
406  g*o*vj*(l11**2)+g*o*vi*(l11**2)))/12
407 
408  c201 = c201 + ni*nj*(((8*pi*(g**3)*(o**2)*uj+8*pi*(g**3)*(o**2)*ui)*l31+6*pi*gsqr*o*(uj**2)-6*pi*gsqr*o*(ui**2))*l33 + &
409  (4*pi*(g**3)*(o**2)*wj+4*pi*(g**3)*(o**2)*wi)*(l31**2)+(12*pi*(g**2)*o*uj*wj-12*pi*(g**2)*o*ui*wi)*l31 + &
410  (2*pi*(g**3)*(o**2)*uj+2*pi*(g**3)*(o**2)*ui)*l21*l23+(pi*(g**3)*(o**2)*wj+pi*(g**3)*(o**2)*wi)*(l21**2) + &
411  (2*pi*(g**3)*(o**2)*uj+2*pi*(g**3)*(o**2)*ui)*l11*l13+(pi*(g**3)*(o**2)*wj+pi*(g**3)*(o**2)*wi)*(l11**2))/12
412 
413  c120 = c120 + ni*nj*(pi*gsqr*o*(4*g*o*uj*(l32**2)+4*g*o*ui*(l32**2)+8*g*o*vj*l31*l32+8*g*o*vi*l31*l32 + &
414  12*uj*vj*l32-12*ui*vi*l32+6*(vj**2)*l31-6*(vi**2)*l31+g*o*uj*(l22**2)+g*o*ui*(l22**2) + &
415  2*g*o*vj*l21*l22+2*g*o*vi*l21*l22+g*o*uj*(l12**2)+g*o*ui*(l12**2)+2*g*o*vj*l11*l12 + &
416  2*g*o*vi*l11*l12))/12
417 
418  c111 = c111 + ni*nj*(pi*gsqr*o*(4*g*o*uj*l32*l33+4*g*o*ui*l32*l33+4*g*o*vj*l31*l33+4*g*o*vi*l31*l33 + &
419  6*uj*vj*l33-6*ui*vi*l33+4*g*o*wj*l31*l32+4*g*o*wi*l31*l32+6*uj*wj*l32-6*ui*wi*l32 + &
420  6*vj*wj*l31-6*vi*wi*l31+g*o*uj*l22*l23+g*o*ui*l22*l23+g*o*vj*l21*l23+g*o*vi*l21*l23 + &
421  g*o*wj*l21*l22+g*o*wi*l21*l22+g*o*uj*l12*l13+g*o*ui*l12*l13+g*o*vj*l11*l13+g*o*vi*l11*l13 + &
422  g*o*wj*l11*l12+g*o*wi*l11*l12))/12
423 
424  c102 = c102 + ni*nj*(gsqr*o*pi*(4*g*o*vj*(l33**2)+4*g*o*vi*(l33**2)+8*g*o*wj*l32*l33+8*g*o*wi*l32*l33 + &
425  12*vj*wj*l33-12*vi*wi*l33+6*(wj**2)*l32-6*(wi**2)*l32+g*o*vj*(l23**2)+g*o*vi*(l23**2)+2*g*o*wj*l22*l23 + &
426  2*g*o*wi*l22*l23+g*o*vj*(l13**2)+g*o*vi*(l13**2)+2*g*o*wj*l12*l13+2*g*o*wi*l12*l13))/12
427 
428  c030 = c030 + ni*nj*(pi*gsqr*o*(vj+vi)*(4*g*o*(l32**2)+6*vj*l32-6*vi*l32+g*o*(l22**2)+g*o*(l12**2)))/4
429 
430  c021 = c021 + ni*nj*(pi*gsqr*o*(8*g*o*vj*l32*l33+8*g*o*vi*l32*l33+6*(vj**2)*l33-6*(vi**2)*l33 + &
431  4*g*o*wj*(l32**2)+4*g*o*wi*(l32**2)+12*vj*wj*l32-12*vi*wi*l32+2*g*o*vj*l22*l23+2*g*o*vi*l22*l23 + &
432  g*o*wj*(l22**2)+g*o*wi*(l22**2)+2*g*o*vj*l12*l13+2*g*o*vi*l12*l13+g*o*wj*(l12**2)+g*o*wi*(l12**2)))/12
433 
434  c012 = c012 + ni*nj*(gsqr*o*pi*(4*g*o*vj*(l33**2)+4*g*o*vi*(l33**2)+8*g*o*wj*l32*l33+8*g*o*wi*l32*l33 + &
435  12*vj*wj*l33-12*vi*wi*l33+6*(wj**2)*l32-6*(wi**2)*l32+g*o*vj*(l23**2)+g*o*vi*(l23**2)+2*g*o*wj*l22*l23 + &
436  2*g*o*wi*l22*l23+g*o*vj*(l13**2)+g*o*vi*(l13**2)+2*g*o*wj*l12*l13+2*g*o*wi*l12*l13))/12
437 
438  c003 = c003 + ni*nj*(pi*gsqr*o*(wj+wi)*(4*g*o*(l33**2)+6*wj*l33-6*wi*l33+g*o*(l23**2)+g*o*(l13**2)))/4
439  END DO
440  END DO
441 
442  d = dp*dp/2 ! Scaling
443  coll(1) = 0
444  coll(2) = 0
445  coll(3) = 0
446  coll(4) = 0
447  coll(5) = d*c200
448  coll(6) = d*c110
449  coll(7) = d*c101
450  coll(8) = d*c020
451  coll(9) = d*c011
452  coll(10) = d*c002
453  coll(11) = d*c300
454  coll(12) = d*c210
455  coll(13) = d*c201
456  coll(14) = d*c120
457  coll(15) = d*c111
458  coll(16) = d*c102
459  coll(17) = d*c030
460  coll(18) = d*c021
461  coll(19) = d*c012
462  coll(20) = d*c003
463 
464  END SUBROUTINE collisions_boltzmann_one_specie
465 
466 
467 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
468 ! C
469 ! Boltzmann collision operator for bidisperse case C
470 ! C
471 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
472 
473  SUBROUTINE collisions_boltzmann_two_species (N1, U1, V1, W1, N2, U2, V2, W2, m1, m2, dp1, dp2, e, Coll)
474 
475  USE constant, only: pi
476 
477  IMPLICIT NONE
478 
479  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: N1, U1, V1, W1, N2, U2, V2, W2
480  DOUBLE PRECISION, INTENT(IN) :: m1, m2, dp1, dp2, e
481  DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NMOM) :: Coll
482 
483  DOUBLE PRECISION, DIMENSION(3) :: a1, a, vr, v1r, v2r, g1, g2
484 
485  DOUBLE PRECISION :: C100, C010, C001, C200, C110, C101, C020, C011
486  DOUBLE PRECISION :: C002, C300, C210, C201, C120, C111, C102, C030, C021, C012, C003
487 
488  DOUBLE PRECISION :: dp, d, g, gsqr, o, mu, ni, ui, vi, wi, nj, uj, vj, wj
489  DOUBLE PRECISION :: cross, denom, L11, L12, L13, L21, L22, L23, L31, L32, L33
490 
491  INTEGER :: i, j
492 
493  a1 = (/0.0462,0.0971,0.8235/) ! chosen at random (cannot be colinear with vr)
494  a = 1/((a1(1))**2 + (a1(2))**2 + (a1(3))**2) * a1;
495 
496  c100 = 0.d0
497  c010 = 0.d0
498  c001 = 0.d0
499  c200 = 0.d0
500  c110 = 0.d0
501  c101 = 0.d0
502  c020 = 0.d0
503  c011 = 0.d0
504  c002 = 0.d0
505  c300 = 0.d0
506  c210 = 0.d0
507  c201 = 0.d0
508  c120 = 0.d0
509  c111 = 0.d0
510  c102 = 0.d0
511  c030 = 0.d0
512  c021 = 0.d0
513  c012 = 0.d0
514  c003 = 0.d0
515 
516  ! Reduced mass
517  mu = m1*m2/(m1+m2)
518  o = -(1+e)*mu/m1
519  ! Collision diameter
520  dp = (dp1 + dp2)/2
521  ! Scaling factor
522  d = (dp**2)/m2
523 
524  DO i = 1, qmomk_nn
525  ui = u1(i)
526  vi = v1(i)
527  wi = w1(i)
528  ni = n1(i)
529 
530  v1r(1) = ui
531  v1r(2) = vi
532  v1r(3) = wi
533 
534  DO j = 1, qmomk_nn
535  IF (i == j) THEN
536  cycle ! vi - vj = 0 always => nothing to add to source term
537  END IF
538 
539  nj = n2(j)
540  uj = u2(j)
541  vj = v2(j)
542  wj = w2(j)
543 
544  v2r(1) = uj
545  v2r(2) = vj
546  v2r(3) = wj
547 
548  vr(1) = v1r(1) - v2r(1)
549  vr(2) = v1r(2) - v2r(2)
550  vr(3) = v1r(3) - v2r(3)
551 
552  g = sqrt(vr(1)*vr(1) + vr(2)*vr(2) + vr(3)*vr(3))
553 
554  IF (g == 0) THEN ! Avoiding division by zero in the transformation matrix
555  cycle
556  END IF
557 
558  g1(1) = vr(1)/g
559  g1(2) = vr(2)/g
560  g1(3) = vr(3)/g
561 
562  cross = a(1)*g1(1) + a(2)*g1(2) + a(3)*g1(3)
563 
564  g2(1) = a(1) - cross*g1(1)
565  g2(2) = a(2) - cross*g1(2)
566  g2(3) = a(3) - cross*g1(3)
567 
568  denom = sqrt( g2(1)*g2(1) + g2(2)*g2(2) + g2(3)*g2(3))
569 
570  l11 = g2(1)/denom
571  l12 = g2(2)/denom
572  l13 = g2(3)/denom
573  l21 = ( g1(3)*a(2)-g1(2)*a(3) )/denom
574  l22 = ( g1(1)*a(3)-g1(3)*a(1) )/denom
575  l23 = ( g1(2)*a(1)-g1(1)*a(2) )/denom
576  l31 = g1(1)
577  l32 = g1(2)
578  l33 = g1(3)
579 
580  gsqr = g**2
581 
582  c100 = c100 + ni*nj*(pi*gsqr*o*l31)/2
583  c010 = c010 + ni*nj*(pi*gsqr*o*l32)/2
584  c001 = c001 + ni*nj*(pi*gsqr*o*l33)/2
585 
586  c200 = c200 + ni*nj*(4*pi*(g**3)*(o**2)*(l31**2)+12*pi*gsqr*o*ui*l31 &
587  + pi*(g**3)*(o**2)*(l21**2)+pi*(g**3)*(o**2)*(l11**2))/12
588  c110 = c110 + ni*nj*((4*pi*(g**3)*(o**2)*l31+6*pi*gsqr*o*ui)*l32+6*pi*gsqr*o*vi*l31+pi*(g**3)*(o**2)*l21*l22 + &
589  pi*(g**3)*(o**2)*l11*l12)/12
590 
591  c101 = c101 + ni*nj*((4*pi*(g**3)*(o**2)*l31+6*pi*gsqr*o*ui)*l33+6*pi*gsqr*o*wi*l31+pi*(g**3)*(o**2)*l21*l23 + &
592  pi*(g**3)*(o**2)*l11*l13)/12
593 
594  c020 = c020 + ni*nj*(4*pi*(g**3)*(o**2)*(l32**2)+12*pi*gsqr*o*vi*l32 &
595  + pi*(g**3)*(o**2)*(l22**2)+pi*(g**3)*(o**2)*(l12**2))/12
596 
597  c011 = c011 + ni*nj*((4*pi*(g**3)*(o**2)*l32+6*pi*gsqr*o*vi)*l33+6*pi*gsqr*o*wi*l32+pi*(g**3)*(o**2)*l22*l23 + &
598  pi*(g**3)*(o**2)*l12*l13)/12
599 
600  c002 = c002 + ni*nj*(4*pi*(g**3)*(o**2)*(l33**2)+12*pi*gsqr*o*wi*l33 &
601  + pi*(g**3)*(o**2)*(l23**2)+pi*(g**3)*(o**2)*(l13**2))/12
602 
603  c300 = c300 + ni*nj*(2*pi*(g**4)*(o**3)*(l31**3)+8*pi*(g**3)*(o**2)*ui*(l31**2) &
604  + (pi*(g**4)*(o**3)*(l21**2)+pi*(g**4)*(o**3)*(l11**2) + 12*pi*gsqr*o*(ui**2))*l31 &
605  + 2*pi*(g**3)*(o**2)*ui*(l21**2)+2*pi*(g**3)*(o**2)*ui*(l11**2))/8
606 
607  c210 = c210 + ni*nj*((6*pi*(g**4)*(o**3)*(l31**2)+16*pi*(g**3)*(o**2)*ui*l31 &
608  + pi*(g**4)*(o**3)*(l21**2)+pi*(g**4)*(o**3)*(l11**2) + 12*pi*gsqr*o*(ui**2))*l32 &
609  + 8*pi*(g**3)*(o**2)*vi*(l31**2)+(2*pi*(g**4)*(o**3)*l21*l22+2*pi*(g**4)*(o**3)*l11*l12 &
610  + 24*pi*(g**2)*o*ui*vi)*l31+4*pi*(g**3)*(o**2)*ui*l21*l22+2*pi*(g**3)*(o**2)*vi*(l21**2) &
611  + 4*pi*(g**3)*(o**2)*ui*l11*l12+2*pi*(g**3)*(o**2)*vi*(l11**2))/(24)
612 
613  c201 = c201 + ni*nj*((6*pi*(g**4)*(o**3)*(l31**2)+16*pi*(g**3)*(o**2)*ui*l31 &
614  + pi*(g**4)*(o**3)*(l21**2)+pi*(g**4)*(o**3)*(l11**2) + 12*pi*gsqr*o*(ui**2))*l33 &
615  + 8*pi*(g**3)*(o**2)*wi*(l31**2)+(2*pi*(g**4)*(o**3)*l21*l23+2*pi*(g**4)*(o**3)*l11*l13 &
616  + 24*pi*gsqr*o*ui*wi)*l31+4*pi*(g**3)*(o**2)*ui*l21*l23+2*pi*(g**3)*(o**2)*wi*(l21**2) &
617  + 4*pi*(g**3)*(o**2)*ui*l11*l13+2*pi*(g**3)*(o**2)*wi*(l11**2))/(24)
618 
619  c120 = c120 + ni*nj*((6*pi*(g**4)*(o**3)*l31+8*pi*(g**3)*(o**2)*ui)*(l32**2)+(16*pi*(g**3)*(o**2)*vi*l31 + &
620  2*pi*(g**4)*(o**3)*l21*l22+2*pi*(g**4)*(o**3)*l11*l12+24*pi*gsqr*o*ui*vi)*l32+(pi*(g**4)*(o**3)*(l22**2) + &
621  pi*(g**4)*(o**3)*(l12**2)+12*pi*gsqr*o*(vi**2))*l31+2*pi*(g**3)*(o**2)*ui*(l22**2)+4*pi*(g**3)*(o**2)*vi*l21*l22 + &
622  2*pi*(g**3)*(o**2)*ui*(l12**2)+4*pi*(g**3)*(o**2)*vi*l11*l12)/(24)
623 
624  c111 = c111 + ni*nj*(((6*pi*(g**4)*(o**3)*l31+8*pi*(g**3)*(o**2)*ui)*l32 &
625  + 8*pi*(g**3)*(o**2)*vi*l31+pi*(g**4)*(o**3)*l21*l22 + pi*(g**4)*(o**3)*l11*l12 &
626  + 12*pi*gsqr*o*ui*vi)*l33+(8*pi*(g**3)*(o**2)*wi*l31+pi*(g**4)*(o**3)*l21*l23 + &
627  pi*(g**4)*(o**3)*l11*l13+12*pi*gsqr*o*ui*wi)*l32+(pi*(g**4)*(o**3)*l22*l23+pi*(g**4)*(o**3)*l12*l13+ &
628  12*pi*gsqr*o*vi*wi)*l31+(2*pi*(g**3)*(o**2)*ui*l22+2*pi*(g**3)*(o**2)*vi*l21)*l23 + &
629  2*pi*(g**3)*(o**2)*wi*l21*l22+(2*pi*(g**3)*(o**2)*ui*l12+2*pi*(g**3)*(o**2)*vi*l11)*l13 + &
630  2*pi*(g**3)*(o**2)*wi*l11*l12)/(24)
631 
632  c102 = c102 + ni*nj*((6*pi*(g**4)*(o**3)*l31+8*pi*(g**3)*(o**2)*ui)*(l33**2)+(16*pi*(g**3)*(o**2)*wi*l31 + &
633  2*pi*(g**4)*(o**3)*l21*l23+2*pi*(g**4)*(o**3)*l11*l13+24*pi*(g**2)*o*ui*wi)*l33 + &
634  (pi*(g**4)*(o**3)*(l23**2)+pi*(g**4)*(o**3)*(l13**2)+12*pi*(g**2)*o*(wi**2))*l31+2*pi*(g**3)*(o**2)*ui*(l23**2)+ &
635  4*pi*(g**3)*(o**2)*wi*l21*l23+2*pi*(g**3)*(o**2)*ui*(l13**2)+4*pi*(g**3)*(o**2)*wi*l11*l13)/(24)
636 
637  c030 = c030 + ni*nj*(2*pi*(g**4)*(o**3)*(l32**3)+8*pi*(g**3)*(o**2)*vi*(l32**2) &
638  + (pi*(g**4)*(o**3)*(l22**2)+pi*(g**4)*(o**3)*(l12**2) + 12*pi*(g**2)*o*(vi**2))*l32 &
639  + 2*pi*(g**3)*(o**2)*vi*(l22**2)+2*pi*(g**3)*(o**2)*vi*(l12**2))/8
640 
641  c021 = c021 + ni*nj*((6*pi*(g**4)*(o**3)*(l32**2)+16*pi*(g**3)*(o**2)*vi*l32 &
642  + pi*(g**4)*(o**3)*(l22**2)+pi*(g**4)*(o**3)*(l12**2) + 12*pi*gsqr*o*(vi**2))*l33 &
643  + 8*pi*(g**3)*(o**2)*wi*(l32**2)+(2*pi*(g**4)*(o**3)*l22*l23+2*pi*(g**4)*(o**3)*l12*l13 &
644  + 24*pi*(g**2)*o*vi*wi)*l32+4*pi*(g**3)*(o**2)*vi*l22*l23+2*pi*(g**3)*(o**2)*wi*(l22**2) &
645  + 4*pi*(g**3)*(o**2)*vi*l12*l13+2*pi*(g**3)*(o**2)*wi*(l12**2))/(24)
646 
647  c012 = c012 + ni*nj*((6*pi*(g**4)*(o**3)*l32+8*pi*(g**3)*(o**2)*vi)*(l33**2)+(16*pi*(g**3)*(o**2)*wi*l32 + &
648  2*pi*(g**4)*(o**3)*l22*l23+2*pi*(g**4)*(o**3)*l12*l13+24*pi*gsqr*o*vi*wi)*l33+(pi*(g**4)*(o**3)*(l23**2) + &
649  pi*(g**4)*(o**3)*(l13**2)+12*pi*gsqr*o*(wi**2))*l32+2*pi*(g**3)*(o**2)*vi*(l23**2)+4*pi*(g**3)*(o**2)*wi*l22*l23 + &
650  2*pi*(g**3)*(o**2)*vi*(l13**2)+4*pi*(g**3)*(o**2)*wi*l12*l13)/(24)
651 
652  c003 = c003 + ni*nj*(2*pi*(g**4)*(o**3)*(l33**3)+8*pi*(g**3)*(o**2)*wi*(l33**2) &
653  + (pi*(g**4)*(o**3)*(l23**2)+pi*(g**4)*(o**3)*(l13**2) + 12*pi*gsqr*o*(wi**2))*l33 &
654  + 2*pi*(g**3)*(o**2)*wi*(l23**2)+2*pi*(g**3)*(o**2)*wi*(l13**2))/8;
655 
656  END DO
657  END DO
658 
659  coll(1) = 0;
660  coll(2) = d*c100;
661  coll(3) = d*c010;
662  coll(4) = d*c001;
663  coll(5) = d*c200;
664  coll(6) = d*c110;
665  coll(7) = d*c101;
666  coll(8) = d*c020;
667  coll(9) = d*c011;
668  coll(10) = d*c002;
669  coll(11) = d*c300;
670  coll(12) = d*c210;
671  coll(13) = d*c201;
672  coll(14) = d*c120;
673  coll(15) = d*c111;
674  coll(16) = d*c102;
675  coll(17) = d*c030;
676  coll(18) = d*c021;
677  coll(19) = d*c012;
678  coll(20) = d*c003;
679 
680  END SUBROUTINE collisions_boltzmann_two_species
681 
682 
683 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
684 ! C
685 ! Calculates the rate of change in the moments due to collisions C
686 ! (1 specie) C
687 ! C
688 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
689 
690  SUBROUTINE solve_boltzmann_collisions_one_specie(M, N, U, V, W, dt, &
691  e, dp, order)
692 
693  IMPLICIT NONE
694 
695  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: M
696  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: N, U, V, W
697  DOUBLE PRECISION, INTENT(IN) :: dt, e, dp
698  INTEGER, INTENT(IN) :: order
699 
700  DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: Coll, Mtmp
701  DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Ntmp, Utmp, Vtmp, Wtmp
702  DOUBLE PRECISION :: h
703 
704 
705  IF (order == 0) THEN ! Euler method - First order
706  h = dt
707  CALL collisions_boltzmann_one_specie (n, u, v, w, dp, e, coll)
708  mtmp = m + h*coll
709  ELSE IF (order == 1) THEN ! Runge - Kutta second order
710  h = dt ;
711  CALL collisions_boltzmann_one_specie (n, u, v, w, dp, e, coll)
712  mtmp = m + 0.5*h*coll
713  CALL eight_node_3d (mtmp, ntmp, utmp, vtmp, wtmp) ! Projection step
714  CALL collisions_boltzmann_one_specie (ntmp, utmp, vtmp, wtmp, dp, e, coll)
715  mtmp = m + h*coll
716  ELSE
717  print *,'QMOMK: Wrong order in Boltzmann collisions'
718  error_stop
719  END IF
720  m= mtmp
721  END SUBROUTINE solve_boltzmann_collisions_one_specie
722 
723 
724 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
725 ! C
726 ! Calculates the rate of change in the moments due to collisions C
727 ! (2 specie) C
728 ! C
729 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
730 
731  SUBROUTINE solve_boltzmann_collisions_two_species(M1, N1, U1, V1, &
732  w1, m2, n2, u2, v2, w2, dt, m_1, m_2, dp1, dp2, e11, e12, order)
733 
734  IMPLICIT NONE
735 
736  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: M1
737  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NMOM) :: M2
738  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: N1, U1, V1, W1, N2, U2, V2, W2
739  DOUBLE PRECISION, INTENT(IN) :: dt, dp1, dp2, m_1, m_2, e11, e12
740  INTEGER, INTENT(IN) :: order
741 
742  DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: Coll, Colltmp, M1tmp
743  DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: N1tmp, U1tmp, V1tmp, W1tmp
744  DOUBLE PRECISION :: h
745 
746  IF (order == 0) THEN ! Euler method - First order
747  h = dt ;
748  ! Collisions inside the same specie
749  CALL collisions_boltzmann_two_species (n1, u1, v1, w1, n1, u1, v1, w1, m_1, m_1, dp1, dp1, e11, colltmp)
750  coll = colltmp
751  ! Collisions between the two species
752  CALL collisions_boltzmann_two_species (n1, u1, v1, w1, n2, u2, v2, w2, m_1, m_2, dp1, dp2, e12, colltmp)
753  coll = coll + colltmp
754  m1tmp = m1 + h*coll
755  ELSE IF (order == 1) THEN ! Runge-Kutta second order method
756  h = dt ;
757  ! Collisions inside the same specie
758  CALL collisions_boltzmann_two_species (n1, u1, v1, w1, n1, u1, v1, w1, m_1, m_1, dp1, dp1, e11, colltmp)
759  coll = colltmp
760  ! Collisions between the two species
761  CALL collisions_boltzmann_two_species (n1, u1, v1, w1, n2, u2, v2, w2, m_1, m_2, dp1, dp2, e12, colltmp)
762  coll = coll + colltmp
763  m1tmp = m1 + 0.5*h*coll ;
764  CALL eight_node_3d (m1tmp, n1tmp, u1tmp, v1tmp, w1tmp) ! Projection step
765  ! Collisions inside the same specie
766 
767  CALL collisions_boltzmann_two_species (n1tmp, u1tmp, v1tmp, w1tmp, &
768  n1tmp, u1tmp, v1tmp, w1tmp, m_1, m_1, dp1, dp1, e11, colltmp)
769 
770  coll = colltmp
771  ! Collisions between the two species
772  CALL collisions_boltzmann_two_species (n1tmp, u1tmp, v1tmp, w1tmp, n2, u2, v2, w2, m_1, m_2, dp1, dp2, e12, colltmp)
773  coll = coll + colltmp
774  m1tmp = m1 + h*coll;
775  ELSE
776  print *,'QMOMK: Wrong order in Boltzmann collisions'
777  error_stop
778  END IF
779  m1 = m1tmp;
780  END SUBROUTINE solve_boltzmann_collisions_two_species
781 
782 
783 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
784 ! C
785 ! C
786 ! C
787 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
788 
789  DOUBLE PRECISION FUNCTION radial_g0(ALPHA)
790 
791  USE param1, only: small_number, large_number
792  USE constant, only: ep_star
793 
794  IMPLICIT NONE
795 
796  DOUBLE PRECISION, INTENT(IN) :: ALPHA
797 
798  DOUBLE PRECISION ALPHA_MAX
799 
800  alpha_max = 1. - ep_star
801 
802  radial_g0 = 1.
803  IF (alpha < alpha_max - small_number) THEN
804  radial_g0 = 1./(1.0-alpha) + 3.0*alpha/(2.*(1.-alpha)**2) + (alpha**2)/(2*(1.0-alpha)**3)
805  ELSE
806  radial_g0 = large_number
807  END IF
808 
809  END FUNCTION radial_g0
810 
811 END MODULE qmomk_collision
qmomk_collision::solve_boltzmann_collisions_one_specie
subroutine, public solve_boltzmann_collisions_one_specie(M, N, U, V, W, dt, e, dp, order)
Definition: qmomk_collision_mod.f:692
qmomk_parameters
Definition: qmomk_parameters_mod.f:1
qmomk_collision::collisions_boltzmann_one_specie
subroutine collisions_boltzmann_one_specie(N, U, V, W, dp, e, Coll)
Definition: qmomk_collision_mod.f:296
qmomk_collision::collisions_bgk
subroutine, public collisions_bgk(M, Dt, tcol, dp, e)
Definition: qmomk_collision_mod.f:121
qmomk_collision::collisions_istantaneous
subroutine, public collisions_istantaneous(M, dp)
Definition: qmomk_collision_mod.f:47
qmomk_quadrature::eight_node_3d
subroutine, public eight_node_3d(mom, n, u, v, w)
Definition: qmomk_quadrature_mod.f:220
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
qmomk_collision::compute_collision_time
subroutine, public compute_collision_time(M, dp, theta, tcol, dt)
Definition: qmomk_collision_mod.f:246
qmomk_collision::radial_g0
double precision function, public radial_g0(ALPHA)
Definition: qmomk_collision_mod.f:790
param1::large_number
double precision, parameter large_number
Definition: param1_mod.f:23
constant::ep_star
double precision ep_star
Definition: constant_mod.f:29
qmomk_collision
Definition: qmomk_collision_mod.f:19
qmomk_quadrature
Definition: qmomk_quadrature_mod.f:16
qmomk_collision::collisions_boltzmann_two_species
subroutine collisions_boltzmann_two_species(N1, U1, V1, W1, N2, U2, V2, W2, m1, m2, dp1, dp2, e, Coll)
Definition: qmomk_collision_mod.f:474
constant::pi
double precision, parameter pi
Definition: constant_mod.f:158
qmomk_parameters::qmomk_nn
integer, parameter qmomk_nn
Definition: qmomk_parameters_mod.f:11
qmomk_collision::solve_boltzmann_collisions_two_species
subroutine, public solve_boltzmann_collisions_two_species(M1, N1, U1, V1, W1, M2, N2, U2, V2, W2, dt, m_1, m_2, dp1, dp2, e11, e12, order)
Definition: qmomk_collision_mod.f:733