[Krao]

Dear OpenFOAM Users,

After successfully completing the 2 wing UAV propeller simulation using MRFSimpleFoam and pimpleDyMFOAM I am now trying to reduce the computation time by carrying out simulation with only 1 wing instead of 2 with MRF method. The procedure which I am following is illustrated below.

1) Generation of separate mesh for static domain and the dynamic domain
2) Merging the two regions by using mergeMeshes
3) Generation of required topoSets
4) Definition of appropriate Initial and boundary conditions
5) Simulation using MRFSimpleFoam

The above steps mentioned were well accomplished while simulating the complete propeller. But with single propeller and half of the domain my, simulation is crashing at the beginning. A close investigation in this direction revealed that the min patch weight at few regions are '0', which means there is no interaction between the static and quasi-rotating meshes. I have attached the mesh below, along with a sample definition of initial and boundary conditions. Also, I have attached the output file at the end.

AMI_Dynamic.png AMI_Static.png Dynamic_Domain.png Static_Domain.png



With my experience so-far with OpenFOAM, I think it is happening because of the marked region which is not present in the stationary domain. Please let me know how to overcome these problems. If at all I have missed any thread with respect to above mentioned problem, please do let me know. I thank the forum members for lending their hands to teach me OpenFOAM this far.

Thank you.

Log_Output

Create time

Create mesh for time = 0


SIMPLE: convergence criteria
field U tolerance 0.0001
field p tolerance 0.001
field "(k|epsilon|omega)" tolerance 0.001

Reading field p

Reading field U

Reading/calculating face flux field phi

AMI: Creating addressing and weights between 80320 source faces and 85398 target faces
AMI: Patch source sum(weights) min/max/average = 0.505603, 1.5061, 1.01059
AMI: Patch target sum(weights) min/max/average = 0, 1.23623, 0.991999
Selecting incompressible transport model Newtonian
Selecting turbulence model type RAS
Selecting RAS turbulence model kOmegaSST
Selecting patchDistMethod meshWave
bounding k, min: 0 max: 0.00375 average: 0.00375
bounding omega, min: 0 max: 3.375 average: 3.375
Creating MRF zone list from MRFProperties
creating MRF zone: MRF1
No finite volume options present


Starting time loop

forces forces:
Not including porosity effects
Time = 0.0001

smoothSolver: Solving for Ux, Initial residual = 1, Final residual = 4.42885e-12, No Iterations 5
smoothSolver: Solving for Uy, Initial residual = 1, Final residual = 4.63258e-12, No Iterations 5
smoothSolver: Solving for Uz, Initial residual = 1, Final residual = 2.84064e-12, No Iterations 5
[1] #0 Foam::error:rintStack(Foam::Ostream&)[2] #0 Foam::error:rintStack(Foam::Ostream&) at ??:?
at ??:?
[1] #1 [2] #1 Foam::sigFpe::sigHandler(int)Foam::sigFpe::sigHand ler(int) at ??:?
[1] #2 ? at ??:?
[2] #2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
[1] #3 Foam::divide(Foam::Field<double>&, double const&, Foam::UList<double> const&) in "/lib/x86_64-linux-gnu/libc.so.6"
[2] #3 Foam::divide(Foam::Field<double>&, double const&, Foam::UList<double> const&) at ??:?
[1] #4 at ??:?
[2] #4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:perator/<Foam::fvPatchField, Foam::volMesh>(Foam::dimensioned<double> const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&)Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:perator/<Foam::fvPatchField, Foam::volMesh>(Foam::dimensioned<double> const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
[1] #5 in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
[2] #5 ?? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
[2] #6 in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
[1] #6 ?? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
[2] #7 __libc_start_main in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
[1] #7 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[2] #8 in "/lib/x86_64-linux-gnu/libc.so.6"

FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 1 -1 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{
inlet
{
type fixedValue;
value uniform (0 0 0);
}
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (0 0 0);
}
leftSideWall
{
type slip;
}
reightSideWall
{
type slip;
}
back
{
type slip;
}
symmetry
{
type symmetry;
}


