Measurement of Aerodynamic Characteristics Using Cinder Models through Free Fall Experiment

: Rocks ejected from a volcanic eruption often cause loss of lives and structures. Aerodynamic characteristics are needed for evaluating motions of volcanic rocks for the reduction of damage. Falling motions of volcanic rock were measured by using models imitated the conﬁguration of cinders collected at the site of the experiment, Sakurajima volcano. Two types, one with sharp edges and one without sharp edges, were selected as representative of cinder and a sphere was selected as reference model. The falling motions of the models dropped down from a drone were recorded by video camera and a stand-alone measuring system that included a pressure sensor, acceleration and angular velocity sensors in the models. The motion, posture, velocity and acceleration of the model were obtained in order to measure the three-dimensional falling trajectory. The drag and the deviation angle between relative wind direction and wind force direction were examined. The variation of the drag coefﬁcient and the deviation angle with Reynolds number was clariﬁed.


Introduction
The eruption of Mt. Ontake in September 2014 caused 63 deaths and missing persons [1]. Most victims died by flying rocks: lapilli and volcanic blocks. The flying rocks not only cause loss of lives but also loss of structures. Shelters have been constructed near volcanic crater to avoid impact damage caused by fall of the rocks. We have to know the magnitude of impact force to design the shelter but there is little information on aerodynamic characteristics of flying rocks to estimate the flying range and the impact force of the rocks. Rocks ejected from a volcanic eruption have various shapes. However, several studies [2][3][4] assume simple shape for the flying rocks and constant drag coefficient for calculating the flying trajectory.
Alatorre-Ibargüengoitia and Delgado-Granados [5] selected real cinders and measured drag force in a subsonic wind tunnel and obtained the variation of drag coefficient with Reynolds number. On the other hand, many wind tunnel studies on simple shaped objects such as spheres, cubes, cylinders and so on have been carried out and aerodynamic characteristics have been clarified. For example, Tachikawa [6] found two-dimensional aerodynamic characteristics for flat plates. Okazaki and Maruyama [7] used a six-component force balance to obtain the three-dimensional static aerodynamic characteristics of square flat plates and rood tiles. Richards et al. [8] used a six-component force balance to obtain the three-dimensional aerodynamic characteristics of rectangular plates and long rods.
Many of these measurements were done by fixing the flying body in the wind tunnel under a uniform flow. However, the objects usually change posture in the fluctuating wind flow with the change of wind direction in case of real volcanic eruption. Matsui et al. [9] measure the unsteady aerodynamic characteristics of falling bodies by 50 m free fall experiment under calm wind condition inside a dome. They used a stand-alone measuring system that assembled acceleration and angular velocity sensors in flying models such as a cube, a rod and a plate.
We arranged a free fall experiment using the advanced stand-alone measuring system that is robust against impact force. This measuring system similar to that used in Matsui's experiment [9] and we added a pressure sensor. We measured the free fall trajectories and the motion of falling bodies in the natural wind. The unsteady aerodynamic characteristics were analyzed and the variation of the drag coefficient and the deviation angle with Reynolds number was clarified. These aerodynamic parameters are essential for the accurate prediction of trajectory and impact velocity of flying body, which is expected to contribute to the design of volcanic shelters, evacuation action plans and so on.

Free Fall Experiment
A free fall experiment was carried out in Kurokamijigokugawara, Sakurajima, Kagoshima Prefecture, Japan ( Figure 1). There was nothing blocking the view because there were no plants or buildings around. The trajectory of the cinder model falling could be easily obtained and safety was guaranteed.
obtain the three-dimensional aerodynamic characteristics of rectangular plates and long rods. Many of these measurements were done by fixing the flying body in the wind tunnel under a uniform flow. However, the objects usually change posture in the fluctuating wind flow with the change of wind direction in case of real volcanic eruption. Matsui et al. [9] measure the unsteady aerodynamic characteristics of falling bodies by 50 m free fall experiment under calm wind condition inside a dome. They used a stand-alone measuring system that assembled acceleration and angular velocity sensors in flying models such as a cube, a rod and a plate.
We arranged a free fall experiment using the advanced stand-alone measuring system that is robust against impact force. This measuring system similar to that used in Matsui's experiment [9] and we added a pressure sensor. We measured the free fall trajectories and the motion of falling bodies in the natural wind. The unsteady aerodynamic characteristics were analyzed and the variation of the drag coefficient and the deviation angle with Reynolds number was clarified. These aerodynamic parameters are essential for the accurate prediction of trajectory and impact velocity of flying body, which is expected to contribute to the design of volcanic shelters, evacuation action plans and so on.

