Qustion about extracting the lift moment of a thin flat

[Winxingtf]

Dear all,

Recently, I want to test a thin flat plate problem both in 2d and 3d in PyFR v1.12.1. The parameters I need are the lift force Fl and lift moment M. The plugin [soln-plugin-fluidforce-name] is helpful, but seems to contain only the forces in two direction. There are two directions I can follow. One is to add the third component to the plugin in solver stage. I’ve checked the fluidforce.py under the plugins folder, which is quite hard for me to understand. I also looked over the Ansys Fluent Theory Guide, which details the calculation of the total force component and total moment:

Fa = a·Fp + a·Fv
Ma = rAB×Fp + rAB×Fv

Where a means specified force vector, Fp and Fv means pressure and viscous force vector, respectively. rAB means moment vector, point from moment center A to force origin point B. So if I want to get the moment, I need a position vector and a force vector. How can I extract the vectors and accomplish it in the plugin? Or in a cumbersome way using the [soln-plugin-sampler] instead, which would burden my work enormously.

The other way is to obtain the moment term with paraview filters in post-process stage. With the filter such as ExtractBlock、ExtractSurface、GenerateSurfaceNormals, etc. several tutorials can be found on YouTube. But in most cases, the output files are in .case format from OpenFOAM. It would be confusing to do the same operation since many filters are gray to pick on both 2d and 3d cases.

Really hope someone can offer me any suggestion, and thanks a lot in advance for your kind help!

Regards, Thatcher

[GonzaloS]

Hi @Winxingtf ,

The fluidforce plugin already computes the two or three-dimensional pressure and viscous forces separately. We cannot directly recover the moments from the forces, since we need to assess the integral M = \int_S \vect{r} \cross (d F_p + d F_v) to compute the moments.

If you feel comfortable with VTK, you could try to interpolate your solution to the STL/STEP description of your surface to interest to then compute the forces/moments. This could also be accomplished by extracting the surfaces and the clipping/selecting the cells/faces of your surface. However, the latter approach might not work for certain geometries.

Implementing the computation of moments directly in the plugin is rather straightforward (i.e. cross product of the the forces at each flux point by the vector r and then assessing the integral using the same quadratures used to compute the forces). You would need to add an extra input token to define the reference point.

[Winxingtf]

Hi Gonzalo S,

Thanks for your explanation and suggestions. I agree with you that I need to do some work with either paraview or plugins-fluidforce stuff. However, both are not easy for me. Since I need to obtain the time-averaged force, I would like to try the second way, which would be more straightforward as you suggested.

After groping through the code, I thought I just need to modify those three lines and add a new one in fluidforce.py. (Assume the reference point is set to (0,0), then the components of vector r would be its coordinates.)

Line123: f = np.zeros(2*ndims if self._viscous else ndims) #Define the force vector
Line147: f[:ndims] += np.einsum('i...,ij,jik', qwts, p, norms) #Do the quadrature to get Fp
Line174: f[ndims:] += np.einsum('i...,klij,jil', qwts, vis, norms) #Do the quadrature to get Fv

Also, I’m supposed to add f[] += np.einsum('',f,r,norms) to get the required M. So now the problem is where can I get the coordinates of each solution points and how to combine the two force components to total force for integration. Actually, I’ve just switched from commercial software to an open-source one, so there are a lot for me to learn. Please correct me if I misconstrue somewhere and hope you can elaborate a bit more.

Regards, Thatcher

[GonzaloS]

The problem is a bit more complex since you need to do the cross product. Something along these lines should work, although this is untested code.

# Get flux points
ele_ploc_fpts = [e.ploc_at_np('fpts') for e in eles]

# TODO: select the flux points belonging to the BC
# see the eidx and fidx loops in fluidforce
bc_fpts = ...

# Moments
mom = np.zeros(2*ndims if self._viscous else ndims)

# Vector from flux points to moment origin
# WARNING: may need to use some numpy magic here
# to compute r (np.newaxis may help)
r = ele_ploc_fpts - origin