AMI1
{
type cyclicAMI;

}

AMI2
{
type cyclicAMI;

}

blades
{
type movingWallVelocity;
value uniform (0 0 0);
}



}

[Swift]

Hi Krao,

As we had discussed (Periodic cyclic AMI + MRF crashing on pressure), try using the tutorial case to get good weightings on your AMI. That should solve that problem. Let me know if that solves your weighting issue. It might be worth trying for the full domain. In that tutorial, the rotating domain and stationary domain are meshed at the same time, so its a little bit easier.

In addition I think you are not using symmetry incorrectly. Your domain is not symmetrical about a plane. Your geometry is, but the flow is not. Think about if your impeller rotates out of the symmetry plane. The symmetry plane will make it as if the impeller rotates into the domain at that point. That is probably not what your are simulating. You would want the impeller to come in at the top half, i.e. 180 degrees about the rotation axis. To model it correctly, instead of using symmetry you need to use cyclicAMI. Note the red:

 boundaryField
 {
    rotor_periodic_1 
    {
        type            cyclicAMI;
        inGroups        1(cyclicAMI);
        matchTolerance  0.1;
        transform       rotational;
        neighbourPatch  rotor_periodic_2;
        rotationAxis    (0 0 1);
        rotationCentre  (0 0 0);
        nFaces          1628;
        startFace       171370;
    }
 }

see https://www.cfdsupport.com/OpenFOAM-...t/node116.html
Thus you will make use of two pairs of cyclicAMI, one pair for 'transform rotational' and one pair for 'transform noOrdering'. So, in your last image, the face closest to the view. Split this face in half along the rotational axis. The upper half will be one of the rotational AMI faces and the lower half, the other. Essentially, you need to use a periodic symmetry. This might be worth googling to see how its done.

It is also possible to use MRF without cyclicAMI altogether and this generates pretty good results as well. I use this technique extensively. Note this models the full domain, but it is possible to use the periodic rotational symmetry.

In my snappyHexMesh, I still call the interface between the rotating domain and stationary domain AMI and in "refinementSurfaces" I use:

AMIBottom
        {
            level       (2 2);
            cellZone    rotatingDomainBottom;
            faceZone    rotatingDomainBottom;
            cellZoneInside  inside;
        }

You will need to add cyclicAMI for the rotational symmetry faces.

Thomas

[Krao]

Hi Thomas,

Thank you very much for the quick reply, I had forgot to mention my boundary file, which I have attached for the future reference. As you have mentioned I am using the symmetry in a wrong way. I will try the new simulation with the suggestions made by you and post my findings.

9
(
blades
{
type wall;
inGroups 1(wall);
nFaces 10110;
startFace 5982954;
}
back
{
type wall;
inGroups 1(wall);
nFaces 5401;
startFace 6015462;
}
leftSideWall
{
type wall;
inGroups 1(wall);
nFaces 3000;
startFace 6020863;
}
reightSideWall
{
type wall;
inGroups 1(wall);
nFaces 3027;
startFace 6023863;
}
outlet
{
type patch;
nFaces 1620;
startFace 6026890;
}
inlet
{
type patch;
nFaces 1603;
startFace 6028510;
}
symmetry
{
type symmetry;
inGroups 1(symmetry);
nFaces 56056;
startFace 6030113;
}
AMI1
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 80320;
startFace 6086169;
matchTolerance 0.0001;
transform noOrdering;
neighbourPatch AMI2;
method faceAreaWeightAMI;
}
AMI2
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 85398;
startFace 6166489;
matchTolerance 0.0001;
transform noOrdering;
neighbourPatch AMI1;
method faceAreaWeightAMI;
}
)

[Krao]

Hi Thomas,

As you have suggested, I revisited my mesh and improved its quality. For that I made use of surfaceAutoPatches. This has solved the problem of AMI weight going to zero. As suggested by you I have changed the boundary conditions of rotating patches which I am attaching for the future users. Only thing I did not understand is the wrong usage of symmetry patches! Which boundary condition, can I use at the center of the domain? Can I use the normal wall boundary condition at the center of the stationary domain?

