Help with force/drag calculation for maximum airspeed calculation

I'm trying to do a rough calculation of the maximum airspeed attainable by a particular aircraft. I've completed some force simulations in other sims before, but not one where I have a set amount of available thrust (40 kN) and plan to increase the lbm_u flow speed until the total drag force in the y-axis equals this amount.

The problem is I can't really seem to get the conversion to/from real life to sim units correct in the code (or in my head) as far as I can tell + Primarily I'm not sure of how to determine the real life airspeed based on the _lbmu flow speed.

(The results are also very unstable and can't get the simulation to converge on an accurate value, but one problem at a time 🤔)

I've read right through these other issue threads (https://github.com/ProjectPhysX/FluidX3D/issues/36, https://github.com/ProjectPhysX/FluidX3D/issues/32) relating to force readouts but I think I'm still missing something here. My code and force results are below:

void main_setup() {

    // setup the volume and objects
    const uint L = 512u;                // size of test volume
    const uint Nx = to_uint(L * 1.1f);      // adjust size of volume x axis
    const uint Ny = to_uint(L * 1.1f);      // adjust size of volume y axis
    const uint Nz = to_uint(L * 0.4f);      // adjust size of volume z axis
    const float size = 0.5f * (float)L;     // scaling of object compared with size of the volume

    // setup aircraft parameters
    const float knots = 300.0f;                     // speed through air (only affects Reynolds number I think?)
    const float AoA = 0.0f;                         // angle of attack (°), rotates object in field on x-axis
    const float3 center = float3(0.5f * Nx, 0.45f * Ny, 0.55f * Nz);    // offset the aircraft position
    const float3x3 rotation = float3x3(float3(1, 0, 0), radians(-AoA));     // set the aircraft rotation, check the specified axis

    // setting SI units
    const float si_x = 18.0f;           // characteristic length (m)
    const float si_u = knots * 0.5144f;     // convert airspeed to SI units (m/s)
    const float si_rho = 1.2f;          // air density 
    const float si_nu = 1.48E-5f;           // kinematic shear viscosity of air (~sea level)

    // setting LBM units
    const float lbm_u = 0.12f;
    const float lbm_rho = 1.0f;         // density in LBM units always 1
    const float lbm_x = L;

    // set unit conversion factors between SI and LBM units
    units.set_m_kg_s(lbm_x, lbm_u, lbm_rho, si_x, si_u, si_rho);
    const float lbm_nu = units.nu(si_nu); // kinematic shear viscosity in LBM units

    // create LBM object, setup mesh and set flags
    LBM lbm(Nx, Ny, Nz, lbm_nu);

    // setup mesh
    Mesh* mesh = read_stl(get_exe_path() + "../stl/aircraft.stl", lbm.size(), center, rotation, size);
    lbm.voxelize_mesh_on_device(mesh);
    lbm.flags.read_from_device();

    const uint N = lbm.get_N();
    for (uint n = 0ull, x = 0u, y = 0u, z = 0u; n < N; n++, lbm.coordinates(n, x, y, z)) {
        if (lbm.flags[n] != TYPE_S) lbm.u.y[n] = lbm_u;
        if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) lbm.flags[n] = TYPE_E;
    }

    // set overlay options
    key_1 = true;
    key_2 = false;
    key_3 = true;
    key_4 = true;

    // setup text file output
    std::ofstream fout(get_exe_path() + "results.txt");
    if (!fout) {
        std::cerr << "Could not open file." << std::endl;
    }

    // run the simulation
    lbm.run(0u);
    while (lbm.get_t() < 30000u) {

        lbm.run(100u);

        // calculate the force on the object
        lbm.calculate_force_on_boundaries();
        lbm.F.read_from_device();
        const float3 force = lbm.calculate_force_on_object(TYPE_S);

        // write the calculated force to file
        fout << to_string(units.si_F(force.y)) << std::endl;
        print_info(to_string(units.si_F(force.y)));
    }

    // close the output file
    fout.close();
    wait(); // wait for a keypress to close the Program
}

My defines header here:

