There are three basic ingredients in the hard-sphere gas simulation that will be dealt with in turn: 1) the propagation of a sphere between bounces, 2) the bounce off a wall, and 3) the bounce off each other.<\/p>\n\n\n\n
Free propagation of each hard sphere between wall and particle interactions is done via simple rectilinear propagation. Assuming that the initial position and velocity of the sphere is ${\\vec r}_0$ and ${\\vec v}$, then its subsequent location is<\/p>\n\n\n\n
\\[ {\\vec r}(t) = {\\vec r}_0 + {\\vec v} \\cdot (t-t_0) \\; . \\]<\/p>\n\n\n\n
The velocity remains unchanged until collision with a wall or another hard sphere.<\/p>\n\n\n\n
When the atoms reach a wall, they bounce elastically with their normal component of velocity being negated.\u00a0 For example, in the example pictured here,<\/p>\n\n\n\n we see a particle approaching a wall in the $y-z$ plane.\u00a0 After the bounce, the $x$-component of the velocity has changed but the transverse components in $y$ and $z$ have not:<\/p>\n\n\n\n \\[\u00a0 {\\vec v} = v_x {\\hat x} + v_y {\\hat y} + v_z {\\hat z} \\; \\rightarrow \\; {\\vec v} = \\, – v_x {\\hat x} + v_y {\\hat y} + v_z {\\hat z} \\; .\\]<\/p>\n\n\n\n Since time in the simulation elapses in discrete steps, collision with the wall is checked for after each time step.\u00a0 If a component of the particle’s position puts it beyond the confining box (i.e., it passed through a wall at $x = L$ so that its current $x$-component of position is $x = L + \\delta_x$ ) then that component of the velocity is negated as above and the distance it had traveled beyond the wall is increment to the appropriate component in the opposite direction (i.e., $x = L – \\delta_x$).<\/p>\n\n\n\n The collision between hard spheres is the most complicated part of this simulation.\u00a0 Like the case with the wall, it will often be the case that at the end of a time step, two particles have their hard body radii $R_1$ and $R_2$ overlapping.<\/p>\n\n\n\n Once this overlapping end condition is detected, the next step is to determine the time $\\delta t$ to back up such that the spheres are just touching. The algorithm to determine the time step is most easily understood in a frame in which one of the particles, say $m_2$ in blue, is at the center of the coordinate system and the relative velocity vector ${\\vec v} = {\\vec v}_1 – {\\vec v}_2$ is along the horizontal. Letting ${\\vec r} = {\\vec r}_1 – {\\vec r}_2$ be the relative position of $m_1$ relative to $m_2$ then the unknown quantity is the vector ${\\vec d}$ between the current relative position and the relative position ${\\vec r}^{\\, *}$ when the two spheres are just touching.<\/p>\n\n\n\n Since ${\\vec r}^* + {\\vec d} = {\\vec r}$ and we know that the length of ${\\vec r}^*$ is $R_1 + R_2 = 2 R$ (i.e., $R$ is the average radius) then a quadratic formula can be developed from<\/p>\n\n\n\n \\[ |{\\vec r}^{\\, *}|^2 = |{\\vec r} – {\\vec d}|^2 \\; . \\]<\/p>\n\n\n\n After some straightforward algebra, we arrive at<\/p>\n\n\n\n \\[ d^2 – 2 r \\cos (\\theta) d + (r^2 – 4R^2) = 0 \\; , \\]<\/p>\n\n\n\n where $d = |{\\vec d}|$, $r = |{\\vec r}|$ and $\\cos (\\theta) = {\\hat r} \\cdot {\\hat V}$, since the relative velocity ${\\vec V}$ is parallel to ${\\vec d}$.<\/p>\n\n\n\n Solving for $d$ yields<\/p>\n\n\n\n \\[ d = r \\cos (\\theta) – \\sqrt{ 4R^2 + r^2 ( \\cos(\\theta)^2 – 1 ) } \\; , \\]<\/p>\n\n\n\n where we picked the negative radical because it guarantees $\\delta t$ to be negative signifying that the contact time is in the past as we know it must be.<\/p>\n\n\n\n The final ingredient is to determine how the velocities of the two particles should alter due to their collision. The direction is straightforward as the component of each particle normal to the tangent plane containing their contact must change sign but the magnitude requires a bit more care.<\/p>\n\n\n\n The analysis begins first by calculating the center-of-mass <\/p>\n\n\n\n via the usual formula<\/p>\n\n\n\n \\[ m_1 {\\vec r}_1 + m_2 {\\vec r}_2 = M {\\vec R} \\; . \\]<\/p>\n\n\n\n Taking a time derivative of this definition<\/p>\n\n\n\n \\[ {\\vec p}_1 + {\\vec p}_2 = {\\vec P}_{tot} = m_1 {\\vec v}_1 + m_2 {\\vec v}_2 = M {\\vec U} \\; \\]<\/p>\n\n\n\n leads to a useful relationship between the center-of-mass kinematics, defined in terms of the total mass $M$ and the velocity of the center-of-mass ${\\vec U}$<\/p>\n\n\n\n \\[ {\\vec U} = \\frac{ {\\vec P}_{tot} }{M} \\; .\\]<\/p>\n\n\n\n The next step is finding the positions, ${\\vec \\rho}_i$, of both masses relative to the center-of-mass using ${\\vec R} + {\\vec \\rho}_i = {\\vec r}_i $ to yield<\/p>\n\n\n\n \\[ {\\vec \\rho}_i = {\\vec r}_i – {\\vec R} \\; . \\]<\/p>\n\n\n\n Multiplying by $m_i$, taking the time derivative, and defining the momentum of each mass relative to the center-of-mass as ${\\vec \\pi}_i = m_i \\frac{d}{dt} {\\vec \\rho}_i$ gives the expression<\/p>\n\n\n\n \\[ {\\vec \\pi}_i = {\\vec p}_i – m_i {\\vec U} = {\\vec p}_i – m_i \\frac{{\\vec P}_{tot}}{M} \\; , \\]<\/p>\n\n\n\n where the last equality follows from the earlier relationship between the total momentum and the velocity of the center-of-mass.<\/p>\n\n\n\n The physical interpretation of the ${\\vec \\pi}_i$ is most clearly revealed by summing them<\/p>\n\n\n\n \\[ {\\vec \\pi}_1 + {\\vec \\pi}_2 = {\\vec p}_1 – m_1 \\frac{{\\vec P}_{tot}}{M} + {\\vec p}_2 – m_2 \\frac{{\\vec P}_{tot}}{M} \\\\ = ({\\vec p}_1 + {\\vec p}_2) – (m_1 + m_2) \\frac{{\\vec P}_{tot}}{M} = 0 \\; .\\]<\/p>\n\n\n\n This relation and the one before supply the physical interpretation for the ${\\vec \\pi}$\u2019s: they are the momenta of each particle in a coordinate frame comoving with the center of mass (i.e., relative to the center-of-mass). The Galilean transform that affects this result is said to have transformed the momenta (or velocity) into the center-of-momentum frame.<\/p>\n\n\n\n The application of the on-contact bounce is trivial in this frame. The total change in momentum is done only in the direction parallel to the line joining the two masses (i.e., normal to the tangent plane containing the contact)<\/p>\n\n\n\n \\[ \\Delta {\\vec p} = ({\\vec \\pi}_i \\cdot {\\hat \\rho}_i ) {\\hat \\rho}_i \\; \\]<\/p>\n\n\n\n and must be equal and opposite for each particle.<\/p>\n\n\n\n Each of the ${\\vec \\pi}_i$ changes to its after bounce value ${\\vec \\pi}’_i$ according to<\/p>\n\n\n\n \\[ {\\vec \\pi}’_i = {\\vec \\pi}_i – 2 \\Delta {\\vec p} \\; .\\]<\/p>\n\n\n\n The change in velocity is then found by first reversing the transformation to get ${\\vec p}’_i$ via<\/p>\n\n\n\n \\[ {\\vec p}’_i = {\\vec \\pi}’_i + m_i \\frac{{\\vec P}_{tot}}{M} \\; , \\]<\/p>\n\n\n\n from which immediately follows<\/p>\n\n\n\n \\[ {\\vec v}’_i = {\\vec p}’_i \/ m_i \\; . \\]<\/p>\n\n\n\n Next month\u2019s column covers the details of how these theoretical results are translated into a simulation and explores the results.<\/p>\n","protected":false},"excerpt":{"rendered":"![\"\"](\"http:\/\/underthehood.blogwyrm.com\/wp-content\/uploads\/2023\/06\/UTH_06Jun_2023_Wall_Bounce.png\")
<\/figure>\n\n\n\n![\"\"](\"http:\/\/underthehood.blogwyrm.com\/wp-content\/uploads\/2023\/06\/UTH_06Jun_2023_overlapping_spheres.png\")
<\/figure><\/div>\n\n\n\n![\"\"](\"http:\/\/underthehood.blogwyrm.com\/wp-content\/uploads\/2023\/06\/UTH_06Jun_2023_deltat_geometry.png\")
<\/figure><\/div>\n\n\n\n![\"\"](\"http:\/\/underthehood.blogwyrm.com\/wp-content\/uploads\/2023\/06\/UTH_06Jun_2023_center_of_mass_geometry.png\")
<\/figure><\/div>\n\n\n\n