Hello, everyone. I’m currently simulating the rotation of a Inclined cylinder. I want to rotate it along the cylinder’s axis, but the axis is tilted at an angle of 5°. I’ve decomposed the rotation direction along the X and Y axes, but I’ve found that it doesn’t end up rotating around the axis. Why is that?
GT.zip (249.7 KB)
Try setting this up as a rotating boundary rather than a rotating internal surface:
I believe there is an error in the computation of rotation for internal surfaces. File get_stl_data.f, lines 1597-1605:
IS_WALL_ANGLE(IS_ID,:) = IS_WALL_ANGLE(IS_ID,:) + IS_WALL_OMEGA(IS_ID,:)*DTSOLID
! OMEGA_CONV_FACT converts the angle to radians.
! The routines building the rotation matrices needs degrees.
Theta_deg(:) = IS_WALL_ANGLE(IS_ID,:) * OMEGA_CONV_FACT * 180.0/PI
CALL BUILD_X_ROTATION_MATRIX(Theta_deg(1), Rx)
CALL BUILD_Y_ROTATION_MATRIX(Theta_deg(2), Ry)
CALL BUILD_Z_ROTATION_MATRIX(Theta_deg(3), Rz)
R = MATMUL(Rz,MATMUL(Ry,Rx))
This is accumulating the Euler angles independently and computing a product matrix. The problem is that Rx and Ry do not commute so this is not actually rotation about a fixed axis! (The composition of two independently-growing Euler angles traces a path on SO(3) that curves away from the one-parameter subgroup exp(t[ω]×). You can see it’s wrong at a specific point: take t such that θx = 90°. Then Rx maps the y-axis to z. Ry has meanwhile rotated by −2.615/29.886·90 ≈ −7.9°. The product rotates points in a way whose equivalent single axis is no longer (0.9962, −0.0872, 0) — the tilt has effectively migrated out of the x–y plane because Rx carried the y-component up toward z.)
However the rotating-BC code at line 1647 of the same file seems to do this correctly:
OP = CENTER(:) - BC_WALL_ROT_CENTER(BC_ID, :)
AVG_VEL(:) = BC_WALL_VEL(BC_ID, :) + CROSS(OMEGA, OP)
AVG_VEL = v + ω × OP is correct, R = Rz·Ry·Rx is broken.
Your cylinder should probably be a BC rather than IC, but if you really need a rotating IC you can patch the code:
Copy get_stl_data.f to the project directory for editing and replace lines 1601-1605 with:
! omega in rad/s
omega(:) = IS_WALL_OMEGA(IS_ID,:) * OMEGA_CONV_FACT
! total angle since start, about the FIXED axis omega
! (accumulate the scalar angle, or use omega*time)
theta = sqrt(dot_product(omega,omega)) * <accumulated time>
if (theta > 0) then
n(:) = omega(:) / sqrt(dot_product(omega,omega))
R = build_rodrigues(n, theta) ! I + sin θ [n]_x + (1-cos θ)[n]_x^2
else
R = identity
end if
After I get a chance to test it, this fix will be included in the next MFiX release.
– Charles
Thank you very much, Charles. If I change the rotation plane to “boundary,” how does MFIX handle this? Does it rotate the STL file directly, or does it apply an additional velocity?
I believe rotating IS rotates the STL, while rotating BC applies additional velocity. You can look at get_stl_data.f for details
I get it, Charles! So if I want to simulate geometry like this involving axial asymmetry, should I still use an IS? I tried using the “Boundary” rotation option for this wall with protrusions, but found that the wall wasn’t actually rotated.
Ah, if the boundary has protrusions the rotating-BC approach probably won’t work. Try patching the IS code as described above
PS. I’m off for the rest of the week, I’ll follow up on this next week.
Another thing you can try is to make the cylinder horizontal and tilt the gravity vector
The moving IS or rotating BC are simple work arounds that have only been tested and meant to be used for simple rotations along one on the major axes (x, y or z). The rotating BC doesn’t actually move the wall, it only applies a tangential velocity. This is meant to be used for perfectly cylindrical walls (smooth walls, no protrusions). Rotating the IS moves the STL geometry, but rotation should only be applied along one major direction. The fluid doesn’t see the BC or IS rotation, so this should only be used for granular flows.
I think you could achieve what you want by keeping the cylinder’s axis along the x-direction, apply the rotation only along the x-axis, but tilt the gravity by 5 degrees. Then you can tilt particles and geometry by 5 degrees when you do the visualization (say in Paraview).
Other notes: You have small particles and a large geometry, you should expect the simulation to take a very long time if you have many particles. It may be a good idea to upgrade mfix (you are using the 24.4 version).
Thank you very much, Charles! I think I understand now.
Thank you very much, Jeff. I’m wondering if upgrading to the new version will help improve computation speed?
You might also consider trying this case in MFIX-Exa, especially if you have access to a GPU
I tested this case with my new code and discovered a small issue.
If you do decide to upgrade to MFIX 26.1.1, due to a change in the “simpleeval” module, the formula containing “pi” throws an error (pi not properly pre-defined). To fix this, replace the file simpleeval_wrapper.py with the file attached below.
simpleeval_wrapper.py (1.7 KB)
(This file goes in $CONDA_PREFIX/lib/python3.13/site-packages/mfixgui/tools/ on Linux or
%CONDA_PREFIX%\Lib\site-packages\mfixgui\tools on Windows)
Note that we will be issuing a 26.1.2 release soon which will fix this and a few other bugs)
The patch I uploaded was incomplete. If you want to do a rotation around a non-coordinate axis, put this file in your project directory and build a custom solver:
get_stl_data.f (91.6 KB)
This uses a Rodrigues matrix to rotate about the specified axis rather than a product of elementary X,Y,Z rotations which leads to drift.
I tested with your case and it works:
Thank you very much, Charles! I’ll give your suggestion a try.
Note that the new rotation code has been added to MFiX version 26.1.2 (just released today)