```cpp #pragma once //#define D2Q9 // choose D2Q9 velocity set for 2D; allocates 53 (FP32) or 35 (FP16) Bytes/node //#define D3Q15 // choose D3Q15 velocity set for 3D; allocates 77 (FP32) or 47 (FP16) Bytes/node //#define D3Q19 // choose D3Q19 velocity set for 3D; allocates 93 (FP32) or 55 (FP16) Bytes/node; (default) #define D3Q27 // choose D3Q27 velocity set for 3D; allocates 125 (FP32) or 71 (FP16) Bytes/node #define SRT // choose single-relaxation-time LBM collision operator; (default) //#define TRT // choose two-relaxation-time LBM collision operator //#define FP16S // compress LBM DDFs to range-shifted IEEE-754 FP16; number conversion is done in hardware; all arithmetic is still done in FP32 //#define FP16C // compress LBM DDFs to more accurate custom FP16C format; number conversion is emulated in software; all arithmetic is still done in FP32 //#define BENCHMARK // disable all extensions and setups and run benchmark setup instead //#define VOLUME_FORCE // enables global force per volume in one direction, specified in the LBM class constructor; the force can be changed on-the-fly between time steps at no performance cost #define FORCE_FIELD // enables computing the forces on solid boundaries with lbm.calculate_force_on_boundaries(); and enables setting the force for each lattice point independently (enable VOLUME_FORCE too); allocates an extra 12 Bytes/node //#define MOVING_BOUNDARIES // enables moving solids: set solid nodes to TYPE_S and set their velocity u unequal to zero #define EQUILIBRIUM_BOUNDARIES // enables fixing the velocity/density by marking nodes with TYPE_E; can be used for inflow/outflow; does not reflect shock waves //#define SURFACE // enables free surface LBM: mark fluid nodes with TYPE_F; at initialization the TYPE_I interface and TYPE_G gas domains will automatically be completed; allocates an extra 12 Bytes/node //#define TEMPERATURE // enables temperature extension; set fixed-temperature nodes with TYPE_T (similar to EQUILIBRIUM_BOUNDARIES); allocates an extra 32 (FP32) or 18 (FP16) Bytes/node #define SUBGRID // enables Smagorinsky-Lilly subgrid turbulence model to keep simulations with very large Reynolds number stable #define INTERACTIVE_GRAPHICS // enable interactive graphics; start/pause the simulation by pressing P; either Windows or Linux X11 desktop must be available; on Linux: change to "compile on Linux with X11" command in make.sh //#define INTERACTIVE_GRAPHICS_ASCII // enable interactive graphics in ASCII mode the console; start/pause the simulation by pressing P //#define GRAPHICS // run FluidX3D in the console, but still enable graphics functionality for writing rendered frames to the hard drive #define GRAPHICS_FRAME_WIDTH 3840 // set frame width if only GRAPHICS is enabled #define GRAPHICS_FRAME_HEIGHT 2160 // set frame height if only GRAPHICS is enabled #define GRAPHICS_BACKGROUND_COLOR 0x000000 // set background color; black background (default) = 0x000000, white background = 0xFFFFFF #define GRAPHICS_U_MAX 0.18f // maximum velocity for velocity coloring in units of LBM lattice speed of sound (c=1/sqrt(3)) (default: 0.15f) #define GRAPHICS_Q_CRITERION 0.0001f // Q-criterion value for Q-criterion isosurface visualization (default: 0.0001f) #define GRAPHICS_BOUNDARY_FORCE_SCALE 100.0f // scaling factor for visualization of forces on solid boundaries if VOLUME_FORCE is enabled and lbm.calculate_force_on_boundaries(); is called (default: 100.0f) #define GRAPHICS_STREAMLINE_SPARSE 47 // set how many streamlines there are every x lattice points #define GRAPHICS_STREAMLINE_LENGTH 256 // set maximum length of streamlines // ############################################################################################################# #define TYPE_S 0b00000001 // (stationary or moving) solid boundary #define TYPE_E 0b00000010 // equilibrium boundary (inflow/outflow) #define TYPE_T 0b00000100 // temperature boundary #define TYPE_F 0b00001000 // fluid #define TYPE_I 0b00010000 // interface #define TYPE_G 0b00100000 // gas #define TYPE_X 0b01000000 // reserved type X #define TYPE_Y 0b10000000 // reserved type Y #if defined(FP16S) || defined(FP16C) #define fpxx ushort #else // FP32 #define fpxx float #endif // FP32 #ifdef BENCHMARK #undef UPDATE_FIELDS #undef VOLUME_FORCE #undef FORCE_FIELD #undef MOVING_BOUNDARIES #undef EQUILIBRIUM_BOUNDARIES #undef SURFACE #undef TEMPERATURE #undef SUBGRID #undef INTERACTIVE_GRAPHICS #undef INTERACTIVE_GRAPHICS_ASCII #undef GRAPHICS #endif // BENCHMARK #ifdef SURFACE // (rho, u) need to be updated exactly every LBM step #define UPDATE_FIELDS // update (rho, u, T) in every LBM step #endif // SURFACE #ifdef TEMPERATURE #define VOLUME_FORCE #endif // TEMPERATURE #if defined(INTERACTIVE_GRAPHICS) || defined(INTERACTIVE_GRAPHICS_ASCII) #define GRAPHICS #endif // INTERACTIVE_GRAPHICS || INTERACTIVE_GRAPHICS_ASCII ```

Here's a snippet of the force results output (does anyone know how to get the output to not be scientific notation? Google Sheets doesn't handle these values and I have to resort to MS Excel):

Force Calculation Results

### Taken each 100 steps ``` 5.45981456E4 7.24956846E4 7.05784989E4 4.87261344E4 5.33241606E4 4.76977920E4 4.39427472E4 4.65264656E4 4.49056340E4 4.15205240E4 4.89547016E4 4.37330200E4 3.41821480E4 5.59037636E4 3.54986192E4 3.03640080E4 6.60130216E4 1.98685264E4 4.37987564E4 5.84950160E4 1.40907776E4 5.30362702E4 4.93470624E4 1.66901088E4 6.21138572E4 3.50296092E4 2.53584744E4 6.05184031E4 2.65506504E4 3.66595888E4 5.78224232E4 2.07870364E4 4.56666420E4 4.79558512E4 2.18530130E4 5.24934960E4 4.20690488E4 2.38719248E4 5.73561760E4 3.04396701E4 3.29020452E4 5.45614432E4 2.86208776E4 3.78154800E4 5.19157028E4 2.27532816E4 4.62905692E4 4.43532228E4 2.78953000E4 4.90866040E4 3.82027696E4 2.95206312E4 5.04203034E4 3.29360580E4 3.73151400E4 4.77963496E4 3.10918498E4 3.99938512E4 4.41328812E4 2.99023864E4 4.72581864E4 3.70907136E4 3.52803944E4 4.51192808E4 3.29675292E4 3.91297720E4 4.56213856E4 3.03644109E4 4.48174240E4 3.93006824E4 3.07123041E4 4.83290144E4 3.47505568E4 3.55810620E4 4.79376888E4 2.96577928E4 3.90675064E4 4.69986776E4 2.75215248E4 4.63786412E4 4.14923382E4 2.75947616E4 4.86272048E4 3.66286016E4 3.10061097E4 5.13841009E4 3.11370969E4 3.47928692E4 5.00028658E4 2.79550312E4 3.98808744E4 4.84784648E4 2.45359184E4 4.63440372E4 4.16207886E4 2.70788336E4 4.90002344E4 3.83687952E4 2.75404432E4 5.26642036E4 3.07034731E4 3.39467932E4 5.09975672E4 2.84062696E4 3.78370336E4 5.01098347E4 2.38102340E4 4.47991560E4 4.50371884E4 2.43556548E4 4.90558720E4 3.96822504E4 2.49006056E4 5.38041876E4 3.20686198E4 3.12109590E4 5.30605126E4 2.81475904E4 3.52483344E4 5.24868106E4 2.32955694E4 4.32908916E4 4.73114536E4 2.31392718E4 4.78467272E4 4.26806164E4 2.36173772E4 5.27146720E4 3.64221096E4 2.72658848E4 5.37688492E4 3.21892882E4 3.15374017E4 5.30127192E4 2.87015728E4 3.51730704E4 5.13758898E4 2.62463404E4 4.04599190E4 4.75241136E4 2.80091504E4 4.12005997E4 4.52993488E4 2.95608160E4 4.34895752E4 4.26450872E4 3.17720150E4 4.10161400E4 4.15983200E4 3.45705508E4 4.11882591E4 4.11889410E4 3.62099792E4 3.64361740E4 4.39449404E4 3.55373500E4 3.74828080E4 4.55044604E4 3.47750972E4 3.48780344E4 4.73938992E4 3.23299218E4 3.82154488E4 4.80530736E4 3.00424647E4 3.91072776E4 4.62211560E4 2.82382848E4 4.31937932E4 4.56638528E4 2.68889168E4 4.58463956E4 4.08407545E4 2.79705000E4 4.83229064E4 3.97819688E4 2.86021256E4 4.94393400E4 3.49910664E4 3.12321615E4 5.04131508E4 3.44256544E4 3.30876518E4 5.02775574E4 2.88833376E4 3.78664400E4 4.77260304E4 3.04172182E4 4.01574278E4 4.58602380E4 2.68378544E4 4.47494124E4 4.18165016E4 3.05465555E4 4.53165436E4 3.94830728E4 2.89034680E4 4.83971688E4 3.50941636E4 3.48705148E4 4.67443224E4 3.39380216E4 3.38679432E4 4.85436392E4 3.05737758E4 4.04101801E4 4.52647828E4 3.02574444E4 4.03174400E4 4.44756936E4 2.88578224E4 4.54932880E4 4.09248829E4 2.99064352E4 4.50767516E4 3.92672752E4 3.06318474E4 4.82059480E4 3.57805896E4 3.25918674E4 4.66934728E4 3.37587928E4 3.48344444E4 4.87301160E4 3.17187286E4 3.70109248E4 4.61975004E4 2.93210528E4 4.05846787E4 4.57513000E4 2.92939016E4 4.27224446E4 4.21224594E4 2.82099104E4 4.62798644E4 3.98973344E4 3.09912205E4 4.61989500E4 3.70167640E4 3.05215764E4 4.95039272E4 3.33795764E4 3.57161308E4 4.69126512E4 3.14432621E4 3.61775900E4 4.83867120E4 2.86830904E4 4.32326270E4 4.34546520E4 2.88929984E4 4.28917074E4 4.26834106E4 2.88632296E4 4.81515312E4 3.75875688E4 3.07526446E4 4.73358968E4 3.54067968E4 3.24639488E4 5.02047062E4 3.08033443E4 3.67542480E4 4.71359352E4 3.03295064E4 ```

