Texts on physics, maths and programming (Несколько текстов для зачёта), страница 30

Fig. 2. Schematic showing the mechanics of the proposed technology.

3. Mathematical model

In this research fractionation of particles is based on density and compressibility differences of fluid and particles rather than on particle size. Employing the above basic principles of physics of particles in an acoustic field, a mathematical model is developed to calculate trajectories of deflected particles subjected to acoustic standing waves. Table 1 gives the properties of the particles and the fluid that are used in this research.

Table 1.

Properties of particles and suspending medium

Description	Solid (SiC)	Medium (DI water)
Density (ρ) (kg/m3)	3217	1000
Frequency of sound in medium (f) (kHz)	–	333
Viscosity of medium (μ) (N s/μm2)	–	9.98E−16
Acoustic energy in medium (J/m3)	–	133
Power in medium (W/m3)	–	56 000
Quality factor (Q) of chamber	–	5000

Rearranging Eqs. Figs. (3) and (5) yields the following equation:

(6)

Simplifying Eq. (6) with values in Table 1 yields:

x″+cx′-ksinβx=0,

(7)

where

x′=ν,

(8)

where, c = 1.42E6, k = E8, and β = 2.78E−3 are constants representing the physical parameters of Eq. (7). This equation will be constantly used during the mathematical derivation. The notation (′) indicates the derivative, d/dt, and x is the position of the migrating particle in the x-direction between a transducer and a reflector separated by one half wavelength (=λ/2) of the resonant sound at the given frequency. The parameters c, k, and β are all positive constants having the following orders (O) of magnitude: c = O(10⁶), k = O(10⁸), β = O(10⁻³).

A study of the behavior of the solutions is discussed with an explanation of the available solution techniques. Solutions of Eq. (7) will be used as a basis for concluding some results during the derivation. The above equation is extremely stiff, so most numerical solution methods—even stiff equation solvers—provide little useful information.

Note that Eq. (7) has the form somewhat like a damped nonlinear spring (or pendulum) equation. Several publications cited in the literature assumed instantaneous viscous relaxation where the inertial term dv/dt or x″ was neglected. This type of singular perturbation approximation; namely:

cx′-ksinβx=0

(9)

has been used to approximate the solution for Eq. (7). It will be shown in this paper that although this approximation produces results that are qualitatively correct, quantitative errors are incurred that can be significant for some applications.

3.1. Preliminary analysis

There is no closed form solution for Eq. (7), so one has to rely on qualitative and approximate analysis techniques that will be discussed in this section. Useful information on the nature of the solutions of Eq. (7) can be obtained by recasting it as a two-dimensional system and performing a careful phase plane analysis. To do this, first set x′ = ν and then rewrite (7) as the following system in the x, ν-plane:

ν′=ksinβx-cν = (x′,ν′)=F(x,ν).

(10)

Eq. (10) in matrix form, with f₁ and f₂ as function of x′ and v′ respectively is:

(11)

The fixed (critical or stationary) points of Eq. (11) are those points in the x, ν-plane satisfying function F(x,ν)=0 ν=0, ksinβx-cν=0:

(12)

Hence, these points are (mπ/β, 0), where m is an integer. To find the local behavior of the solutions to Eq. (11) near the fixed points, the derivative of F, F′, is computed at these fixed points:

(13)

Consequently, when , and when m is even:

(14)

and when m is odd:

(15)

To compute the spectral properties of Eqs. Figs. (14) and (15), one needs to calculate the eigenvalues of the derivative matrices. The eigenvalues (μ) of Eq. (14) are the roots of the equation:

(16)

Hence, there is one negative eigenvalue ( ) and one positive eigenvalue ( ) given, respectively, by

(17)

(18)

where the approximate formulas follow from Eq. (7). The eigenspaces (E) associated to these eigenvalues are calculated as follows:

(19)

(20)

Similarly, the eigenvalues of Eq. (15) are the roots of the equation:

(21)

Both of these roots are negative. One of these has a large magnitude and the other has small absolute value; these are given as:

(22)

(23)

The corresponding eigenspaces (E) are:

(24)

(25)

It is concluded that the fixed points of the form (2lπ/β, 0) are all saddle points having local phase plane behavior as shown in Fig. 3, where l is any integer 0. Note that the angle θ between eigenvector (w^u) and x axis is .

Fig. 3. Phase plane near (2lπ/β, 0).

On the other hand, the fixed points of the form ((2l + 1)π/β, 0) are all sinks with local phase plane behavior shown below in Fig. 4. Fig. 5 shows the entire phase plane of Eq. (11) obtained by collecting the above results.

Fig. 4. Phase plane near ((2l + 1)π/β, 0).

Fig. 5. Phase portrait of Eq. (11).

Useful information on the solutions of Eq. (7) (and equivalently (11)) can be extracted from the above analysis. The relevant portion of the phase plane for our research is the region R = {(x, v):0   x    8π/β}. Clearly any trajectory θ_t(x₀, v₀) = (x(t), v(t)) starting at (x₀, v₀)   R, except for points on the stable manifold, has the property:

(26)

Moreover,

(27)

In fact, all of these trajectories converge very rapidly to the separatrix G (which is just the portion of the unstable manifold for the origin, , that is contained in R). This separatrix appears to be quite close to the graph of v = k/c sinβx, and so one can infer the following result that shall be proved in the sequel:

(28)

One can also see from the above analysis that the often-used singular perturbation approximation Eq. (8) is not really adequate to represent the solutions of Eq. (7) ( (11)). This is because it is just a single curve in the phase plane region R beginning at origin, O = (0, 0), and terminating at (π/β, 0); albeit a curve to which (almost) all of the solutions (trajectories, orbits) of Eq. (7) converge to rather rapidly as t → ∞(x → π/β).

Since Eq. (7) ( (11)) has no closed form solution, it is natural to try to compute numerical approximate solutions using a standard ODE solver such as the Runge–Kutta method. However, an attempt to apply the Runge–Kutta method to Eq. (7) directly leads to rather surprising results that are completely unsatisfactory. To explain what happens, first note that the system form of Eq. (7), namely (11), is such that there is an enormous discrepancy between the magnitude of the derivative of the first component f₁: = v and the corresponding magnitudes of the second component f₂:=ksinβx-cv; in fact, the ratios of the second set of magnitudes to the first is at least O(10⁸). Such systems of differential equations are called stiff. In the presence of such enormous magnitudes, one might expect to encounter some problems with conventional ordinary differential equation (ODE) solvers. More specifically, in order to obtain good accuracy using the Runge–Kutta method one should take a time step of size (h), no more than about h = 0.1 in order to benefit from the global approximation error of O(h⁴) of the Runge–Kutta method. If one uses this increment, it is found that there is a little movement in the trajectory after several hundred steps, so if one starts in R near x = 0, very little progress is made in moving across R to the point (π/β, 0) even after hundreds of steps. Furthermore, if one continues this for thousands of steps, round-off errors accumulate to the point where the approximate numerical solution starts behaving in ways that are not consistent with the properties of the solution that are deduced from the above analysis.