11
(
blades
{
type wall;
inGroups 1(wall);
nFaces 9973;
startFace 5970050;
}
back
{
type wall;
inGroups 1(wall);
nFaces 5097;
startFace 5980023;
}
leftSideWall
{
type wall;
inGroups 1(wall);
nFaces 2744;
startFace 5985120;
}
reightSideWall
{
type wall;
inGroups 1(wall);
nFaces 2771;
startFace 5987864;
}
outlet
{
type patch;
nFaces 1456;
startFace 5990635;
}
inlet
{
type patch;
nFaces 1437;
startFace 5992091;
}
symmetry
{
type symmetry;
inGroups 1(symmetry);
nFaces 52678;
startFace 5993528;
}
AMI1
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 101643;
startFace 6046206;
matchTolerance 0.0001;
transform noOrdering;
neighbourPatch AMI2;
method faceAreaWeightAMI;
}
AMI2
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 84970;
startFace 6147849;
matchTolerance 0.1;
transform noOrdering;
neighbourPatch AMI1;
method faceAreaWeightAMI;
}
AMI3
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 13281;
startFace 6232819;
matchTolerance 0.1;
transform rotational;
neighbourPatch AMI4;
rotationAxis (0 0 1);
rotationCentre (0 0 0);
method faceAreaWeightAMI;
}
AMI4
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 13285;
startFace 6246100;
matchTolerance 0.0001;
transform rotational;
neighbourPatch AMI3;
rotationAxis (0 0 1);
rotationCentre (0 0 0);
method faceAreaWeightAMI;
}
)

[Swift]

Hi Krao,

Can you upload a sketch of your geometry and how you have named the patches? This might make it easier to help.

Thomas

[Krao]

Hi Thomas,

thank you once again for your reply. I have attached the image of my computational domain with the names of the patches. Also I have edited few of my boundary patches.

1) As you have pointed the wrong use of symmetry, I came to know that the usage of 'symmetry patch' as shown in the figure is improper. I think I need to split this patch (in gray) into two halves and maybe name it like AMI 5 & AMI 6. Before proceeding with that, I would like to know whether the AMI (rotation transformation) can be applied to a stationary domain, which is my bounding box(computation area).

2) I am using noOrdering for AMI 1 & AMI 2. This is the half cylinder portion which is the MRF region, where the propeller blade is present. For the face of the same half-cylinder region marked AMI 3 & AMI 4, I am using rotational transformation. This usage of patches are due to the fact that the half cylinder region in static and dynamic domain have had difference in geometry, as shown in the last two figures. The static region does not have any face and contains only half cylinder, because of the snapping effect. The dynamic region have the half cylinder with a face. Therefore I have split this face into 2 patches and given AMI with rotational transformation.

For meshing I am using cfMesh and for combining the static and dynamic regions I am using mergeMesh utility.

Domain_Nomenclature.jpg AMI_Static.png AMI_Dynamic.png

In the above figure the highlighted circular region demonstrates the MRF region. The relevant boundary file is as follows.

11
(
blades
{
type wall;
inGroups 1(wall);
nFaces 9973;
startFace 5970050;
}
back
{
type wall;
inGroups 1(wall);
nFaces 5097;
startFace 5980023;
}
leftSideWall
{
type wall;
inGroups 1(wall);
nFaces 2744;
startFace 5985120;
}
reightSideWall
{
type wall;
inGroups 1(wall);
nFaces 2769;
startFace 5987864;
}
outlet
{
type patch;
nFaces 1456;
startFace 5990633;
}
inlet
{
type patch;
nFaces 1439;
startFace 5992089;
}
symmetry
{
type symmetry;
inGroups 1(symmetry);
nFaces 52678;
startFace 5993528;
}
AMI1
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 101643;
startFace 6046206;
matchTolerance 0.0001;
transform noOrdering;
neighbourPatch AMI2;
method faceAreaWeightAMI;
}
AMI2
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 84970;
startFace 6147849;
matchTolerance 0.1;
transform noOrdering;
neighbourPatch AMI1;
method faceAreaWeightAMI;
}
AMI3
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 13281;
startFace 6232819;
matchTolerance 0.1;
transform rotational;
neighbourPatch AMI4;
rotationAxis (0 0 1);
rotationCentre (0 0 0);
method faceAreaWeightAMI;
}
AMI4
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 13285;
startFace 6246100;
matchTolerance 0.0001;
transform rotational;
neighbourPatch AMI3;
rotationAxis (0 0 1);
rotationCentre (0 0 0);
method faceAreaWeightAMI;
}
)

Please let me know the mistakes I am doing, thank you. I would be pleased to share more information for the above problem.

Regards,
Krao

[Swift]

Hi Krao,

1) As you have pointed the wrong use of symmetry, I came to know that the usage of 'symmetry patch' as shown in the figure is improper. I think I need to split this patch (in gray) into two halves and maybe name it like AMI 5 & AMI 6. Before proceeding with that, I would like to know whether the AMI (rotation transformation) can be applied to a stationary domain, which is my bounding box(computation area).

This is exactly what you need to do. You could also group what you would call AMI5 and AMI6 with AMI3 and AMI4. You will therefore have no symmetry plane.

Look at this post Periodic cyclic AMI + MRF crashing on pressure
I was modelling only a quarter of the domain, whereas you are modelling half. In my case, each side of the wedge was a rotational AMI (or according to your naming: AMI5 and AMI6). The only difference is that your rotational planes are actually in the same plane verses mine that were 90 degrees apart.

2) I am using noOrdering for AMI 1 & AMI 2. This is the half cylinder portion which is the MRF region, where the propeller blade is present. For the face of the same half-cylinder region marked AMI 3 & AMI 4, I am using rotational transformation. This usage of patches are due to the fact that the half cylinder region in static and dynamic domain have had difference in geometry, as shown in the last two figures. The static region does not have any face and contains only half cylinder, because of the snapping effect. The dynamic region have the half cylinder with a face. Therefore I have split this face into 2 patches and given AMI with rotational transformation.

Its a little hard to follow what you have done. The last two pictures you sent confused me. AMI1 and AMI2 should both look like image 70079. AMI1 and AMI2 should both have almost the same number of faces, but they appear to be quite different from your boundary file you attached. See red below:

AMI1
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 101643;
startFace 6046206;
matchTolerance 0.0001;
transform noOrdering;
neighbourPatch AMI2;
method faceAreaWeightAMI;
}
AMI2
{
type cyclicAMI;
inGroups 1(cyclicAMI);
nFaces 84970;
startFace 6147849;
matchTolerance 0.1;
transform noOrdering;
neighbourPatch AMI1;
method faceAreaWeightAMI;
}

The faces in image 70078 that aligns up with your symmetry plane should be split about the rotation axis into AMI3 and AMI4. These two patches should be the only patches in this region. Correct me if I am wrong, but does part of either AMI1 or AMI2 overlap with AMI3 or AMI4? If not, then you have it right.

The amount of faces in AMI3 and AMI4 appear to be more similar, I think you are OK there.

For meshing I am using cfMesh and for combining the static and dynamic regions I am using mergeMesh utility.

I have only briefly tried using cfMesh for creating the static and rotating domains and merging them. It will be interesting to see what AMI weights you get in the end.

This should get your simulation going. Let me know how it goes.

Regards,
Thomas

[Krao]

Hi Thomas,

Thank you very much for your time and clearing all my queries.

The faces in image 70078 that aligns up with your symmetry plane should be split about the rotation axis into AMI3 and AMI4. These two patches should be the only patches in this region. Correct me if I am wrong, but does part of either AMI1 or AMI2 overlap with AMI3 or AMI4? If not, then you have it right.

Yes the AMI 1 and AMI 2 will not overlap with AMI 3 and AMI 4.

I have only briefly tried using cfMesh for creating the static and rotating domains and merging them. It will be interesting to see what AMI weights you get in the end.

I am using cfMesh from last 2 yeras and used it for dynamic simulation with pimpleDyMFoam and it works really well. I should cross check the refinements.

[Krao]

Hi Thomas,

Greetings!

I have successfully completed my propeller simulation and I would like to discuss few results here.

The AMI weights of different patch-pairs are as follows

AMI: Creating addressing and weights between 8361 source faces and 6796 target faces
AMI: Patch source sum(weights) min/max/average = 0.997204, 1.0174, 1.00002
AMI: Patch target sum(weights) min/max/average = 0.986969, 1.00041, 0.999951
AMI: Creating addressing and weights between 2673 source faces and 2673 target faces
AMI: Patch source sum(weights) min/max/average = 0.739323, 1, 0.999619
AMI: Patch target sum(weights) min/max/average = 0.964354, 1, 0.99995
AMI: Creating addressing and weights between 22322 source faces and 21947 target faces
AMI: Patch source sum(weights) min/max/average = 0.999699, 1, 1
AMI: Patch target sum(weights) min/max/average = 0.99962, 1, 1

The results obtained for the half propeller simulation are as follows,

sum of forces:
pressure : (0.291062 3.70328 22.952)
viscous : (0.000151899 0.251535 -0.0315308)
porous : (0 0 0)
sum of moments:
pressure : (-0.0905238 5.44235 -0.793572)
viscous : (-0.00136577 -0.00609683 -0.0587384)
porous : (0 0 0)

To compare the above results with the full propeller simulation, the results obtained with full propeller are stated below

sum of forces:
pressure : (0.000162248 -0.00117233 45.3323)
viscous : (4.234e-05 1.75508e-05 -0.0776876)
porous : (0 0 0)
sum of moments:
pressure : (-0.000181585 0.00397288 -1.50671)
viscous : (-3.22819e-06 -1.70228e-05 -0.158295)
porous : (0 0 0)

As it can be observed, the results in the Z-direction for the single propeller are almost half of the results from the full propeller. Interestingly, with the single propeller simulation, I am getting force and moment in Y-direction, which should be zero. Apart from computational domain, I have not changed anything from the full propeller case. What would be the possible reason for this error, or should I just ignore the forces in Y-direction and consider only the results in Z-direction. Your insight would be invaluable.

Regards,
Krao

[Swift]

Hi Krao,

Congratulations on getting it to work! Well done.

Your AMI weightings look good.

As it can be observed, the results in the Z-direction for the single propeller are almost half of the results from the full propeller. Interestingly, with the single propeller simulation, I am getting force and moment in Y-direction, which should be zero. Apart from computational domain, I have not changed anything from the full propeller case. What would be the possible reason for this error, or should I just ignore the forces in Y-direction and consider only the results in Z-direction. Your insight would be invaluable.

Why should the force in the y direction be zero. What force do you expect is creating the moment about the y-y axis?


You don't need to be concerned about your results but you do need to understand where they come from, they are not errors and you can't ignore them.

Imagine you had modeled your full case, but instead of modelling your propeller as one part, you modeled it as two parts, one for each blade. You would then get force results for each blade that look similar to the half-domain case that you just modeled. Think about what forces and moments results you would get.

Firstly your moment about y-y axis. For each half impeller you would get a moment and the total would some to the total moment you have already calculated.

But now think about the force required to produce a moment about y-y. if your blade is in the x direction, the force for the one blade would be positive y, and the force on the other blade would be in the negative y direction. The sum of forces would be close to zero, which is what you are seeing in your complete case, but if you only look at half the case, then you will have a force.

If you applied the same thinking to all the forces and moments, you will understand the origin of each force.

For my interest, did you use cfMesh's hexahedral or polyhedral mesh?

Thomas

[Krao]

Hi Thomas,

Thank you very much for the detailed explanation. Got a better understanding of how the forces function works as well. Just for curiosity I will try the simulation for the second half whenever I get some time. It should yield same result as expected.

For my interest, did you use cfMesh's hexahedral or polyhedral mesh?

I have used Hexahedral(cartesianMesh) in cfMesh.

I appreciate your time and invaluable inputs throughout the simulation process.

Kind regards,

Krao