So the forces seem in the ballpark of 40 kN, but what would I say the real-life airspeed of the simulation is here?

Screenshots:

Plenty of volume around the aircraft: ![Bare](https://user-images.githubusercontent.com/94650926/215614326-f0ba8596-0992-428c-b6e3-b49dceb52a07.png) Early in the sim (force lines start flailing around further into the simulation): ![Underway](https://user-images.githubusercontent.com/94650926/215614334-e040e10c-021b-4185-a8b7-dfbd4f79f7cd.png)

The results are all over the place and I'm completely stuck now. Can anyone see where I'm going wrong and/or how I can get the results I'm looking for? 😟

Greetings!

To calculate the drag force of 40kN, you can start by using simple geometry such as a sphere or a vector cylinder in 2D or 3D. You will need the drag coefficient, air density, and the geometry to calculate the air forces. Once you have calculated the air forces, you can convert them to LBM units.

If you need an example, let me know and I can show you some. Also, if you have an STL file, you can use it to help with the calculation of dry forces.

Regarding your air wing, you mentioned that it is four instead of two which is an interesting point. It might take some time to calculate, but you can give it a try.

Let me know if there's anything else I can help with!

@randomwangran An example would be great! Regarding the number of number of wings, etc. I'm primarily (at this stage) interested only with calculating the frictional (parasite) drag as the induced drag decreases significantly at higher airspeeds, so the performance of the wing isn't of much concern, only the amount of frictional drag it creates... Unless, are you saying that a C_d which includes induced drag needs to be calculated first? I'm sure if I see an example that it'll clarify the workflow for me. Thanks so much!

RE: https://github.com/ProjectPhysX/FluidX3D/issues/55#issuecomment-1412738676 @randomwangran Would it be possible to see your example?

I think there might be multiple things going on here. 1.) Your mach number is about 0.44 correct? Making the assumption that air is incompressible (which we are with FluidX3D) with this mach number might be causing issues. I think you can go up to 0.3, maybe even 0.4 in some situations, but a lower mach number would be better. I would try a lower airspeed to lower your mach number. 2.) You're using a domain with about 64M voxels right? Looking at your screenshot the voxels flagged as solids make up a small portion of this domain and the mesh resolution looks pretty low. So, a low resolution cartesian mesh will present many flat surfaces orthogonal to the airflow. This could cause very unsteady flow due to vortex shedding off of these flat surfaces, which makes your lift and drag coeffs time dependant and jump around. Its kind of like having a very high angle of attack on an airfoil. For example look at Figure 12 of this paper here. This shows how variable the drag coefficient (and thus the drag force) can be at high angles of attack due to vortex shedding. So, you could try increasing the relative size of the airplane compared to your domain size where its a tradeoff between boundary effects and resolving the object, and/or increase the number of voxels in the domain. Also its worth saying that drag on an entire aircraft moving through the air will be unsteady and jump around a lot if it has a reynolds number that is in the transitional or turbulent flow regime. So, an unsteady drag force is to be expected and isn't too surprising.

