modified: src/electrostatics.f90

	modified:   src/particlesTransport.f90
	modified:   src/var.f90

src/electrostatics.f90

@@ -291,7 +291,7 @@
 ! pipe:
 !        ut2=up(n)**2
 !        un2=vp(n)**2+wp(n)**2
-
+!
       if (direction.eq.1) then
         ut2=vp(n)**2+wp(n)**2
         un2=up(n)**2
@@ -304,9 +304,10 @@
       endif
 
 ! John (1979)
-      alpha1=(0.625_pr*pi*rhop*(1._pr+restRatio)*un2* &
-            ((1._pr-nyw**2)/Ew+(1._pr-nyp**2)/Ep))**(0.4_pr)
-      AoAtot=alpha1/4._pr
+!      alpha1=(0.625_pr*pi*rhop*(1._pr+restRatio)*un2* &
+!            ((1._pr-nyw**2)/Ew+(1._pr-nyp**2)/Ep))**(0.4_pr)
+!      AoAtot=alpha1/4._pr
+      AoAtot=1._pr
 
       dqp=Qaccfactor*AoAtot*(qpmax-q_el(n))
       q_el(n)=min(dqp+q_el(n),qpmax)
@@ -331,17 +332,21 @@
       integer :: n1,n2
       character(70) :: filename
 
-      urel2=(up(n1)-up(n2))**2+(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2
-
-      if (urel2.eq.0._pr) then
-        dqp=0.5_pr*(q_el(n2)-q_el(n1))
-      else
-! Soo (1971)
-        alpha1=radp(n1)*radp(n2)*(0.625_pr*pi*rhop*(1._pr+restRatio)*urel2 &
-               *(radp(n1)+radp(n2))**(0.5_pr)/(radp(n1)**3+radp(n2)**3) &
-               *(1._pr-nyp**2)/Ep)**(0.4_pr)
-        dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1))
-      endif
+!      urel2=(up(n1)-up(n2))**2+(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2
+!
+!      if (urel2.eq.0._pr) then
+!        dqp=0.5_pr*(q_el(n2)-q_el(n1))
+!      else
+!! Soo (1971)
+!        alpha1=radp(n1)*radp(n2)*(0.625_pr*pi*rhop*(1._pr+restRatio)*urel2 &
+!               *(radp(n1)+radp(n2))**(0.5_pr)/(radp(n1)**3+radp(n2)**3) &
+!               *(1._pr-nyp**2)/Ep)**(0.4_pr)
+!        dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1))
+!      endif
+
+      AoAtot=1._pr
+
+      dqp=Qaccfactor*AoAtot*(q_el(n2)-q_el(n1))
 
       q_el(n1)=q_el(n1)+dqp
       q_el(n2)=q_el(n2)-dqp

src/particlesTransport.f90

@@ -393,7 +393,7 @@
       ppcollnum(n1)=ppcollnum(n1)+1
       ppcollnum(n2)=ppcollnum(n2)+1
       numcol=numcol+1  
-!      if (q_el(n1).ne.q_el(n2)) call chargeParticleParticle(n1,n2)
+      if (q_el(n1).ne.q_el(n2)) call chargeParticleParticle(n1,n2)
 
 ! fictitious contact point
       xpc(n1)=xp(n1)+fracdt*up(n1)*dt

src/var.f90

@@ -27,7 +27,7 @@
       Ep=1.e8_pr, &                 !< particle Young's modulus (kg/s**2/m)
       nyp=0.4_pr , &                !< particle Poisson ratio
       nyw=0.28_pr, &                !< duct Poisson ratio
-      Qaccfactor= 100._pr           !< artificially accelerate the charging rate
+      Qaccfactor= 0.1_pr            !< artificially accelerate the charging rate
       
       character(70) :: version='pafiX v1.1.0 (Copyright 2019 by H. Grosshans)'