From 12e61731dfe96733a027235127fa423035e1823e Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Tue, 23 Mar 2021 16:43:11 +0100 Subject: [PATCH] Ionization boundary ready to testing. Fixed an issue in which some particles in the corner were interacting with the axis boundary. Now the axis acts as a reflective boundary in case a particle is wrongly assigned to it. --- src/modules/mesh/moduleMeshBoundary.f90 | 29 +++++++++++++++++-------- src/modules/moduleBoundary.f90 | 15 ++++++++----- src/modules/moduleInput.f90 | 9 +++++--- 3 files changed, 36 insertions(+), 17 deletions(-) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 7f5fb88..7553733 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -111,6 +111,7 @@ MODULE moduleMeshBoundary USE moduleSpecies USE moduleMesh USE moduleRefParam + USE moduleRandom IMPLICIT NONE CLASS(meshEdge), INTENT(inout):: edge @@ -118,6 +119,7 @@ MODULE moduleMeshBoundary REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass REAL(8):: ionizationRate INTEGER:: ionizationPair, p + REAL(8):: v0(1:3) !random velocity of neutral TYPE(particle), POINTER:: newElectron TYPE(particle), POINTER:: newIon @@ -128,21 +130,26 @@ MODULE moduleMeshBoundary eRel = mRel*vRel**2*5.D-1 IF (eRel > bound%eThreshold) THEN - !TODO: Check units - ionizationRate = bound%n0*bound%crossSection%get(eRel) + ionizationRate = part%weight*bound%n0*bound%crossSection%get(eRel)*vRel - ionizationPair = INT(ionizationRate*tauMin*ti_ref/species(bound%sp)%obj%weight) + !Rounds the number of particles up + ionizationPair = NINT(ionizationRate*bound%effectiveTime/species(bound%sp)%obj%weight) !Create the new pair of particles DO p = 1, ionizationPair ALLOCATE(newElectron) ALLOCATE(newIon) + !Assign random velocity to the neutral + v0(1) = bound%v0(1) + bound%vTh*randomMaxwellian() + v0(2) = bound%v0(2) + bound%vTh*randomMaxwellian() + v0(3) = bound%v0(3) + bound%vTh*randomMaxwellian() + newElectron%sp = part%sp newIon%sp = bound%sp - newElectron%v = 10.D0*bound%v0 + (1.D0 + bound%deltaV/NORM2(bound%v0)) - newIon%v = bound%v0 + newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0)) + newIon%v = v0 newElectron%r = edge%randPos() newIon%r = newElectron%r @@ -174,16 +181,18 @@ MODULE moduleMeshBoundary END DO - !Removes ionizing electron - part%n_in = .FALSE. - END IF END SELECT + !Removes ionizing electron regardless the number of pair created + part%n_in = .FALSE. + END SUBROUTINE ionization - !Symmetry axis. Dummy function + !Symmetry axis. Reflects particles. + !Although this function should never be called, it is set as a reflective boundary + !to properly deal with possible particles reaching a corner and selecting this boundary. SUBROUTINE symmetryAxis(edge, part) USE moduleSpecies IMPLICIT NONE @@ -191,6 +200,8 @@ MODULE moduleMeshBoundary CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part + CALL reflection(edge, part) + END SUBROUTINE symmetryAxis !Points the boundary function to specific type diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 66a98df..320bdad 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -35,9 +35,10 @@ MODULE moduleBoundary !Ionization boundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization - REAL(8):: m0, n0, v0(1:3), T0 !Properties of background neutrals. + REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals. INTEGER:: sp !Ion species TYPE(table1D):: crossSection + REAL(8):: effectiveTime REAL(8):: eThreshold REAL(8):: deltaV @@ -99,7 +100,7 @@ MODULE moduleBoundary END SUBROUTINE initWallTemperature - SUBROUTINE initIonization(boundary, m0, n0, v0, T0, speciesID, crossSection, eThreshold) + SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, speciesID, effTime, crossSection, eThreshold) USE moduleRefParam USE moduleSpecies USE moduleCaseParam @@ -107,8 +108,10 @@ MODULE moduleBoundary IMPLICIT NONE CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary - REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 + REAL(8), INTENT(in):: me !Electron mass + REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties INTEGER:: speciesID + REAL(8):: effTime CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection REAL(8), INTENT(in):: eThreshold @@ -117,13 +120,15 @@ MODULE moduleBoundary SELECT TYPE(boundary) TYPE IS(boundaryIonization) boundary%m0 = m0 / m_ref - boundary%n0 = n0 + boundary%n0 = n0 * Vol_ref boundary%v0 = v0 / v_ref - boundary%T0 = T0 / T_ref + boundary%vTh = DSQRT(kb*T0/m0)/v_ref boundary%sp = speciesID + boundary%effectiveTime = effTime / ti_ref CALL boundary%crossSection%init(crossSection) CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2) + boundary%deltaV = DSQRT(boundary%eThreshold/me) END SELECT diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 96debdd..b7d8bc4 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -600,9 +600,9 @@ MODULE moduleInput CHARACTER(:), ALLOCATABLE:: object, bType REAL(8):: Tw, cw !Wall temperature and specific heat !Neutral Properties - REAL(8):: m0, n0 + REAL(8):: m0, n0, T0 REAL(8), DIMENSION(:), ALLOCATABLE:: v0 - REAL(8):: T0 + REAL(8):: effTime REAL(8):: eThreshold !Energy threshold INTEGER:: speciesID CHARACTER(:), ALLOCATABLE:: speciesName, crossSection @@ -647,12 +647,15 @@ MODULE moduleInput IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary') CALL config%get(object // '.neutral.temperature', T0, found) IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary') + CALL config%get(object // '.effectiveTime', effTime, found) + IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for neutrals in ionization", 'readBoundary') CALL config%get(object // '.energyThreshold', eThreshold, found) IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary') CALL config%get(object // '.crossSection', crossSection, found) IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary') - CALL initIonization(boundary(i)%bTypes(s)%obj, m0, n0, v0, T0, speciesID, crossSection, eThreshold) + CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, & + speciesID, effTime, crossSection, eThreshold) CASE('wallTemperature') CALL config%get(object // '.temperature', Tw, found)