I think there might be multiple things going on here. 1.) Your mach number is about 0.44 correct? Making the assumption that air is incompressible (which we are with FluidX3D) with this mach number might be causing issues. I think you can go up to 0.3, maybe even 0.4 in some situations, but a lower mach number would be better. I would try a lower airspeed to lower your mach number. 2.) You're using a domain with about 64M voxels right? Looking at your screenshot the voxels flagged as solids make up a small portion of this domain and the mesh resolution looks pretty low. So, a low resolution cartesian mesh will present many flat surfaces orthogonal to the airflow. This could cause very unsteady flow due to vortex shedding off of these flat surfaces, which makes your lift and drag coeffs time dependant and jump around. Its kind of like having a very high angle of attack on an airfoil. For example look at Figure 12 of this paper here. This shows how variable the drag coefficient (and thus the drag force) can be at high angles of attack due to vortex shedding. So, you could try increasing the relative size of the airplane compared to your domain size where its a tradeoff between boundary effects and resolving the object, and/or increase the number of voxels in the domain. Also its worth saying that drag on an entire aircraft moving through the air will be unsteady and jump around a lot if it has a reynolds number that is in the transitional or turbulent flow regime. So, an unsteady drag force is to be expected and isn't too surprising.

@trparry Thanks for the reply!

1.) Regarding the mach number, I must be mistaken on how it's represented in the code. From what I can make out in defines.hpp, local Mach 1.0 is an lbm speed of $1/\sqrt{3}$ (≈0.577) so the mach number of the lbm speed I set (0.12) would be a Mach of $0.12/0.577≈0.21$ I must have this wrong though? I understand that this method of simulation doesn't cover compressibility, and I'm not expecting that this airframe will be able to approach speeds anywhere near those with compressibility effects. So if the mach number I'm trying to induce here is infact 0.44 then that would certainly explain why I'm getting enormous drag results.

2.) I included screenshots/code in the original post from smaller domains of 64 million voxels, just to quickly get some screenshots and a list of results - but I was getting similar results from domains that would fill up the 12 GB of memory on my 3080 Ti, but I'll refine the relative size of the aircraft in the domain and the domain itself once I resolve the other issues I'm having first.

So I guess firstly I'll need clarification on the speed of sound/mach number in the sim? Thanks in advance

@SirWixy

1.) I am simply taking 300 knots and dividing by the speed of sound in air to get 0.44 mach number. 300 knots = 154 m/s. 154 m/s / 343 m/s = 0.44 Mach.

2.) Yeah I think this is one of the biggest challenges with FluidX3D. Without being able to locally refine the mesh, to sufficiently resolve the boundary layer near the body you have to globally increase the total number of cells by quite a bit to get sufficient accuracy. This offsets some of the other gains in speed that FLuidX3D is known for, if you were to compare it to a software that can use far fewer cells with local mesh refinement to get the same accuracy. For example, look near the end of issue #32, where they used 500k cells in Fluent and 150M in FluidX3D to get the same results. If FluidX3D is able to implement local mesh refinement in the future, then this gets really interesting because not only would FluidX3D be fast, but it would also be efficient in the number of cells required to achieve useful accuracy. For the same level of accuracy, I think even without local mesh refinement FluidX3D is still faster than FVM based code for some applications because a GPU on millions of cells is sometimes still faster than a CPU working on only a few hundred thousand. This demonstrates the pure muscle of a modern GPU. However, local mesh refinement would definitely push this software into sort of "the holy grail" of CFD for incompressible flows.

To be able to sufficiently resolve complex geometry such as an entire airplane, you need more cells. This is why some of the test cases use billions of cells. So, you might consider using a more capable GPU hosted on Amazon AWS, Azure, Google cloud, paperspace, etc. Paperspace has some easy to use pre-configured virtual machines that you can get up and running pretty quickly.

