Skip to content

Commit d210323

Browse files
RemiLeheRemiLehe
committed
Modify splitting of energy for ionized particles
1 parent 010fcbe commit d210323

File tree

1 file changed

+32
-88
lines changed

1 file changed

+32
-88
lines changed

‎Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H‎

Lines changed: 32 additions & 88 deletions
Original file line numberDiff line numberDiff line change
@@ -393,94 +393,38 @@ public:
393393
uz_p2 -= uCOM_z;
394394

395395
if (mask[i] == int(ScatteringProcessType::IONIZATION)) {
396-
// calculate kinetic energy of the collision (in eV)
397-
const amrex::ParticleReal E1 = (
398-
0.5_prt * m1 * (ux1*ux1 + uy1*uy1 + uz1*uz1) / PhysConst::q_e
399-
);
400-
const amrex::ParticleReal E2 = (
401-
0.5_prt * m2 * (ux_p2*ux_p2 + uy_p2*uy_p2 + uz_p2*uz_p2) / PhysConst::q_e
402-
);
403-
const amrex::ParticleReal E_coll = E1 + E2;
404-
405-
// subtract the energy cost for ionization
406-
const amrex::ParticleReal E_out = (E_coll - ionization_energy) * PhysConst::q_e;
407-
408-
// Momentum and energy division after the ionization event
409-
// is done as follows:
410-
// Three numbers are generated that satisfy the triangle
411-
// inequality. These numbers will be scaled to give the
412-
// momentum for each particle (satisfying the triangle
413-
// inequality ensures that the momentum vectors can be
414-
// arranged such that the net linear momentum is zero).
415-
// Product 2 is definitely an ion, so we first choose its
416-
// proportion to be n_2 = min(E2 / E_coll, 0.5 * R), where
417-
// R is a random number between 0 and 1.
418-
// The other two numbers (n_0 & n_1) are then found
419-
// by choosing a random point on the ellipse characterized
420-
// by the semi-major and semi-minor axes
421-
// a = (1 - n_0) / 2.0
422-
// b = 0.5 * sqrt(1 - 2 * n_0)
423-
// The numbers are found by randomly sampling an x value
424-
// between -a and a, and finding the corresponding y value
425-
// that falls on the ellipse: y^2 = b^2 - b^2/a^2 * x^2.
426-
// Then n_0 = sqrt(y^2 + (x - n_2/2)^2) and
427-
// n_1 = 1 - n_2 - n_0.
428-
// Next, we need to find the number C, such that if
429-
// p_i = C * n_i, we would have E_0 + E_1 + E_2 = E_out
430-
// where E_i = p_i^2 / (2 M_i), i.e.,
431-
// C^2 * \sum_i [n_i^2 / (2 M_i)] = E_out
432-
// After C is determined the momentum vectors are arranged
433-
// such that the net linear momentum is zero.
434-
435-
const amrex::ParticleReal n_2 = std::min(E2 / E_coll, 0.5_prt * amrex::Random(engine));
436-
437-
// find ellipse semi-major and minor axis
438-
const amrex::ParticleReal a = 0.5_prt * (1.0_prt - n_2);
439-
const amrex::ParticleReal b = 0.5_prt * std::sqrt(1.0_prt - 2.0_prt * n_2);
440-
441-
// sample random x value and calculate y
442-
const amrex::ParticleReal x = (2._prt * amrex::Random(engine) - 1.0_prt) * a;
443-
const amrex::ParticleReal y2 = b*b - b*b/(a*a) * x*x;
444-
const amrex::ParticleReal n_0 = std::sqrt(y2 + x*x - x*n_2 + 0.25_prt*n_2*n_2);
445-
const amrex::ParticleReal n_1 = 1.0_prt - n_0 - n_2;
446-
447-
// calculate the value of C
448-
const amrex::ParticleReal C = std::sqrt(E_out / (
449-
n_0*n_0 / (2.0_prt * m1) + n_1*n_1 / (2.0_prt * m_ioniz_product1) + n_2*n_2 / (2.0_prt * m_ioniz_product2)
450-
));
451-
452-
// Now that appropriate momenta are set for each outgoing species
453-
// the directions for the velocity vectors must be chosen such
454-
// that the net linear momentum in the current frame is 0.
455-
// This is achieved by arranging the momentum vectors in
456-
// a triangle and finding the required angles between the vectors.
457-
const amrex::ParticleReal cos_alpha = (n_0*n_0 + n_1*n_1 - n_2*n_2) / (2.0_prt * n_0 * n_1);
458-
const amrex::ParticleReal sin_alpha = std::sqrt(1.0_prt - cos_alpha*cos_alpha);
459-
const amrex::ParticleReal cos_gamma = (n_0*n_0 + n_2*n_2 - n_1*n_1) / (2.0_prt * n_0 * n_2);
460-
const amrex::ParticleReal sin_gamma = std::sqrt(1.0_prt - cos_gamma*cos_gamma);
461-
462-
// choose random theta and phi values (orientation of the triangle)
463-
const amrex::ParticleReal Theta = amrex::Random(engine) * 2.0_prt * MathConst::pi;
464-
const amrex::ParticleReal phi = amrex::Random(engine) * MathConst::pi;
465-
466-
const amrex::ParticleReal cos_theta = std::cos(Theta);
467-
const amrex::ParticleReal sin_theta = std::sin(Theta);
468-
const amrex::ParticleReal cos_phi = std::cos(phi);
469-
const amrex::ParticleReal sin_phi = std::sin(phi);
470-
471-
// calculate the velocity components for each particle
472-
ux1 = C * n_0 / m1 * cos_theta * cos_phi;
473-
uy1 = C * n_0 / m1 * cos_theta * sin_phi;
474-
uz1 = -C * n_0 / m1 * sin_theta;
475-
476-
ux_p1 = C * n_1 / m_ioniz_product1 * (-cos_alpha * cos_theta * cos_phi - sin_alpha * sin_phi);
477-
uy_p1 = C * n_1 / m_ioniz_product1 * (-cos_alpha * cos_theta * sin_phi + sin_alpha * cos_phi);
478-
uz_p1 = C * n_1 / m_ioniz_product1 * (cos_alpha * sin_theta);
479-
480-
ux_p2 = C * n_2 / m_ioniz_product2 * (-cos_gamma * cos_theta * cos_phi + sin_gamma * sin_phi);
481-
uy_p2 = C * n_2 / m_ioniz_product2 * (-cos_gamma * cos_theta * sin_phi - sin_gamma * cos_phi);
482-
uz_p2 = C * n_2 / m_ioniz_product2 * (cos_gamma * sin_theta);
483-
396+
// calculate kinetic energy of the collision
397+
const amrex::ParticleReal p_non_target_in = m1 * std::sqrt(ux1*ux1 + uy1*uy1 + uz1*uz1);
398+
const amrex::ParticleReal E_coll = 0.5_prt * p_non_target_in * (1.0_prt/m1 + 1.0_prt/m2);
399+
// subtract the energy cost for ionization (converted from eV to J)
400+
const amrex::ParticleReal E_out = E_coll - ionization_energy * PhysConst::q_e;
401+
402+
// Momentum and energy division after the ionization event is done as follows:
403+
// we assume that, in the center of mass frame, the two products (the ion and the electron)
404+
// have equal velocity in and that their velocity vector is aligned with that of the target particle
405+
406+
// With these assumptions, conservation of momentum in the center of mass frame allows to express the momenta
407+
// of the products, as a function of the outcoming momentum of the non-target particle:
408+
// p_product1 = - p_non_target_out * m_ioniz_product1 / m2
409+
// p_product2 = - p_non_target_out * m_ioniz_product2 / m2
410+
// (with m2 = m_ioniz_product1 + m_ioniz_product2, by conservation of mass)
411+
412+
// Amplitude of the momentum of the non-target particle
413+
const amrex::ParticleReal p_non_target_out = std::sqrt( 2.0_prt * E_out / (1.0_prt / m1 + 1.0_prt / m2 ) );
414+
415+
// Reduce the velocity of the non-target particle accordingly
416+
const amrex::ParticleReal reduce_factor = p_non_target_out / p_non_target_in;
417+
ux1 *= reduce_factor;
418+
uy1 *= reduce_factor;
419+
uz1 *= reduce_factor;
420+
// Obtain the other velocity vectors by conservation of momentum
421+
const amrex::ParticleReal mass_ratio = m1 / m2;
422+
ux_p1 = mass_ratio * m_ioniz_product1 * ux1;
423+
uy_p1 = mass_ratio * m_ioniz_product1 * uy1;
424+
uz_p1 = mass_ratio * m_ioniz_product1 * uz1;
425+
ux_p2 = mass_ratio * m_ioniz_product2 * ux1;
426+
uy_p2 = mass_ratio * m_ioniz_product2 * uy1;
427+
uz_p2 = mass_ratio * m_ioniz_product2 * uz1;
484428
}
485429
else {
486430
amrex::Abort("Unknown scattering process.");

0 commit comments

Comments
 (0)