## 2 Mathematical Model
### Governing Equations  
According to [PERSON]'s second law \(F=ma\), the control equation of ice crystal movement can be deduced as follows [20]:  
\[m_{ic}\frac{d^{2}\overset{\rightharpoonup}{r_{ic}}}{dt^{2}}=\overset{ \rightharpoonup}{P}+\overset{\rightharpoonup}{M_{a}}+\overset{ \rightharpoonup}{G}+\overset{\rightharpoonup}{B}+\overset{ \rightharpoonup}{D}\, \tag{1}\]  
where \(\overset{\rightharpoonup}{D}\) is the viscous drag force of air, \(\overset{\rightharpoonup}{P}\) is the force due to the difference in pressure gradient between front and back, \(\overset{\rightharpoonup}{M}\) is the apparent mass force, \(\overset{\rightharpoonup}{G}\) is gravity, \(\overset{\rightharpoonup}{B}\) is the unsteady force, \(m_{ic}\) is the mass of the ice crystal, \(d^{2}\overset{\rightharpoonup}{r_{ic}}/dt^{2}\) is the acceleration of motion, and the position of the ice crystal is \(\overset{\rightharpoonup}{r_{ic}}=(x,y)\). Since the density of an ice crystal is much greater than the density of air, the apparent mass force and the force due to the difference in pressure gradients between the front and the back of the ice crystal can be neglected. The unsteady force indicates that the ice crystal movement process is unsteady flow, and the steady drag coefficient is used in the calculation, so there is a certain error between them. In order to correct this error, an additional force \(\overset{\rightharpoonup}{B}\) is introduced, which is usually ignored in the calculation. Therefore, the control equation of ice crystal movement can be simplified to the following form:  
\[m_{w}\frac{d^{2}\overset{\rightharpoonup}{r_{w}}}{dt^{2}}=\overset{ \rightharpoonup}{D}+\overset{\rightharpoonup}{G}\, \tag{2}\]  
Gravity can be calculated using the following equation:  
\[\overset{\rightharpoonup}{G}=(\rho_{ic}-\rho_{a})v_{ic}\overset{ \rightharpoonup}{g}\, \tag{3}\]  
where \(\rho_{a}\) is the air density, \(\rho_{ic}\) is the ice crystal density, \(v_{ic}\) is the ice crystal volume, and \(\overset{\rightharpoonup}{g}\) is the acceleration of gravity. The viscous drag can be expressed using the following equation:  
\[D=C_{d}\times 0.5\rho_{a}\cdot U^{2}\cdot A_{ic}, \tag{4}\]  
where \(A_{ic}\) is the windward area of the ice crystal, \(A_{ic}=\pi r^{2}\). \(U\) is the vector difference between the ice crystal velocity and air velocity, \(U=\left|\overset{\rightharpoonup}{u}_{a}-\overset{\rightharpoonup}{u}_{ic}\right|\). \(C_{d}\) is the drag coefficient. Therefore, Equation (2) can also be expressed as:  
\[m_{ic}\overset{\rightharpoonup}{a}_{ic}=(\rho_{ic}-\rho_{a})v_{ic}\overset{ \rightharpoonup}{g}+0.5\times C_{d}\rho_{a}A_{ic}U^{2}\, \tag{5}\]  
In the calculation of drag coefficient \(C_{d}\), the drag coefficient and the relative Reynolds number are usually written as a function of drag: \(f=C_{d}\text{Re}/24\). When the surrounding flow satisfies Stokes' law, it can be expressed as:  
\[C_{d}\text{Re}/24=1, \tag{6}\]  
It is generally accepted that when relative Reynolds numbers \(\text{Re}<1\), the above equation can be used. When the surrounding flow no longer satisfies Stokes' law, it can be assumed that the value \(f\) only varies with Re, and the value can be determined through relevant experiments. In this paper, when calculating the ice crystal motion trajectory, it is considered that the ice crystal is spherical, and the spherical drag coefficient equation can be used. The common drag coefficient equations are as follows:  
\[C_{d}\text{Re}/24=1+0.197\text{Re}^{0.63}+2.6\times 10^{-4}\text{Re}^{1.38}\, \tag{7}\]  
\[C_{d}\text{Re}/24=1+0.15\text{Re}^{0.687}\, \tag{8}\]\[C_{d}\text{Re}/24=1+1/6\text{Re}^{2/3}\, \tag{9}\]  
where the relative Reynolds number can be expressed as:  
\[\text{Re}=\rho_{a}\left|\overset{\rightarrow}{u}_{a}-\overset{\rightarrow}{u}_{ ic}\right|MVD/\mu_{a}, \tag{10}\]  
This equation is decomposed in each direction and can be written as:  
\[