Free Fall Experiment
A free fall experiment was carried out in Kurokamijigokugawara, Sakurajima, Kagoshima Prefecture, Japan ( Figure 1). There was nothing blocking the view because there were no plants or buildings around. The trajectory of the cinder model falling could be easily obtained and safety was guaranteed.
The falling motion of volcanic rocks were measured by using models which imitated the configuration of cinders collected at the site of the experiment. The cinder models were lifted up to an altitude of about 150 m by a drone and dropped as shown in Figure 2. Falling trajectories of the cinder models were recorded by video cameras. Falling motions: acceleration, angular velocity and pressure, were measured by the stand-alone measuring system.  The falling motion of volcanic rocks were measured by using models which imitated the configuration of cinders collected at the site of the experiment. The cinder models were lifted up to an altitude of about 150 m by a drone and dropped as shown in Figure 2. Falling trajectories of the cinder models were recorded by video cameras. Falling motions: acceleration, angular velocity and pressure, were measured by the stand-alone measuring system.

Cinder Models
Several shapes were observed in the corrected cinders in the experimental field. After visual inspection, two representative types of shapes were selected for the cinder model. The models imitated the shape of the randomly selected cinders. One has rough texture without sharp edges and its overall shape is relatively smooth, as model A shown in Figure 3. The other, model B, has sharp edges on the surface. We used 3D scanning data of actual cinders to make cinder models. The instantaneous (unaveraged) projected area to the wind direction will be changed when the falling cinder model rotates. This means that our experimental results with six models are equivalent to those using cinder models with various shapes. Cinder models were made of polyurethane foam and were 20 cm to 30 cm in diameter. A sphere, made of polystyrene foam with 30 cm in diameter, was used for comparative verification because lots of information on aerodynamic characteristics were experimented in the previous study. Table 1 shows the specification of the models: model mass, projected area and equivalent diameter, for six test cases. The projected area A was obtained as an averaged value of projected area for overall angle of view. The equivalent diameter L is obtained as 2 .

Cinder Models
Several shapes were observed in the corrected cinders in the experimental field. After visual inspection, two representative types of shapes were selected for the cinder model. The models imitated the shape of the randomly selected cinders. One has rough texture without sharp edges and its overall shape is relatively smooth, as model A shown in Figure 3. The other, model B, has sharp edges on the surface. We used 3D scanning data of actual cinders to make cinder models. The instantaneous (unaveraged) projected area to the wind direction will be changed when the falling cinder model rotates. This means that our experimental results with six models are equivalent to those using cinder models with various shapes. Cinder models were made of polyurethane foam and were 20 cm to 30 cm in diameter. A sphere, made of polystyrene foam with 30 cm in diameter, was used for comparative verification because lots of information on aerodynamic characteristics were experimented in the previous study. Table 1 shows the specification of the models: model mass, projected area and equivalent diameter, for six test cases. The projected area A was obtained as an averaged value of projected area for overall angle of view. The equivalent diameter L is obtained as 2 A π .

Wind Speed
The vertical distribution of wind speed was measured by Doppler lidar (Mitsubishi in Tokyo, Japan, DIABREZZA) during the free fall experiment of all cases as shown in Figure 4. The 20 s mean wind speeds show a gentle breeze condition. These wind speeds are used to calculate the drag coefficient and the deviation angle in Section 3.2. The maximum absolute wind speed of 4 m/s can be found in case B1.

Wind Speed
The vertical distribution of wind speed was measured by Doppler lidar (Mitsubishi in Tokyo, Japan, DIABREZZA) during the free fall experiment of all cases as shown in      Figure 5 shows the stand-alone measuring system in the cinder model. The system was composed of 3-axis acceleration and angular velocity sensors (Inven Sense in California, USA, MPU-9250), a pressure sensor (Bosch Sensortec in Reutlingen, Germany, GY-BMP280-3.3), a microcomputer board (Raspberry Pi in Cambridge, England, Model A+), an SD memory card, a switch giving an ON/OFF signal of measurement and a battery. Specification of acceleration and angular velocity sensors were shown in Table 2. These measuring devices were put in a box. Rods were attached to the xy-axis of the box so that the measuring devices were fixed in the cinder model. Cotton was stuffed around the box. A hook was attached to the cinder model so that the drone could be easily lifted.     The cinder model was lifted up by a drone (Luce Search in Hiroshima, Japan, SPIDER CS-6) as shown in Figure 6. The drone was equipped with a device for falling model and a smartphone with GPS for recording location. A flight record of drones during the whole day of the experiment shows the altitude of dropping cinder models at about 150 m from the ground as shown in Figure 7.    The pressure sensor record shown in Figure 8 indicates an altitude of about 150 m for the falling test of case A2 as shown in Figure 9. The records of acceleration and angular velocity sensors of case A2 shown in Figures 10 and 11.