@trparry Oh you're looking at the 300 Knots figure, yes that'd be correct Mach for that airspeed. But before I start adjusting that si_u speed (it's just set to 300 knots as a placeholder for now) to find the drag equilibrium, I'll need to figure out how to set the lattice speed, lbm_u, don't I? Changing the lbm_u value has quite an effect of the results of the sim, so I presume it's not a trivial setting (i.e. not like the lbm_rho value which cancels out anyway) and this is what I was trying to ask in the original post. What should the lbm_u value be set to?

@SirWixy the value of lbm_u can be somewhere in approximately 0.0001-0.2. Accuracy depends on tau, the LBM relaxation time, which is set by converting si_nu to lbm_nu=(tau-0.5)/3 during unit conversion, and that depends on the choice of lbm_u. LBM generally is most accurate at tau=1 (lbm_nu=1/6), and you can work out the unit conversion in reverse to see which lbm_u that would be equivalent to. But achieving tau=1 is not feasible in most cases, as lbm_u would be way too small or way too large then, which leads to excessive floating-point round-off errors (u too small) or instability (u too large).

@ProjectPhysX Okay I might be close now, but what is the unit conversion in reverse of lbm_u? Say if I make lbm_nu = 1/6 (which I may have to adjust, to get an lbm_u in the range of 0.0001-0.2 as you said), then how do I know what lbm_u that is equivalent to? I couldn't find a "float u_from_nu" in units.hpp or any mention in the comments of what's going on 😟

@SirWixy crash course unit conversion:

There is 2 unit sytems, SI and LBM units. The conversion factors are the base units themselves, for example you can say meter=5, kg=3, second=4, and then convert all units derived from the three base units.

By default, there is the set_m_kg_s(const float x, const float u, const float rho, const float si_x, const float si_u, const float si_rho) function that calculates the 3 base units from length, velocity and density in both unit systems:

    const float si_L = 18.0f; // length in [m]
    const float si_u = 300.0f*0.5144f; // velocity in [m/s]
    const float si_nu = 1.48E-5f; // air kinematic shear viscosity in [m^2/s]
    const float si_rho = 1.225f; // air density in [kg/m^3]

    const float lbm_L = 512.0f; // length in LBM units
    const float lbm_u = 0.1f; // velocity in LBM units
    float lbm_nu; // kinematic shear viscosity in LBM units (calcualte this later)
    const float lbm_rho = 1.0f; // density in LBM units (always 1)

    set_m_kg_s(lbm_L, lbm_u, lbm_rho, si_L, si_u, si_rho); // find base units [m], [kg], [s]

    lbm_nu = units.nu(si_nu); // now that base units [m], [kg], [s] are known, any other quantity can be converted

However you don't have to use a length, velocity and density. Any three linearly independent quantities will work, such as length, kinematic shear viscosity and density. Once you have the base units, you can feed them to units.set_m_kg_s(const float m, const float kg, const float s):

    const float si_L = 18.0f; // length in [m]
    const float si_u = 300.0f*0.5144f; // velocity in [m/s]
    const float si_nu = 1.48E-5f; // air kinematic shear viscosity in [m^2/s]
    const float si_rho = 1.225f; // air density in [kg/m^3]

    const float lbm_L = 512.0f; // length in LBM units
    float lbm_u; // velocity in LBM units (calcualte this later)
    const float lbm_nu = 1.0f/6.0f; // kinematic shear viscosity in LBM units, such that 1=tau=3*nu+0.5
    const float lbm_rho = 1.0f; // density in LBM units (always 1)

    { // find base units [m], [kg], [s]
        const float m = si_L/lbm_L; // length si_L = lbm_L*[m] --> [m] = si_L/lbm_L
        const float kg = si_rho/lbm_rho*cb(m); // density si_rho = lbm_rho*[kg/m^3] --> [kg] = si_rho/lbm_rho*[m^3]
        const float s = lbm_nu/si_nu*sq(m); // kinematic shear viscosity si_nu = lbm_nu*[m^2/s] --> [s] = lbm_nu/si_nu*[m^2]
        units.set_m_kg_s(m, kg, s);
    }

    lbm_u = units.u(si_u); // now that base units [m], [kg], [s] are known, any other quantity can be converted
    print_info("for tau=1, lbm_u="+to_string(lbm_u, 8u));

I am getting values of tens of thousands for lbm_u now, and there's no reasonable change I can make to lbm_nu to bring lbm_u into the range of 0.0001 - 0.2.

lbm_u = 30547.933:

```cpp void main_setup() { // setup the volume and objects const uint lbm_L = 512u; // size of test volume const uint Nx = to_uint(lbm_L * 1.1f); // adjust size of volume x axis const uint Ny = to_uint(lbm_L * 1.1f); // adjust size of volume y axis const uint Nz = to_uint(lbm_L * 0.4f); // adjust size of volume z axis const float size = 0.75f * lbm_L; // scaling of object compared with size of the volume // setup aircraft parameters const float knots = 150.0f; // speed through air const float AoA = 0.0f; // angle of attack (°), rotates object in field on x-axis const float3 center = float3(0.5f * Nx, 0.45f * Ny, 0.55f * Nz); // offset the aircraft position const float3x3 rotation = float3x3(float3(1, 0, 0), radians(-AoA)); // set the aircraft rotation, check the specified axis // setting SI units const float si_L = 18.0f; // characteristic length const float si_u = knots * 0.5144f; // convert airspeed to SI units const float si_rho = 1.225f; // air density (~sea level) const float si_nu = 1.48E-5f; // kinematic shear viscosity of air (~sea level) // setting LBM units float lbm_u; const float lbm_rho = 1.0f; // density in LBM units always has to be 1 const float lbm_nu = 1.0f/6.0f; // kinematic shear viscosity in LBM units // set unit conversion factors between SI and LBM units const float m = si_L/lbm_L; // length si_L = lbm_L*[m] --> [m] = si_L/lbm_L const float kg = si_rho/lbm_rho*cb(m); // density si_rho = lbm_rho*[kg/m^3] --> [kg] = si_rho/lbm_rho*[m^3] const float s = lbm_nu/si_nu*sq(m); // kinematic shear viscosity si_nu = lbm_nu*[m^2/s] --> [s] = lbm_nu/si_nu*[m^2] units.set_m_kg_s(m, kg, s); lbm_u = units.u(si_u); // now that base units [m], [kg], [s] are known, any other quantity can be converted print_info("desired lbm_u = " + to_string(lbm_u, 8u)); // create LBM object LBM lbm(Nx, Ny, Nz, lbm_nu); // setup mesh and set flags Mesh* mesh = read_stl(get_exe_path() + "../stl/aircraft.stl", lbm.size(), center, rotation, size); lbm.voxelize_mesh_on_device(mesh); lbm.flags.read_from_device(); const uint N = lbm.get_N(); for (uint n = 0ull, x = 0u, y = 0u, z = 0u; n < N; n++, lbm.coordinates(n, x, y, z)) { if (lbm.flags[n] != TYPE_S) lbm.u.y[n] = lbm_u; if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) lbm.flags[n] = TYPE_E; } std::ofstream fout(get_exe_path() + "results.txt"); if (!fout) { std::cerr << "Could not open file." << std::endl; } while (lbm.get_t() < 30000u) { lbm.run(100u); // calculate the force on the object lbm.calculate_force_on_boundaries(); lbm.F.read_from_device(); const float3 force = lbm.calculate_force_on_object(TYPE_S); // write the calculated force to file fout << to_string(units.si_F(force.y)) << std::endl; print_info(to_string(units.si_F(force.y))); } // close the output file fout.close(); } ```

I think that's time for me to leave this software alone now, I'm just not familiar enough with CFD or LBM to create any rational simulations. Thanks for all the time helping me out with various aspects of the software, it's much appreciated.

@SirWixy this is expected for lbm_nu=1/6. To get lbm_u in a reasonable range, lbm_nu has to be tiny, like 0.00001. lbm_nu=1/6 would be ideal but is not possible in most cases; I just wanted to show you why. There is a very wide range of lbm_nu that works reasonably accurate, even if it has to be tiny. Generally it's simpler to set/tune lbm_u, and calculate lbm_nu per unit conversion, as that is more flexible.

LBM has trouble for too large Mach numbers. In your initial post, Ma=0.46; air speed is too large. For high Mach numbers, the air compresses significantly in front of the airplane. At even higher Mach, shock waves form. But LBM can only handle the regime where the air does not compress significantly, which generally is Ma<~0.3. For very high airspeed, LBM is not the right tool unfortunately.

ProjectPhysX / FluidX3D

Help with force/drag calculation for maximum airspeed calculation #55