Microdrone for Atmospheric Sensing -- CFD Animation

Today, I am going to share with you how to simulate the freefall motion in CFD solver ANSYS FLUENT. The example given is a rotor, a prototype of atmospheric sensor. The domain is divided into the central area, surrounding the geometry, and the peripheral area. We are going to define a sliding mesh interface between the central and peripheral area:

Turbulence model is set to be k epsilon with an enhanced wall treatment. The gravity is set to be zero: it will be replaced by an external force added to the microdrone.

Now, the only thing left is to craft a UDF that automatically updates the angular velocity of the drone and its falling velocity in the air. This UDF is somewhat tricky: 1) It requires the use of a in built function ‘Compute_Force_And_Moment’ which has no documentation; 2) Some data transmission has to be allowed between two different macros. So here, the source code is provided with a step-by-step illustration:

1) Include ‘udf.h’ and define constants

As usual, udf.h is included at the beginning. We also define three constants here to increase the readability.

#include "udf.h"

#define WEIGHT 0.3434335

#define IZ 1.744e-4

#define MASS 35e-3

2) Initialize three static variables

Most UDF scripts are executed after each time step. So if you want to update a variable based on its previous value, it is highly recommended to claim them as static. Static variables would only be initialized once, keeping its value unless it is updated explicitly.

static real wind_speed = 0.01; /*Falling speed, always positive*/

static real w = 0.0; /*Angular velocity*/

static real lift = 0.0; /*Lift produced by air*/

3) Update the rotation of drone

Here we call DEFINE_ZONE_MOTION to modify the angular velocity of central area. To do so, we use Compute_Force_And_Moment. The force and the moment will be stored in f_drone and m_drone separately. The angular velocity is updated by Euler explicit scheme.

/*Updating the rotation of drone*/

DEFINE_ZONE_MOTION(drone_rotation, omega, axis, origin, velocity, time, dtime)

{

Domain *domain = Get_Domain(1); /*Get the pointer to fluid domain*/

Thread *tf = Lookup_Thread(domain, 14); /*Thread pointing to the drone, 14 is the ID of drone wall in this case*/

real x_cg[ND_ND];

real f_drone[ND_ND];

real m_drone[ND_ND];

N3V_D(x_cg, =, 0.0, 0.0, 0.0);

Compute_Force_And_Moment(domain, tf, x_cg, f_drone, m_drone, TRUE);

lift = f_drone[2];

/*Updating the force applied on drone, gravity included*/

real force;

lift = f_drone[2];

force = lift - WEIGHT;

/*Updating the motion of frame*/

w = w + m_drone[2] / IZ * dtime;

N3V_S(omega, =, w);

N3V_D(velocity, =, 0.0, 0.0, 0.0);

N3V_S(origin, =, 0.0); /* default values, line could be omitted */

N3V_D(axis, =, 0.0, 0.0, 1.0); /* default values, line could be omitted */

/*Updating the wind speed*/

wind_speed = wind_speed - force / MASS * dtime;

/*Writing data to file*/

FILE*fp;

fp = fopen("data.txt", "a");

fprintf(fp, "%f,%f,%f,%f,%f\n", time, w, wind_speed, lift,m_drone[2]);

fclose(fp);

return;

}

4) Update the wind speed

We calculated the wind speed in function DEFINE_ZONE_MOTION because it is used here to update the speed at the velocity inlet. The speed is also determined by Euler explicit scheme.

/*Updating the wind speed*/

DEFINE_PROFILE(velocity_inlet, t, i)

{

face_t f;

begin_f_loop(f, t)

{

F_PROFILE(f, t, i) = wind_speed;

}

end_f_loop(f,t)

}

One thing that needs to mention is that instead of changing the velocity of microdrone, we are actually adjusting the velocity of the air. This helps to save a lot of computational resources: if we let the drone fall in a static domain, we have to create a cylindrical domain that is of the same height as initial releasing point (10 km in our case). By putting the reference frame on the drone, we are allowed to shorten the domain and make our case calculable on a desktop.

An animation of the case has been made and shown in the Design Expo at University of Michigan: