[hoemmaria]

Hello!

I am using OpenFoam 4.1 and I am trying to implement an actuator disk. The problem I am having is that velocity and pressure is not calculated. At the first time step it says "No iterations = 0". The initial velocity and pressure files 0/U and 0/p are copied to all the following time steps. This is obviously not what I want. Here are my files and process:

I want to use the simpleFoam solver and I introduce the actuation disk in a fvOptions file as follows:

constant/fvOptions:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

disk1
{
    type            actuationDiskSource;
    active          on;

    actuationDiskSourceCoeffs
    {
        fields      (U);

        selectionMode   cellSet;
        cellSet         actuationDisk1;
        diskDir         (1 0 0);    // Orientation of the disk
        Cp              0.4;
        Ct              0.5;
        diskArea        0.16;
        upstreamPoint   (4.9 0.5 0.5);
    }
}

// ************************************************************************* //

0/p :

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
      type          zeroGradient;

    }
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
}

0/U:

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            pressureInletVelocity;
        value           uniform (0 0 0);
    }
    outlet
    {
        type            inletOutlet;
        inletValue      uniform (0 0 0);
        value           uniform (0 0 0);
    }
    walls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

I am using the simpleFoam solver and I have not made any changes to the pEqn.H, UEqn.H, simpleFoam.C or createFields.H.

My blockMeshDict only includes inlet, outlet and walls.

I use topoSet with the following file:

system/topoSetDict

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      topoSetDict;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

actions
(
    // actuationDisk1
    {
        name    actuationDisk1CellSet;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (4.9 0.1 0.1) (5.1 0.9 0.9);
        }
    }
    {
        name    actuationDisk1;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo
        {
            set actuationDisk1CellSet;
        }
    }
);

Both epsilon and k are computed as normal when run the code, but nothing happens to neither the velocity nor the pressure. I have also tried the turbineSiting tutorial where the actuationDiskSource is being used, and it runs without problem. But I do not seem to find the cause of my case failing compared to the turbineSiting tutorial.

I any of you have tried something similar before or have any suggestions to what I might be doing wrong, I am most grateful for help.


Kind regards,
Maria Hoem

[lebc]

Hi Maria,

I use OF3.0.1, so maybe there are a couple of changes that I'm not aware, but maybe I can help you (I'm using actuationDiskSource intensively... )

In 0/U file, try to change the inlet by this (attention to the "Uinlet" line):

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

Uinlet (10 0 0);

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

internalField uniform $flowVelocity;

boundaryField
{

inlet
{
type fixedValue;
value uniform $Uinlet;
}

About the outlet, I kept the same as the turbine siting tutorial.

I believe the 0/p file is Ok, just try to modify the U file and tell me if it is working now!

By your U file it seems that you are using U=0 at the inlet, so there is no fluid flowing!

To run my cases I have modified the turbine siting tutorial (one step at each time), maybe it is easier if you try it!

Just by curiosity, what are you using actuationDiskSource for?!

Best Regards,

Luis

[hoemmaria]

Hi, Luis.

Thank you for your answer!

I have now tried to change the inlet velocity to be non-zero like you suggested, and modified the turbineSiting tutorial step by step. I solved the new case using both SimpleFoam and my own solver actuationDiskSourceSimpleFoam where actuationDiskSource.C/.H and actuationDiskSourceTemplates.C is included. They both give me the same solution, but it does not seem like the actuator disk is doing anything.
There is no pressure jump over the actuator disk like I want there to be. I have attached a photo of the velocity and pressure through the tunnel. Any further suggestions?

I am using actuationDiskSource for simulating a simplified wind turbine in a wind tunnel for my master thesis. And I am going to compare the results with measured data from the tunnel.

Kind regards,
Maria

[lebc]

Hi Maria,

Could you send me your case files?!

You said that you have created a new solver, but it is not needed to use the actuationDiskSource (unless you did any other changes...)

Nice, I'm also using it for my master thesis, I want to simulate a wind farm using these kind of models! Maybe we can exchange some knowledge! Take my email, we can compare our results and help each other (I say that because I'm having a hard time using OF for this!): luis.bonicruz@gmail.com

If you are using the atmBoundarylayer (not sure if this is the name...) you have also to change the wind speed there.

Best Regards,
Luis

[louisgag]

Hello Maria and Luis,
I can confirm that for me as well on OF-4.x there is no pressure jump when using the actuationDiskSource method. I think the rotorDisk does however work. I am currently investigating...
Regards,
-Louis

[hoemmaria]

Thank you for your response, Luis There are way too many files to upload here in the forum, so I have made a google drive link where you, and anyone else who wishes to, can view my case files: https://drive.google.com/open?id=0B0...2JJd0lqMEhhWms

I guess I will just stick to SimpleFoam then

And I tried changing the ABLConditions by changing the velocity and zref to be within my mesh. It did not do any large changes (I think), only lowering the pressure in total, but still no jump.

Hello Maria and Luis,
I can confirm that for me as well on OF-4.x there is no pressure jump when using the actuationDiskSource method. I think the rotorDisk does however work. I am currently investigating...
Regards,
-Louis

Thank you Louis for your input. Looking forward to your investigations

[louisgag]

Hello again,

Well my problem was simply that I was trying to simulate a rotor using the actuationDiskSource method, I understood this after having tried different settings and looked at the code. So, in the meantime I was also able to successfully run a test case of a wind turbine with the actuationDiskSource, getting both pressure and velocity jumps across the disk using pimpleFoam.

I looked at your case Maria, although I didn't run it because all the files are in the same directory. First thing I note is, as Luis pointed out, you have a null velocity. This means the wind turbine won't do anything... Also, make sure diskDir points downstream behind the turbine, which may already be your case...

I would also suggest using a different inletValue velocity at the outlet than the one used at the actual inlet. The details on this boundary conditions are given here: https://cfd.direct/openfoam/user-guide/boundaries/

Hope this helps, otherwise please upload the case with the actual file/folder structure ;-)

Best Regards,


-Louis

[hoemmaria]

I am now trying with a non-zero inlet velocity and a InletValue velocity at the outlet twice as big as the actual inlet velocity. I am running it with pimpleFoam, however it takes quite a long time on my computer. Will know if it worked later. Or do you mean that I should have a non-zero internalField velocity?

Otherwise I have uploaded the case files in folders - of course I should have done this right away. They are available at: https://drive.google.com/open?id=0B0...2JJd0lqMEhhWms

Thank you for taking the time to look into the problem!

- Maria

[lebc]

Hi Maria,

I've checked your files, and they seem to be all ok... Basically your actuationDisk is not "being seen" by the solver. Unfortunately I didn't find an answer to it, but I'll keep trying!

@louis

The must be a pressure jump, otherwise it would not have an actuator disk behaviour...

You should expect something like the picture below.

My advices for now:
- Keep as close as possible from the tutorial, and check your U file on 0 folder, it is a little bit different from the tutorial.
- Try a bigger domain, just to avoid the wall boundary layer (for now) and make sure your model is doing what it was supposed to do. After you can start improving it.
- For simplicity, use simpleFoam (it should work, I have only used it here!)

I wish I could work on it more time, but during the day I can't analyze it at work...

Try these and let me know if it worked!

Best Regards,
Luis

[louisgag]

Hello Maria,

I've got your case running fine by making two modifications:

1: as Luis suggested, remove the boundary layer from the walls. I did this by setting the walls to slip:

    walls
    {
        type            slip;
    }

2: 1 m/s is small, should work, but try putting something strong, like 20 m/s, if you want to see marked pressure and velocity jumps.

Uinlet (20 0 0);

For internalField velocity, you can make it equal to the inlet but it won't change much as it is only the initial condition. This case should run fast, not sure why it would be slow on your computer... seems to converge before your convergence criteria are met, you might want to make them bigger to have the solve stop way before 5000 iterations...

@Luis, I am modeling a powered rotor, such as those for helicopters, so this class is not for me ;-) thanks for you input.

Regards,


-Louis

[hoemmaria]

Thank you to both of you

I have now changed all the initial files to have wall boundary condition as slip and changed Uinlet to (20 0 0) as you have suggested. I do indeed get a promising pressure jump, however I am not sure about the velocity, it barely changes along the symmetry line (see attached figure). The actuation disk seem to be acting equally in both x-directions (positive and negative) from the midpoint at x=4.5, and there is no wake as shown in your illustrative figure, Luis.

Did you get the same when you ran it, Louis?

[louisgag]

Yes, I rechecked and the same velocities occur on my version of OF. Perhaps the larger domain that Luis suggested would help. As far as I am concerned I can't help you much more as I have no experience with wind turbines.
Good luck,
-Louis

[lebc]

Hi Maria and Louis,

@Maria

Could you send your new files? About the files size, you don't need to send the mesh files generated by blockMesh, just the blockMeshDict , and this way the files will be very small.

Additionally, you can try running parallel, its very useful! My mesh is not so fine, and the convergence is after about 1000 iterations...

The results I got with your initial files were the same as that if no wind turbine was present.

@Louis

I'm using the slip conditions for all boundaries, except for inlet and outlet, and a big domain. My internalField is set o (0 0 0) and inlet velocity is 10m/s.

I was also testing the rotorDiskSource, in order to check if it is suitable to represent a wind turbine, but as it is a little bit more complicated I left the investigation for further steps...

Best Regards,
Luis

[hoemmaria]

Thank you Louis for your effort. I am grateful for your help.

@Luis

Yeah, if the system do not see any actuator disk, that is a problem. I have added my files in a -zip file.
Run blockMesh, topoSet and then simpleFoam.

[lebc]

Hi Maria,

I didn't have much time to work on your problem, so I'm sharing one case that I have successfully run some months ago, and you can use it as a base to adjust your case.

Some advises:

1- Keep working on your case, use mine just as reference, don't turn it into yours, because this way you can acquire some knowledge;

2- Try to modify file by file, and if it works you will know where the mistake was;

3- In case you have any cp or ct curve for small turbines as yours (~1m diameter), try to compare them with the values you are using. Although I work with wind turbines, I only have knowledge on the MW scale (~100m diameter), so it is possible to have a different behavior for cp and ct values for your scale.

I will keep analyzing your case files in my free time, please share your results if you have any progress!

Best Regards,
Luis

[hoemmaria]

@Luis

I think I solved it with the help of your files

Things I changed and observed:

1. I had too few cells for the actuator disk to be seen, it was smaller than the cell i x-direction. Thus, adding snappyHexMesh in that area and make the mesh finer at the disk area helped. Or I could have made the whole mesh finer.
2. I made the actuator disk region definition in topoSetDict a cylinder (like in you example) instead of a box.
3. I made the blockMeshDict to consist of only one box, not two (it was two boxes in the example I had gotten the file from). Not sure if this did any difference on its own, but its more tidy at least.
4. I have added the postprocessing option residuals to controlDict to be able to see if the solution converges. The file that follows is placed in "system". To view the residuals, type: foamMonitor -l postprocessing/residuals/0/residuals.dat.

I have also added my final files if anyone wants to look at them or have similar problems in the future.

Thank you so much for all the help!

Kind regards,
Maria Hoem

[lebc]

Hi Maria,

Very nice!!!

In addition, you can adapt the "Allrun" file I sent within my case, so you can run your simulations faster, and avoid having the information on terminal (instead, everything will be written to log files, which you can open anytime during the simulation).

Also, try to use several refinement regions, you probably will get more accurate results on a more advanced stage of your research.

Could you post a picture of a slice from your results? I have run your case here, but as I could run just some time steps, maybe your results are more representative of what is happening then mine!

Best Regards,
Luis

[kokab]

Hi dears
I also work with AD, my geometry is channel at x as streamwise, z as spanwise and y is the height direction of my channel. I want to implement oblique actuator disk with angle of 30 to wind direction ( which is in streamwise ). I have made some changes in diskDir but there is not any differences in a result. can anyone help me with that?
---------------------------------------------------------------------------------------
my fvOption:
disk1
{
type actuationDiskSource;
active on;
selectionMode cellSet;
cellSet actuationDisk1;

actuationDiskSourceCoeffs
{
fieldNames (U);
diskDir (-0.866025403784 0 -0.5);
Cp 0.53;
Ct 0.58;
diskArea 1;
upstreamPoint (3.2 1 1.5);
}
-------------------------------------------------------------------------------------------------
and topoSetFieldDict:

name actuationDisk1;
type cellSet;
action new;
source boxToCell;
sourceInfo
{
box (3 0.5 1) (3.1 1.5 2);
}
}
{
name actuationDisk1;
type cellZoneSet;
action new;
source setToCellZone;
sourceInfo
{
set actuationDisk1;
}

[lebc]

Hi Kokab,

Have you tried to see your mesh on Paraview? It is possible to see the disk region set by topoSet, this way you can try to figure something out.

Also, when you use "source boxToCell", the box faces will be aligned with the (x y z) axis, so I believe you will not see a box rotate 30° (I didn't try it! ). As I don't know which version of OF you are using, I will assume that it is the lastest, so you can take a look at the following link and check the source options on topoSetDict.

https://github.com/OpenFOAM/OpenFOAM...et/topoSetDict

Best Regards,
Luis

[kokab]

Hi dear lebec
thanks a lot for your answer, how can I see the disk region? ( I know it is embarrassing question ) I will try on toposet and tell you the result soon, maybe I need your help more. My OF version is 2.4.0 ( I did work with 3.0.1 too )