# Pressure moments
mom[:ndims] += np.einsum('i...,ij,jik', qwts, p, np.cross(r, norms))

# Viscous moments
...

BTW, running your simulation with a Python debugger (VSCode + Pylance for example) will help you better understand the array indexing used to store data in PyFR.

[Winxingtf]

Hi Gonzalo S,

Thanks for your prompt and kind reply! This would be a good start for me to work on the modification and debugging. I installed PyFR directly using pip install, avoiding much detailed and complicated compilation. So maybe I’ll need to compile the solver from source and figure out how the solver works in these days.

Regards, Thatcher

[fdw]

Check out the source code from git (or download a .tar.gz/.zip file and extract). Then change into this directory and run pip install -e . This will install PyFR but does so in such a way that any modifications you make to the source are automatically visible when you run pyfr (so there is no need to reinstall between changes).

Regards, Freddie.

[Winxingtf]

Hi GonzaloS and fdw,

Sorry for my late response, I was building a new compiling environment these days, and it can finally run with CUDA backend. I tried to debug with the first line as following

@GonzaloS

ele_ploc_fpts = [e.ploc_at_np('fpts') for e in eles]
print(ele_ploc_fpts)

These were included in the for loop: for etype, eidx, fidx, flags in mesh[bc].astype('U4,i4,i1,i2'): for clarification of eles. But the error came up like the following

Traceback (most recent call last):
  File "G:\software\envs\pyfr_tf\Scripts\pyfr-script.py", line 33, in <module>
    sys.exit(load_entry_point('pyfr', 'console_scripts', 'pyfr')())
  File "g:\pyfr\pyfr\pyfr\__main__.py", line 117, in main
    args.process(args)
  File "g:\pyfr\pyfr\pyfr\__main__.py", line 269, in process_restart
    _process_common(args, mesh, soln, cfg)
  File "g:\pyfr\pyfr\pyfr\__main__.py", line 232, in _process_common
    solver = get_solver(backend, rallocs, mesh, soln, cfg)
  File "g:\pyfr\pyfr\pyfr\solvers\__init__.py", line 16, in get_solver
    return get_integrator(backend, systemcls, rallocs, mesh, initsoln, cfg)
  File "g:\pyfr\pyfr\pyfr\integrators\__init__.py", line 36, in get_integrator
    return integrator(backend, systemcls, rallocs, mesh, initsoln, cfg)
  File "g:\pyfr\pyfr\pyfr\integrators\dual\phys\controllers.py", line 8, in __init__
    super().__init__(*args, **kwargs)
  File "g:\pyfr\pyfr\pyfr\integrators\dual\phys\base.py", line 23, in __init__
    self.completed_step_handlers = proxylist(self._get_plugins())
  File "g:\pyfr\pyfr\pyfr\integrators\base.py", line 75, in _get_plugins
    plugins.append(get_plugin(name, self, cfgsect, suffix))
  File "g:\pyfr\pyfr\pyfr\plugins\__init__.py", line 17, in get_plugin
    return subclass_where(BasePlugin, name=name)(*args, **kwargs)
  File "g:\pyfr\pyfr\pyfr\plugins\fluidforce.py", line 78, in __init__
    ele_ploc_fpts = [e.ploc_at_np('fpts') for e in eles]
TypeError: 'ACNavierStokesElements' object is not iterable

So I went to check the solver/acnavstokes/elements.py and found nothing relates to fpts. Meanwhile, I could not truly understand the definition of bc_fpts. To calculate the total force, in my view, the location and the flux value of the solution points belonging to the boundary were needed primly. So maybe just a simple exchange of the norms with np.cross(r, norms)? Correct me if I was wrong.

Regards, Thatcher

[WillT]

I am not an expert on this part of the code but I don’t think that the physical location for the flux points is stored. Hence why this isn’t working. The function ploc_at_np isn’t defined in the solver specific elements class, but in the base class. See pyfr/solver/base/elements.py. From there I think the function you want is plocfpts().

FYI, if you don’t already use it, the terminal tool ack is very useful for searching source code.

[GonzaloS]

Hi @Winxingtf could you please try Add moments computation to fluidforce by GonzaloSaez · Pull Request #239 · PyFR/PyFR · GitHub and see if it works for 2D and 3D cases? I’ll try to clean up a bit and fix some possible bugs soon.

[Winxingtf]

Hi @GonzaloS

Thanks for the PR and it’s really helpful. I had tested the PR in PyFR.1.12.2 both in 2d and 3d simulation. There existed two bugs, and I tried to solve the first one by making some modification.

1. Output .csv files only contains header. The following lines worked.

# Write 
if self.morigin:
      print(intg.tcurr, *f, *m, sep=',', file=self.outf)
else:
      rint(intg.tcurr, *f, sep=',', file=self.outf)

2. When I ran a 3d simulation, an error came out as following

Traceback (most recent call last):
  File "G:\software\envs\pyfr_tf\Scripts\pyfr-script.py", line 33, in <module>
    sys.exit(load_entry_point('pyfr', 'console_scripts', 'pyfr')())
  File "g:\pyfr\pyfr\pyfr\__main__.py", line 117, in main
    args.process(args)
  File "g:\pyfr\pyfr\pyfr\__main__.py", line 269, in process_restart
    _process_common(args, mesh, soln, cfg)
  File "g:\pyfr\pyfr\pyfr\__main__.py", line 243, in _process_common
    solver.run()
  File "g:\pyfr\pyfr\pyfr\integrators\base.py", line 103, in run
    self.advance_to(t)
  File "g:\pyfr\pyfr\pyfr\integrators\dual\phys\controllers.py", line 47, in advance_to
    self._accept_step(self.pseudointegrator._idxcurr)
  File "g:\pyfr\pyfr\pyfr\integrators\dual\phys\controllers.py", line 29, in _accept_step
    self.completed_step_handlers(self)
  File "g:\pyfr\pyfr\pyfr\util.py", line 53, in __call__
    return proxylist(x(*args, **kwargs) for x in self)
  File "g:\pyfr\pyfr\pyfr\util.py", line 53, in <genexpr>
    return proxylist(x(*args, **kwargs) for x in self)
  File "g:\pyfr\pyfr\pyfr\plugins\fluidforce.py", line 221, in __call__
    m[self._mcomponents:] += np.einsum(
  File "<__array_function__ internals>", line 5, in einsum
  File "G:\software\envs\pyfr_tf\lib\site-packages\numpy\core\einsumfunc.py", line 1359, in einsum
    return c_einsum(*operands, **kwargs)
ValueError: operand has more dimensions than subscripts given in einstein sum, but no '...' ellipsis provided to broadcast the extra dimensions.

Maybe an if statement is required? All in all, the function was quite practical, and the calculation went all well. I was planning to keep testing to see if the moment was calculated correctly.

Regards, Thatcher

[GonzaloS]

I think that I fixed the first error in my PR yesterday. Could please try again with the current PR commit?

The second may be an issue with np.cross. Could you send a me very small 3D case to test this please?

Regards,

Gonzalo

[Winxingtf]

I revisited the Github and saw two PR posted, I might mistake the two, and I’ll retry with the newer one. Also, there is a 3d square cylinder cases in the attachment, about 31,000 elements.

Regards, Thatcher
The 3d square cylinder

[GonzaloS]

I just pushed to the PR a bugfix for 3D computations. I see NaNs in the forces after 1 step so something must be wrong in the config file I suppose. Let me know if the moments that are output sound reasonable

[Winxingtf]

Hi Gonzalo,

Sorry for the NaNs caused by the config file. I just tested the 3d case setting the dt and dtau to 1/10 of the origin, and it went all well. No errors broke the process.

I was running another cases to test the correction of the moment coefficient, up to now, the time-average value was almost the same as the experimental value.

Thanks again for the helpful PR.

Regards, Thatcher

[GonzaloS]

Awesome, that’s great news! Let’s see if we can get this PR merged soon.

Thanks for testing this! Regards,

Gonzalo