From fea4f21e177dc12b1e8b7571844653ee57d59d46 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 28 Sep 2024 20:55:03 +0200 Subject: [PATCH] Partial case with Poisson equation It reproduces Diko's peak! --- plasmaExpansion.f90 | 41 +++++++++++++++++++++-------------------- 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90 index 4e906d5..1877a88 100644 --- a/plasmaExpansion.f90 +++ b/plasmaExpansion.f90 @@ -21,7 +21,7 @@ module referenceValues implicit none real(dp):: L_ref, t_ref, n_ref, u_ref, Temp_ref ! Reference values - real(dp):: phi_ref, p_ref ! Reference values + real(dp):: phi_ref ! Reference values end module referenceValues @@ -291,7 +291,8 @@ program plasmaExpansion use omp_lib implicit none - real(dp), parameter:: m_i = 1.9712258e-25_dp ! Tin aton mass in kg + real(dp), parameter:: m_i = 1.9712258e-25_dp ! Tin atom mass in kg + real(dp), parameter:: m_e = 9.1093837e-31_dp ! Electron mass in kg real(dp), parameter:: gam = 1.0_dp ! Adiabatic coefficient real(dp):: r0, rf @@ -330,7 +331,7 @@ program plasmaExpansion real(dp), allocatable, dimension(:):: phi, phi_old, E real(dp):: phiConv - real(dp):: phi0, phiF, JF + real(dp):: phi0, phiF integer:: k real(dp), allocatable, dimension(:):: fCum_i @@ -347,20 +348,19 @@ program plasmaExpansion u_ref = sqrt(kb * Temp_ref / m_i) L_ref = u_ref * t_ref phi_ref = kb * Temp_ref / qe - p_ref = n_ref * kb * Temp_ref ! Set position to calculate cumulative sum of f (non-dimensional units) rCum = 5.0e-3 / L_ref ! Set input parameters (remember these have to be in non-dimensional units) Temp0 = 60.0_dp * eV_to_K / Temp_ref - TempF = 60.0_dp * eV_to_K / Temp_ref + TempF = 5.0_dp * eV_to_K / Temp_ref n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref c_s = sqrt(T_to_Z(Temp0) * gam * Temp0) - u_bc0 = sqrt(Temp0) + u_bc0 = c_s!sqrt(Temp0) u_bcF = sqrt(TempF) n_bc0 = n_ecr / T_to_Z(Temp0) - n_bcF = n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0) + n_bcF = n_ecr*1.0e-1 / T_to_Z(Temp0) ! Set domain boundaries (non-dimensional units) r0 = 200.0e-6_dp / L_ref @@ -399,8 +399,8 @@ program plasmaExpansion nt = nint((tf - t0) / dt) dt = (tf - t0) / float(nt) - t_bc0 = nint( 20.0e-9_dp / t_ref / dt) - t_bcF = nint( 25.0e-9_dp / t_ref / dt) + t_bc0 = nint(100.0e-9_dp / t_ref / dt) + t_bcF = nint(105.0e-9_dp / t_ref / dt) everyOutput = nint(1.0e-9_dp/t_ref/dt) if (everyOutput == 0) then @@ -474,7 +474,6 @@ program plasmaExpansion ! phi(1) = phi0 ! Dirichlet phi0 = phi(1) ! Neumann phiF = 0.0_dp ! Dirichlet - JF = 0.0_dp ! Current at the edge for floating potential allocate(f0(j0:nv)) f0 = 0.0_dp @@ -547,7 +546,7 @@ program plasmaExpansion !$omp end parallel do ! Solve Poission (maximum number of iterations, break if convergence is reached before) - do k = 1, 100 + do k = 1, 500 ! Store previous value phi_old = phi @@ -579,23 +578,25 @@ program plasmaExpansion ! Iterate system call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info) phi = phi_old + Res - phi(1) = phi(2) ! Neumann + phi(1) = phi(2) ! Ensures smooth transition in Neumann boundary phi0 = phi(1) ! Neumann - ! phi(nr) = phi(nr-5) ! Dirichlet quasi-floating - - ! Calculate updated current - JF = Zave(nr)*n_i(nr)*u_i(nr) - n_e(nr)*sqrt(qe*T_i(1)*Temp_ref/9.1e-31_dp)/u_ref - - ! Iterate floating potential - phiF = phi(nr-5)! + 1.0e-2_dp*JF ! Check if the solution has converged phiConv = maxval(abs(Res),1) - if (phiConv < 1.0e-2_dp) then + if (phiConv < 1.0e-1_dp) then exit end if + ! Calculate new potential to ensure 0 current at the edge + if (n_i(nr) > 1.0e-10_dp) then + phiF = phi0 + T_i(1) * log((2.0_dp*sqrt(pi)*Zave(nr)*n_i(nr)*u_i(nr)) / (Zave(1)*n_i(1)*sqrt(m_i*T_i(1)/m_e))) + + else + phiF = phi(nr-5) + + end if + end do ! Calculate electric field