Video Motion Analysis
The falling motion of the cinder model was recorded by four video cameras with pixels of 3840 × 2160 and frame rate of 29.97 fps (Panasonic in Osaka, Japan, HC-VX985M). The video cameras were arranged at four points A, B, C and N as shown in Figure 12a by motion analysis software (DITECT in Tokyo, Japan, Dipp-Motion V). The Direct Linear Transformation method [10] was used as the image analysis algorithm. Previous to the falling experiment, the drone flew to record the reference points shown as red symbols • in Figure 12a

Motion Analysis
Acceleration and angular velocity in the body coordinate system on the falling model were obtained by the sensors. An example of case A2 was shown in Figures 10 and 11. The bias of angular velocity was corrected. The acceleration and angular velocity were filtered with a 5Hz Butterworth filter.
We converted acceleration and angular velocity into the value in a right-handed absolute coordinate system aligned as X-axis directing to the east, Y-axis to the north and Zaxis to the vertical. The quaternion and the rotation matrix were used for the conversion

Motion Analysis
Acceleration and angular velocity in the body coordinate system on the falling model were obtained by the sensors. An example of case A2 was shown in Figures 10 and 11. The bias of angular velocity was corrected. The acceleration and angular velocity were filtered with a 5Hz Butterworth filter.
We converted acceleration and angular velocity into the value in a right-handed absolute coordinate system aligned as X-axis directing to the east, Y-axis to the north and Z-axis to the vertical. The quaternion and the rotation matrix were used for the conversion [11]. Finally, the falling velocity, the displacement and the angle of the falling body were obtained by integration with trapezoidal rule.
The initial posture cannot be determined because the start time of the model falling is uncertain. The initial posture has been optimized so that the trajectory integrated from the record of sensor best matches the trajectory obtained by video analysis.
The trajectories obtained by sensor and video analysis in the absolute system are shown in Figure 13. Error of three-axis displacement by sensor and video analysis of all cases are shown in Figure 14. The maximum absolute value of error is about 2.5 m. The error is within 2% with respect to the fall height of 150 m. There are good agreements which ensure the accuracy of the motion analysis. The cinder model went out of the camera view in case C1, so the trajectory by video analysis was broken at the halfway height. the record of sensor best matches the trajectory obtained by video analysis.
The trajectories obtained by sensor and video analysis in the absolute system are shown in Figure 13. Error of three-axis displacement by sensor and video analysis of all cases are shown in Figure 14. The maximum absolute value of error is about 2.5 m. The error is within 2% with respect to the fall height of 150 m. There are good agreements which ensure the accuracy of the motion analysis. The cinder model went out of the camera view in case C1, so the trajectory by video analysis was broken at the halfway height.

Falling Motion
The displacement, acceleration, velocity and angular velocity in the absolute coordi-

Falling Motion
The displacement, acceleration, velocity and angular velocity in the absolute coordinate system were obtained. An example of case A2 is shown in Figure 15. The displacement calculated by sensor shows good agreement with video analysis as shown in Figure 15a. After the falling, the vertical acceleration a Z : z component increased and was around 0 m/s 2 after 3 s as shown in Figure 15b. This means that the wind resistance on the cinder model increased and balanced with the gravity, and the velocity reached terminal velocity: −17 m/s, at 2.5 s as shown in Figure 15c. The magnitude of angular velocity of the vertical component ω Z gradually increased after 4 s as shown in Figure 15d. The periodic variation of X and Y components of angular velocity ω X and ω Y showed the precession along the Z axis.

Falling Motion
The displacement, acceleration, velocity and angular velocity in the absolute coordinate system were obtained. An example of case A2 is shown in Figure 15. The displacement calculated by sensor shows good agreement with video analysis as shown in Figure  15a. After the falling, the vertical acceleration : z component increased and was around 0 m/s 2 after 3 s as shown in Figure 15b. This means that the wind resistance on the cinder model increased and balanced with the gravity, and the velocity reached terminal velocity: −17 m/s, at 2.5 s as shown in Figure 15c. The magnitude of angular velocity of the vertical component gradually increased after 4 s as shown in Figure 15d. The periodic variation of X and Y components of angular velocity and showed the precession along the Z axis.

Drag Coefficient and Deviation Angle
We derived the wind resistance force from the acceleration = ( , , ) as where m is the mass and g is the gravitational acceleration as 9.794 m/s 2 . The drag coefficient was derived as

Drag Coefficient and Deviation Angle
We derived the wind resistance force → F from the acceleration → a = (a X , a Y , a Z ) as where m is the mass and g is the gravitational acceleration as 9.794 m/s 2 . The drag coefficient C d was derived as where A is the projected area of the falling body and the dynamic pressure q is defined as where ρ is the air density as 1.2 kg/m 3 and → u = (u X , u Y , u Z ) is the relative velocity defined as the difference of the wind speed → U and the velocity of falling body As shown in Figure 16. The Doppler lidar cannot measure wind speed → U below 50 m. So the wind speed at 50 m was used for the calculation below 50 m. The deviation angle | | between wind direction and wind resistant force direction was derived as Figure 17 shows the variation of drag coefficient with Reynolds number = / for different test cases, where is the representative length shown in Table 1  and is the air viscosity as 1.8 10 −5 Pa•s. varies with the configuration of the cinder models in the range of low . Also varies with the falling motion because the instantaneous projected area to the wind direction is changed. The dispersion of decreases with and may converge to a constant value of around 0.4 that is a little bit larger than that of sphere.
The deviation angle | | varies with configuration of the cinder models and the instantaneous projected area to the wind direction as shown in Figure 18. The magnitude and dispersion of | | decrease with the Reynolds number . In the range of high , | | of the cinder model varies within 22 degrees, independent of . Otherwise | | of the sphere varies within 31 degrees. The deviation angle |θ| between wind direction and wind resistant force direction was derived as (5) Figure 17 shows the variation of drag coefficient C d with Reynolds number R e = ρL| → U|/µ for different test cases, where L is the representative length shown in Table 1 and µ is the air viscosity as 1.8 × 10 −5 Pa·s. C d varies with the configuration of the cinder models in the range of low R e . Also C d varies with the falling motion because the instantaneous projected area to the wind direction is changed. The dispersion of C d decreases with R e and may converge to a constant value of around 0.4 that is a little bit larger than that of sphere. larger than that of sphere.
The deviation angle | | varies with configuration of the cinder models and the instantaneous projected area to the wind direction as shown in Figure 18. The magnitude and dispersion of | | decrease with the Reynolds number . In the range of high , | | of the cinder model varies within 22 degrees, independent of . Otherwise | | of the sphere varies within 31 degrees. The deviation angle |θ| varies with configuration of the cinder models and the instantaneous projected area to the wind direction as shown in Figure 18. The magnitude and dispersion of |θ| decrease with the Reynolds number R e . In the range of high R e , |θ| of the cinder model varies within 22 degrees, independent of R e . Otherwise |θ| of the sphere varies within 31 degrees. Figure 18. Variation of deviation angle |θ| with Reynolds number. Note that unit of θ in Equation (5) is the radian. Unit of y-axis θ in Figure 18 is the degree.

Discussion
To do ballistic simulation, we need to know magnitude and direction of wind force of the cinder. There are various shapes of actual cinders. The magnitude and direction of wind force varies with configuration of the shapes and the instantaneous projected area of the wind. The projected area will be changed when the falling cinder model rotates. This means that our experimental results with six models include results of used cinder models with various shapes. The magnitude and direction of wind force can be calculated by the drag coefficient and deviation angle obtained in our study.
The drag coefficient C d and deviation angle |θ| varies with the cinder configuration and the instantaneous projected area with rotation of the cinder. The magnitude and dispersion of C d and |θ| decrease with the Reynolds number R e as obtained in this study.
The results of C d and |θ| in this study include the effect of cinder model rotation. Therefore, the ballistic simulation with them includes the effect of rotation. We can calculate the trajectory without calculating the dynamics of rotation.
The variations of C d and |θ| do not show obvious dispersion at the same R e , so we may use the constant frequency distribution for C d and |θ|. We will try to have further study on this point.

Conclusions
We conducted the free fall experiment from a drone using cinder models. Falling motion of the cinder models were recorded by video camera and a stand-alone measuring system that was assembled with a pressure sensor, acceleration and angular velocity sensors in the cinder models. Three-dimensional falling trajectory and motion, including displacement, acceleration, velocity and angular velocity, were calculated. The variation of the drag coefficient and the deviation angle with Reynolds number were obtained.
The drag coefficient varies with the configurations of falling body and the falling motion when the Reynolds number is small. The drag coefficient converges into 0.4 that is a little bit larger than that of sphere's value when the Reynolds number is large. The deviation angle of cinder model is within 22 degrees independent of the configurations; that of the sphere varies within 31 degrees when Reynolds number is large.
To do ballistic simulation, it is necessary to determine the frequency distribution of the drag coefficient and deviation angle at the same Reynolds number. We may use the constant frequency distribution for the drag coefficient and deviation angle at the same Reynolds number because the variations of drag coefficient and deviation angle show no obvious dispersion.
The ballistic simulation with the drag coefficient and deviation angle obtained include effect of rotation. It will be convenient to calculate ballistic trajectory without calculating the dynamics of rotation.