Skip to content
Snippets Groups Projects
Commit 4846249a authored by Gizem Özler's avatar Gizem Özler
Browse files

updated charging model

parent 10591331
Branches
No related tags found
1 merge request!5Master
Hide whitespace changes
Inline Side-by-side
input/input_example.dat
+ 2
− 2
View file @ 4846249a
......@@ -26,7 +26,7 @@ particle radius min, max, step (m,m,m), material density (kg/m**3), restitution
10.E-6,140.E-6,10.E-6,2370.,1.0
particle charge: initial, equilibrium (C,C) =
0.,1.E-12
charging model [(r)andom,(c)ondenser,(p)atch], initial TCS density (#/m**2), charging acceleration (-) =
p,1.E13,10.
charging model [(r)andom,(c)ondenser,(p)atch,(s)ingle], initial TCS density (#/m**2), charging acceleration (-) =
s,1.E13,10.
gravity vector in x,y,z-dir. (m/s**2,m/s**2,m/s**2) =
9.81,0,0
src/electrostatics.f90
+ 107
− 0
View file @ 4846249a
......@@ -356,3 +356,110 @@
if (Acmax.gt.(4*pi*radp(n2)**2)) print *, 'contact area larger than the area of the n2 particle'
end
!####################################################################
!> @author Gizem Ozler
!> @brief compute particle-wall charging, single particle charging model
subroutine chargeParticleWallSingle(n,direction)
use var
real(kind=pr) :: ut2,un2,dqp,random(2),random_normal,Efp,alpha_rat,A_rat,Ac_rat,Nw_rat,mu,sigma
integer :: n,direction
character(70) :: filename
if (first_call_particle(n)) then ! first collision of the particle "n"
c_rat(n) = 1._pr ! c1/c0 = 1
c_rat_prev(n) = 1._pr ! c0/c0 = 1
first_call_particle(n) = .false.
else
c_rat_prev(n)= c_rat(n)
c_rat(n)= 1._pr - term_prev(n) ! eqn. 30
endif
if (direction.eq.1) then
ut2=vp(n)**2+wp(n)**2
un2=up(n)**2
elseif (direction.eq.2) then
ut2=up(n)**2+wp(n)**2
un2=vp(n)**2
elseif (direction.eq.3) then
ut2=up(n)**2+vp(n)**2
un2=wp(n)**2
endif
Efp=8._pr*pi*eps_el*q_elNext(n)/radp(n)**2 ! eqn. 36
alpha_rat=1-Efp/Esat ! eqn. 38
A_rat=((radp(n)/rad0)**5*sqrt(un2)/un0)**0.4_pr ! eqn. 24
Ac_rat=(5._pr*pi*rhop*un2*(1._pr-nyp**2)/Ep/128._pr)**(0.4_pr) ! eqn. 33
Nw_rat=A_rat*alpha_rat ! eqn. 22
Np_rat(n)=A_rat*alpha_rat*c_rat(n) ! eqn. 22
mu=mu0*sqrt(Np_rat(n)) ! eqn. 9
sigma=dq0_min*Nw_rat+(sigma0-dq0_min)*Np_rat(n) ! eqn. 18
call random_number(random)
random_normal=cos(2._pr*pi*random(1))*sqrt(-2._pr*log(random(2)+1.e-12_pr)) ! Box-Muller method
dqp=mu*random_normal+sigma ! ...(mu*random_normal+sigma)
q_elNext(n)=max(min(q_elNext(n)+dqp,qpmax),qp0)
term_prev(n)=term_prev(n) + c_rat_prev(n)*Ac_rat*dqp/mu ! eqn. 30
write(filename,'(a,i3.3,a)') 'results/particlesWall_p',myid,'.dat'
open(access='append',unit=10,file=filename)
write(10,77) nt,t,sqrt(ut2),sqrt(un2),dqp,q_elNext(n),c_rat(n)
77 format(i7,6(x,es11.3e2))
close(10)
end
!####################################################################
!> @author Gizem Ozler
!> @brief compute particle-particle charging, single particle charging model
subroutine chargeParticleParticleSingle(n1,n2)
use var
real(kind=pr) dqp,random(2),random_normal,Efp,alpha_rat,urel2,A_rat,Ac_rat,mu,sigma
integer :: n1,n2
character(70) :: filename
if (first_call_particle(n1)) then ! first collision of particles "n1, n2"
first_call_particle(n1) = .false.
c_rat(n1) = 1._pr ! c1/c0 = 1
c_rat_prev(n1) = 1._pr ! c0/c0 = 1
else
c_rat_prev(n1)= c_rat(n1)
c_rat(n1)= 1._pr - term_prev(n1) ! eqn. 30
endif
if (first_call_particle(n2)) then
first_call_particle(n2) = .false.
c_rat(n2) = 1._pr ! c1/c0 = 1
c_rat_prev(n2) = 1._pr ! c0/c0 = 1
else
c_rat_prev(n2)= c_rat(n2)
c_rat(n2)= 1._pr - term_prev(n2) ! eqn. 30
endif
Efp=4._pr*pi*eps_el*(q_elNext(n2)/radp(n2)**2-q_elNext(n1)/radp(n1)**2) ! eqn. 37
alpha_rat=1-Efp/Esat ! eqn. 38
urel2=(up(n1)-up(n2))**2+(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2
A_rat=(radp(n1)**5*radp(n2)**5*sqrt(urel2)/(rad0**5)/((radp(n1)+radp(n2))**2)/(radp(n1)**3+ &
radp(n2)**3)/un0)**0.4_pr ! eqn. 27
Ac_rat=(5._pr*pi*rhop*urel2*radp(n1)**5*radp(n2)**5*(1._pr-nyp**2)/Ep/128._pr/((radp(n1)+ &
radp(n2))**2)/(radp(n1)**3+radp(n2)**3))**(0.4_pr) ! eqn. 34
Np_rat(n1)=A_rat*alpha_rat*c_rat(n1) ! eqn. 22
Np_rat(n2)=A_rat*alpha_rat*c_rat(n2) ! eqn. 22
mu=mu0*sqrt(Np_rat(n1)+Np_rat(n2)) ! eqn. 15
sigma=(sigma0-dq0_min)*(Np_rat(n2)-Np_rat(n1)) ! eqn. 19
call random_number(random)
random_normal=cos(2._pr*pi*random(1))*sqrt(-2._pr*log(random(2)+1.e-12_pr)) ! Box-Muller method
dqp=mu*random_normal+sigma ! ...(mu*random_normal+sigma)
q_elNext(n1)=max(min(q_elNext(n1)+dqp,qpmax),qp0)
q_elNext(n2)=max(min(q_elNext(n2)-dqp,qpmax),qp0)
term_prev(n1)=term_prev(n1) + c_rat_prev(n1)*Ac_rat*dqp/mu ! eqn. 30
term_prev(n2)=term_prev(n2) + c_rat_prev(n2)*Ac_rat*dqp/mu ! eqn. 30
write(filename,'(a,i3.3,a)') 'results/particlesPart_p',myid,'.dat'
open(access='append',unit=10,file=filename)
write(10,77) nt,t,q_elNext(n1),q_elNext(n2),dqp,c_rat(n1), c_rat(n2)
77 format(i7,5(x,es11.3e2))
close(10)
end
src/particles.f90
+ 23
− 0
View file @ 4846249a
......@@ -112,6 +112,10 @@
call allocateParticleArray(nppTot,fz_l)
call allocateParticleArray(nppTot,partn)
call allocateParticleArray(nppTot,cSpecies)
call allocateParticleArray(nppTot,term_prev)
call allocateParticleArray(nppTot,c_rat)
call allocateParticleArray(nppTot,c_rat_prev)
call allocateParticleArray(nppTot,Np_rat)
call allocateParticleArrayI(nppTot,ip)
call allocateParticleArrayI(nppTot,jp)
......@@ -119,6 +123,7 @@
call allocateParticleArrayI(nppTot,wcollnum)
call allocateParticleArrayI(nppTot,ppcollnum)
call allocateParticleArrayI(nppTot,nGlob)
call allocateParticleArrayL(nppTot, first_call_particle)
contains
......@@ -158,6 +163,24 @@
end
subroutine allocateParticleArrayL(nppTot, myvar)
use var
logical, allocatable, dimension(:) :: myvar, tmyvar
integer :: nppTot
if (nppTot .ge. 1) then
allocate(tmyvar(maxnp))
tmyvar(1:npp) = myvar(1:npp)
endif
if (allocated(myvar)) deallocate(myvar)
allocate(myvar(maxnp))
myvar = .true.
if (nppTot .ge. 1) then
myvar(1:npp) = tmyvar(1:npp)
endif
end
end
!####################################################################
......
src/particlesTransport.f90
+ 5
− 0
View file @ 4846249a
......@@ -200,6 +200,7 @@
if (qpmax.ne.qp0) then
if (Qmodel.eq.'r') call chargeParticleWallRandom(n,1)
if (Qmodel.eq.'c') call chargeParticleWallCondenser(n,1)
if (Qmodel.eq.'s') call chargeParticleWallSingle(n,1)
endif
endif
elseif (bcx.eq.'i') then
......@@ -227,6 +228,7 @@
if (qpmax.ne.qp0) then
if (Qmodel.eq.'r') call chargeParticleWallRandom(n,2)
if (Qmodel.eq.'c') call chargeParticleWallCondenser(n,2)
if (Qmodel.eq.'s') call chargeParticleWallSingle(n,2)
endif
endif
endif
......@@ -246,6 +248,7 @@
if (qpmax.ne.qp0) then
if (Qmodel.eq.'r') call chargeParticleWallRandom(n,3)
if (Qmodel.eq.'c') call chargeParticleWallCondenser(n,3)
if (Qmodel.eq.'s') call chargeParticleWallSingle(n,3)
endif
endif
endif
......@@ -375,6 +378,8 @@
if (q_el(n1).ne.q_el(n2)) call chargeParticleParticleCondenser(n1,n2)
elseif (Qmodel.eq.'p') then
call chargeParticleParticlePatch(n1,n2)
elseif (Qmodel.eq.'s') then
call chargeParticleParticleSingle(n1,n2)
endif
! fictitious contact point
xpc(n1)= xp(n1)+fracdt*up(n1)*dtNext
......
src/var.f90
+ 14
− 0
View file @ 4846249a
......@@ -26,6 +26,12 @@
urfp= 4.0_pr, &
C_s= 0.17_pr, & !< Smagorinsky constant
elementaryCharge= -1.60e-19_pr! elementary charge of electron
mu0=1.0_pr, & !< mu
rad0=1.0_pr, & !< particle radius of reference impact condition
un0=1.0_pr, & !< normal impact velocity of reference impact condition
Esat=1.0_pr, & !< electric field at saturation
dq0_min=1.0_pr, & !< the smallest measured impact charge
sigma0=1.0_pr !< sigma
real(kind=pr), parameter :: &
Ew=1.e11_pr, & !< duct Young's modulus (kg/s**2/m)
Ep=1.e8_pr, & !< particle Young's modulus (kg/s**2/m)
......@@ -198,4 +204,12 @@
phi_el, & !< electric potential (V)
Ex_el,Ey_el,Ez_el !< electric field strength (V/m)
! single particle charging model
logical, allocatable, dimension(:), save :: first_call_particle
real(kind=pr), allocatable, dimension(:) :: &
term_prev, &
c_rat, &
c_rat_prev, &
Np_rat
end module var
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment