Flap force simulation

Hi Laurence & all,

so it took a couple days more with a few failed test runs, and on the core subject of stall detection i must admit defeat;

SnappyHexMesh isn't so good at wrapping boundary layers around coarsely resolved small radii, and if i increase the resolution on the foil surface, at some point the cells height becomes so low that it can't create thick enough boundary layers within them to reach a high enough y+ for efficient solving. So the approach that was good enough for my softer, tidal-turbine style foils with twice the chord just keeps generating collapsing boundary layers at the leading edge for this tiny&thin section with its rather sharp leading edge. You'll see in the pic below how the boundary layers collapse at the LE, and that the y+ on much of the suction side is below 30, which should still work, but differently than the parts that actually rely on the wall function ... not so reassuring, and not exactly fast to solve either. It's probably good enough for L/D calculations at decent angles, but with the mess at the leading edge, i wouldn't trust it for stall inception.

Of course, that doesn't mean it can't be done (anything can with more manual work, more structured meshing and more calculation time. I don't have a cfMesh+ license anymore but that could certainly have tackled it, and it would even be scriptable - but require some babysitting nevertheless). But my aim here was to see if the free and automated workflow does the job.

Maybe it can still be useful for cross-code-validation at the calmer AoA's.
Here's what i got at the 'neutral' angle Laurence gave me:
All in N, all for a half T-Foil !

At 0.5s (flow still developing):
Total = Foil, Winglet, Bulb, Strut
Drag Sideforce Lift m/s
pressure 17.86 183.8 330.6 10.00
viscous 3.60 L/D kn
sum 21.5 15.4 19.4

Wing = Foil, Winglet
Drag Sideforce Lift m/s
pressure 14.65 2.6 331.1 10.00
viscous 2.14 L/D kn
sum 16.8 19.7 19.4

Foil without winglet
Drag Sideforce Lift m/s
pressure 14.06 7.3 322.9 10.00
viscous 2.06 L/D kn
sum 16.1 20.0 19.4
-----------

and at 1.0s:
(notice there is almost no lift left, which may be correct for the neutral AoA)

Total = Foil, Winglet, Bulb, Strut
Drag Sideforce Lift m/s
pressure 23.48 180.8 91.4 10.00
viscous 2.53 198% L/D kn
sum 26.0 #REF! 3.5 19.4

Wing = Foil, Winglet
Drag Sideforce Lift m/s
pressure 18.24 4.9 97.8 10.00
viscous 1.61 5% L/D kn
sum 19.8 #REF! 4.9 19.4

Foil without winglet
Drag Sideforce Lift m/s
pressure 17.68 7.8 92.7 10.00
viscous 1.53 8% L/D kn
sum 19.2 #REF! 4.8 19.4

Here's how it looks after 0.5s of solving:

"Surface LIC" on the inner slice still looks unrealistically smooth over the whole chord:


And this is after 1.0s:

Surface LIC on the inner slice starts to show a bit of separation at the trailing edge:


The workflow in the .zip below was primarily intended to investigate semi-foiling bigger boats. It uses blockMesh to create a background mesh with near-equilateral cubes around the hull geometry (for snappyHexMesh and mesh motion to work properly), while keeping cell count reasonably low by adapting the shape to the diagonality of the wake and grading the sizes towards the far field. It does a first run on a coarser mesh (about 3M cells here) without boundary layers to get an initial solution quickly, then maps it on a finer mesh with boundary layers (around 7M cells here) in a subdirectory called "BLrun". Build-up of the forces is captured in postProcessing/forcesTotal/0.1/forces.dat (or forcesWinglet etc) for each run. A .pvsm to look at it via paraview is also provided. (It can also handle mesh motion for sink and trim, just edit the constant/dynamicMeshDict for CG and rotation point and set "foamDictionary constant/dynamicMeshDict -entry mover/accelerationRelaxation -set 0.0" to "0.1" in the last run in the Allrun script.)

Basic usage would be:
- install openFoam 13, i use it on Ubuntu 24.04
- unzip the attached archive into a folder in your working directory
- make the "Allrun", "Allclean" "BoundaryLayerRemesh" "Restart" and "tail" scripts executable using chmod +x filename
- edit the "Allrun" script to adapt
"foamDictionary system/decomposeParDict -entry numberOfSubdomains -set 8" to your physical core count (NOT "threads" or hyperthreading count. Hyperthreading is only marketing BS as far as CFD/HPC is concerned). (If you have only similar cores, set it to 1 less than your actual core count, else the system might become unresponsive. If you have power cores and economy cores like Laurence's i9, set it to the number of power cores.)
- prepare your geometry: nice, watertight and smooth high resolution STL (when in doubt and on Rhino, shrink-wrap it to 0.0002m resolution), bow towards positive X, positive Y = side used in the sim, positive Z = up. The geometry shoud stick out at least very slightly towards the negative Y.
- replace the geometry files called dummy*.stl in constant/geometry with your geometry, names ending on:
***Foil.stl
***Strut.stl
***Winglet.stl
***Bulb.stl
- edit speed in "0/U" (it's in m/s): UMean 10.0;
- In a separate terminal, or using the "screen" utility, cd into your case directory and launch the Allrun script : ./Allrun
- Check on the meshing and calculations via another terminal.
- to look what your computer is doing and how much memory is used: top
- to look at the beginning of a log file:nano log.snappyHexMesh
- to look at the end of a log file: tail -n 5000 log.snappyHexMesh
- to look at a log file mid-flight: tail -f log.snappyHexMesh
- once the meshing went well, you can launch the "tail" script to see a less cluttered version of the log.foamRun
...watch out for the Courant numbers - if they go up above like 20 and keep raising,
problems are usually past a little spray and your calculation will blow up shortly.
- once the solving gets past the short initial runs, forces can be read up in postProcessing/forcesTotal/0.1/forces.dat (or forcesWinglet etc)
- when the initial run without boundary layers is complete, switch into the BLrun subdirectory where the higher resolution version will calculate itself.
- to visualise in paraview, open a new terminal, get into your case folder (and then into the BLrun folder if we're already that far). Type "paraview" & Enter in the terminal, paraview should launch. Hit File/LoadState, it's supposed to open directly into your case folder. Click on "thisCase.pvsm", then in the "Load State Options" window, choose "Search Files under Specified Directory", point it to your case folder (or the corresponding BLrun subdirectory if you're already that far), and tick the "Only use files in data directory" option before hitting Ok. Hopefully, the visualisation will open as above.

Water density and kinematic Viscosity can be adjusted in constant/physicalProperties.water
Air density and kinematic Viscosity can be adjusted in constant/physicalProperties.air
Surface roughness can be adjusted in 0/nut , parameter Ks for each of the patches. (As far as i understood, it is a "sand grain height" expressed in meters)

Performance was about 10h to reach 0.5s on the BLrun (second run with boundary layers) on my twin Epyc "Genoa" generation server with all 24 DDR-4 slots populated. Perhaps 16hours or so to solve to 1s (that's alredy several chord lengths in this case, the flow should be well established by then).

Feedback welcome !
Would be nice to include the files with corrections in case.
Cheers!

 

Maximum flap deflection on my Mach2 foils is less than half of that. I don't know for sure about more modern designs, but I seriously doubt if they are much more.

 

I'm not sure about wether the mesh is appropiate or not. In my opinion it should stay refined to a greater distance of the foil, maybe even more refined.

Use a 2D setup for a foil for which benchmark data are available NACA. Then evaluate the section of @laurencet 's foil in that 2D setup. There is in first instance really no need for hours of running a simulation which basically says notthing.

 

Sure it would be neater with more resolution and smoother boundary layers (typically i'd use 3 or 4 BL's with a growth ratio of 1.2- 1.4). I already wrote why i don't get that with the existing workflow for such a thin tiny foil. I'm also avoiding refinement boxes in openFoam because they need manual setup for each geometry, and because instead of cutting a cell into 8 smaller cubes, refineMesh and refineHexMesh often create pyramids that may create problems later on when they get squashed in mesh motion. On the larger tidal-turbine-style sections i did some grid independence testing and found that the usual high refinement all around the foil up to something like half a chord length with fairly progressive sizing of the "buffer cells" didn't matter too much in practice as long as the first boundary layer or two were reasonably smooth and operating within y+ 50-200. On Laurence's foil i can't validate anything because we have no real-world data. I'd be happy to throw in the Hungry Beaver if someone has the patience to model it

I'd actually be totally happy to use that if someone knowledgable both about practical hydrofoil usage and openFoam puts in the necessary day or two of work to tune a workflow and validate, and then actually shares it in a usable format. It's just not something i have laying around. (I did 2D sections a while ago on openFoam, meshing with cfMesh+ and HexPress, but that involves licenses, and worse, manual work each time, so it doesn't scale). The necessary tools have been freely available for more than a decade and that kind of calculation can be done on any laptop nowadays. Yet noone is doing it because it takes a couple months to get familiar with openFoam initially, and the few folks who did, and went through the trouble of acctually applying it to their field, are typically busy engineering, tweaking, and monetizing.

(Multiphase/marine is a tiny, poor and complicated niche of CFD. Contrary to the mainstream aero business of CFD there are few big money players like military, car or airplane companies, F1 etc to pay for developments. It's also politically and technically complicated for commercial CFD vendors to charge vastly lower license fees to our segment (where a naval architect would maybe spend 10-20% of their time working with it) than for mainstream applications done by specialists in larger teams, where they can simply stress out how their black box workflows save the client one or two 100k/year engineer salaries.)

(Perhaps Naca 63-412 could be a starting point for validation Foiler Design https://www.boatdesign.net/threads/foiler-design.2447/page-6#post-17400 and NACA 63-412 AIRFOIL (n63412-il) http://airfoiltools.com/airfoil/details?airfoil=n63412-il , but the other question is how to get the turbulence levels and smoothness about right in both openFoam and the panel code: XFoil Re and NCrit help https://www.boatdesign.net/threads/xfoil-re-and-ncrit-help.37693/#post-457528 )

From there, the next question would be to which extend the "infinite wing" 2D section data can be trusted to predict stall behaviour accurately. Personally i'd rather expect a stall to start around the tip, foil/winglet transition, bulb fairing, or at the sides of the deployed flap: All the areas where there is a 3D mess.

What i don't quite understand is the "no need to put in compute hours" vibe. Sure, preliminaries could be done quicker and more efficiently once you're set up for that ... but more systemic stuff not so much. Calculation power is a commodity for most tech students, with tens of thousands of HPC nodes rotting away in universities and public-funded cluster access programmes. Same for professionals, with cloud access at fairly reasonable prices nowadays, or powerful servers or workstations available below the price of an ageing second-hand car. Hopefully someone is going to resolve the automated meshing challenge for real one of these days - it should be a much easier task for AI than it was for more rigid algorithms, and quite a suitable one for massive parrallelization (GPU's). Even the energy consumption for an overnight CFD solver run (say 10kWh on older hardware) is below what most professionals burn on their commute to work, and mine runs pretty much entirely on solar energy. It's the work input that is the bottleneck, but that can be reduced to minutes per run with automation.

 

Thanks again—that analysis provided some great insight!
I'm sold on the results and will be installing OpenFOAM to continue the work myself.

The next physical step is swapping to a foil with a rounded leading edge, hopefully make the model easier to solve .

The intersection issues were eye-opening, I'm reluctant to dive into deep optimization until the FEA model is ready.
It's very easy to come up with a great CFD solution that won't survive the forces.

I'm looking forward to the initial build, set up a tow test with some strain guages and establish a performance baseline.


Is this just the induced drag, It might be worth putting a twist in the wing to reduce the lift at the tip.
Do you think the wing tip dihedral is making it worse?

 

That makes alot of sense, if the foil is stalling at 15 degrees, what's the point in operating there ..

 

Precisely! I think you will find that the desirable operating range gets ever smaller as the flap deflection angle increases. Here's a plot of XFOIL results that shows that trend pretty clearly. (The plot was made about 10 years ago as part of preliminary studies for a successful foil design project with US Moth superstar Matt Struble.)

 

Oh wow, i thought my 'slightly below confidence level' results would be rather discouraging

I'd encourage you to stay very sceptical then.
For example :

It might well be that the one you already have is more efficient than anything more rounded you can come up with. The published big-nosed sections often rely heavily on aft-loading, which, as you said, isn't trivial to build strong enough.

Wait, by intersection, do you mean the bulb-foil discontinuity in y+ ?
There is a cut between the two, i meshed them a bit differently by error (already corrected in the case attached above). Light blue and salmon near the bulb means y+ ~50-150, where the wall function works fine.
Deep blue means y+ more like below 30, where use of the wall function isn't valid, but openFoam's "kOmegaSST" implementation should revert to direct solving.
Both are supposed to be valid, and they do actually give a continuity of pressure results as can be seen on the big p-rgh (relative pressure) pic up left.
Mixing the two methods with hard transitions in a critical area isn't optimal though, and it would all be faster to solve if i had gotten a mesh with more uniform first layer thickness and shape between Y+ 30 and 300.
So so far i'd assume that nothing is wrong with your geometry in that area !

These streamlines are only there to show the inevitable wing tip vortex. Coloration is vertical velocity of the local flow(Uz), that takes just a click. 3D velocity magnitude isn't helpful for coloring because the 10m/s of Ux get so dominant that you don't see much. With some paraview filtering, you could certainly find a way to color them by velocity magnitude in the y/z plane.

Induced drag isn't going to be so straightforward to isolate in a single CFD run; even with your wing at a neutral angle (no lift), there is going to be some pressure drag and some viscous drag. As far as i understand, you'd need to subtract the drag figures at that neutral angle from the ones of another solver run with an AoA where the foil actually generates lift. Careful, you might get variations with immersion depth

I'll get back to you on PM; perhaps you'd like to have a closer look at the already solved results ? It's gonna be too large to transfer but can be done on Teamviewer remote control.