diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 567bd68..659cbd7 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -56,6 +56,15 @@ MODULE moduleMesh1DCart END SUBROUTINE transparent + MODULE SUBROUTINE wallTemperature(edge, part) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + + END SUBROUTINE wallTemperature + END INTERFACE TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DCart @@ -189,6 +198,9 @@ MODULE moduleMesh1DCart TYPE IS(boundaryTransparent) self%fBoundary(s)%apply => transparent + TYPE IS(boundaryWallTemperature) + self%fBoundary(s)%apply => wallTemperature + CLASS DEFAULT CALL criticalError("Boundary type not defined in this geometry", 'initEdge1DCart') diff --git a/src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90 b/src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90 index 7bf0f7c..443c7cc 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90 @@ -63,4 +63,38 @@ SUBMODULE (moduleMesh1DCart) moduleMesh1DCartBoundary END SUBROUTINE transparent + SUBROUTINE wallTemperature(edge, part) + USE moduleSpecies + USE moduleBoundary + USE moduleConstParam, ONLY: PI + IMPLICIT NONE + + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + REAL(8):: x, y + + !Modifies particle velocity according to wall temperature + SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj) + TYPE IS(boundaryWallTemperature) + x = 0.D0 + DO WHILE (x == 0.D0) + CALL RANDOM_NUMBER(x) + + END DO + + CALL RANDOM_NUMBER(y) + + part%v(1) = part%v(1) + bound%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + + END SELECT + + SELECT TYPE(edge) + TYPE IS(meshEdge1DCart) + part%v(1) = -part%v(1) + part%r(1) = 2.D0*edge%x - part%r(1) + + END SELECT + + END SUBROUTINE wallTemperature + END SUBMODULE moduleMesh1DCartBoundary diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 35240c0..9693319 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -56,6 +56,15 @@ MODULE moduleMesh1DRad END SUBROUTINE transparent + MODULE SUBROUTINE wallTemperature(edge, part) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + + END SUBROUTINE wallTemperature + END INTERFACE TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DRad @@ -190,6 +199,9 @@ MODULE moduleMesh1DRad TYPE IS(boundaryTransparent) self%fBoundary(s)%apply => transparent + TYPE IS(boundaryWallTemperature) + self%fBoundary(s)%apply => wallTemperature + CLASS DEFAULT CALL criticalError("Boundary type not defined in this geometry", 'initEdge1DRad') diff --git a/src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90 b/src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90 index 9dbc32f..a9cdca8 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90 @@ -65,4 +65,38 @@ SUBMODULE (moduleMesh1DRad) moduleMesh1DRadBoundary END SUBROUTINE transparent + SUBROUTINE wallTemperature(edge, part) + USE moduleSpecies + USE moduleBoundary + USE moduleConstParam, ONLY: PI + IMPLICIT NONE + + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + REAL(8):: x, y + + !Modifies particle velocity according to wall temperature + SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj) + TYPE IS(boundaryWallTemperature) + x = 0.D0 + DO WHILE (x == 0.D0) + CALL RANDOM_NUMBER(x) + + END DO + + CALL RANDOM_NUMBER(y) + + part%v(1) = part%v(1) + bound%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + + END SELECT + + SELECT TYPE(edge) + TYPE IS(meshEdge1DRad) + part%v(1) = -part%v(1) + part%r(1) = 2.D0*edge%r - part%r(1) + + END SELECT + + END SUBROUTINE wallTemperature + END SUBMODULE moduleMesh1DRadBoundary diff --git a/src/modules/mesh/2DCyl/moduleMeshCyl.f90 b/src/modules/mesh/2DCyl/moduleMeshCyl.f90 index 6adc792..c4f2b55 100644 --- a/src/modules/mesh/2DCyl/moduleMeshCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMeshCyl.f90 @@ -55,6 +55,15 @@ MODULE moduleMeshCyl END SUBROUTINE absorption + MODULE SUBROUTINE wallTemperature(edge, part) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + + END SUBROUTINE wallTemperature + MODULE SUBROUTINE transparent(edge, part) USE moduleSpecies IMPLICIT NONE @@ -247,6 +256,9 @@ MODULE moduleMeshCyl TYPE IS(boundaryTransparent) self%fBoundary(s)%apply => transparent + TYPE IS(boundaryWallTemperature) + self%fBoundary(s)%apply => wallTemperature + TYPE IS(boundaryAxis) self%fBoundary(s)%apply => symmetryAxis diff --git a/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 b/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 index 825699d..556abf7 100644 --- a/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 +++ b/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 @@ -100,6 +100,66 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary END SUBROUTINE transparent + !Wall with temperature + SUBROUTINE wallTemperature(edge, part) + USE moduleSpecies + USE moduleBoundary + USE moduleConstParam, ONLY: PI + IMPLICIT NONE + + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2) + INTEGER:: i + REAL(8):: x, y + + !Modifies particle velocity according to wall temperature + SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj) + TYPE IS(boundaryWallTemperature) + DO i = 1, 3 + x = 0.D0 + DO WHILE (x == 0.D0) + CALL RANDOM_NUMBER(x) + + END DO + + CALL RANDOM_NUMBER(y) + + part%v(i) = part%v(i) + bound%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + + END DO + + END SELECT + + !Reflects particle in the edge + SELECT TYPE(edge) + TYPE IS(meshEdgeCyl) + edgeNorm = DSQRT((edge%r(2)-edge%r(1))**2 + (edge%z(2)-edge%z(1))**2) + cosT = (edge%z(2)-edge%z(1))/edgeNorm + sinT = DSQRT(1-cosT**2) + + rp(1) = part%r(1) - edge%z(1); + rp(2) = part%r(2) - edge%r(1); + + rpp(1) = cosT*rp(1) - sinT*rp(2) + rpp(2) = sinT*rp(1) + cosT*rp(2) + rpp(2) = -rpp(2) + + vpp(1) = cosT*part%v(1) - sinT*part%v(2) + vpp(2) = sinT*part%v(1) + cosT*part%v(2) + vpp(2) = -vpp(2) + + part%r(1) = cosT*rpp(1) + sinT*rpp(2) + edge%z(1); + part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%r(1); + part%v(1) = cosT*vpp(1) + sinT*vpp(2) + part%v(2) = -sinT*vpp(1) + cosT*vpp(2) + + END SELECT + + part%n_in = .TRUE. + + END SUBROUTINE wallTemperature + !Symmetry axis. Dummy function SUBROUTINE symmetryAxis(edge, part) USE moduleSpecies diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index dec8e83..0f54ee9 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -24,6 +24,14 @@ MODULE moduleBoundary END TYPE boundaryTransparent + !Transparent boundary + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature + !Thermal velocity of the wall: square root(Wall temperature X specific heat) + REAL(8):: vTh + CONTAINS + + END TYPE boundaryWallTemperature + !Symmetry axis TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis CONTAINS @@ -65,4 +73,17 @@ MODULE moduleBoundary END FUNCTION getBoundaryId + SUBROUTINE initWallTemperature(boundary, T, c) + USE moduleRefParam + IMPLICIT NONE + + CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary + REAL(8), INTENT(in):: T, c !Wall temperature and specific heat + REAL(8):: vTh + + vTh = DSQRT(c * T) / v_ref + boundary = boundaryWallTemperature(vTh = vTh) + + END SUBROUTINE initWallTemperature + END MODULE moduleBoundary diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index ab770a3..fde68b8 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -343,6 +343,7 @@ MODULE moduleInput INTEGER:: i, s CHARACTER(2):: istring, sString CHARACTER(:), ALLOCATABLE:: object, bType + REAL(8):: Tw, cw !Wall temperature and specific heat LOGICAL:: found INTEGER:: nTypes @@ -372,6 +373,14 @@ MODULE moduleInput CASE('transparent') ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj) + CASE('wallTemperature') + CALL config%get(object // '.temperature', Tw, found) + IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary') + CALL config%get(object // '.specificHeat', cw, found) + IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary') + + CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw) + CASE('axis') ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj)