[Neb]

Hi at all,
I created a newBC on the outlet for the pressure using codedfixedvalue.
How can i import the flow rate at time 0 of the outlet as the initial value in my new code (Q_0)?



code
#{

const scalar &t = this ->db().time().value(); // time reading

const surfaceScalarField& phi = db().lookupObject<surfaceScalarField>("phi"); //flowrate reading in the outlet
const fvsPatchField<scalar>& phip=patch().patchField<surfaceScalarField,scalar> (phi);

const scalar &dt = this ->db().time().deltaTValue(); //delta t reading for integration

const scalar C = 4.82e-9; //m4*s2/kg
const scalar Rp= 2.79245; // kg*s/m4
const scalar Rd= 1.7; //kg/s*m4
const scalar p_0=1.07; //pressure t=0 m2/s2
const scalar Q_0=0; //flow rate t=0 m3/s
const scalar e(M_E); //nepero number

operator==(Rp*gSum(phip)+(p_0+Rp*Q_0)*(pow(e,(-t/(Rd*C))))+((1/C)*(gSum(pow(e,(t-(t/(Rd*C))))*(phip)))*dt));
#};

thanks!!

[crubio.abujas]

Hi Neb,

Some things to consider:

• I think there is a small typo in the formula. Just a minus sign in the pressure corrector.

operator==(Rp*gSum(phip)+(p_0-Rp*Q_0)*(pow(e,(-t/(Rd*C))))+((1/C)*(gSum(pow(e,(t-(t/(Rd*C))))*(phip)))*dt));

• There are some redundant parenthesis in your formula. Probably not doing any harm, but still I think readability matters.

• OF have tools to manage the dimensions, so maybe is a good idea to use them:

const dimensionedScalar C   (dimensionSet(-1, 4, 2, 0, 0, 0, 0), 4.82e-9);
const dimensionedScalar Rp  (dimensionSet(1, -4, -1, 0, 0, 0, 0), 2.79245);
const dimensionedScalar Rd  (dimensionSet(1, -4, -1, 0, 0, 0, 0), 1.7);

• As far as I know gSum will return you a scalar with the summation of the field (the patch in this case), but you need to impose a scalarField, so there is a value for each face of the path. So I think p_0, and Q_0 shall be scalarFields, and you shall remove the gSum and handle all the patch values.

• There is a problem with your approach with the integral. You need a summation along time, not just in an instant. I first thought about having a static variable initialized to 0 in the beginning and adding up the integrand on the subsequent steps. The problem is that the weight factor of the volFlow (the exponential inside the integral) depend on the integral limits, and thus the relative importance of a given instant changes along the simulation. So I think you will need to carry all the previous values of Q to solve that integral. That or find any ingenious analytical way to handle it.

And there are my two cents on the problem. I hope that helps you handling the problem.

[Neb]

Thanks for your answer, you have been very helpful!!

Unfortunately I recently started with OpenFoam and I am in difficulty. I don't even know the language used..
I figured there were these problems on the integral and the flow rate calculation in the outlet.

• As far as I know gSum will return you a scalar with the summation of the field (the patch in this case), but you need to impose a scalarField, so there is a value for each face of the path. So I think p_0, and Q_0 shall be scalarFields, and you shall remove the gSum and handle all the patch values.

I can't understand what I have to do here. If I leave gSum (phip) does it not calculate the sum of the flow rate on my patch (outlet) at each instant of calculation? I am sorry for this question that may seem simple, but I searched a lot and found nothing that could clarify my ideas. Do I have to do a cycle for every t?

[crubio.abujas]

Ok, I think I know what are you trying. I was thinking that you want that each of the faces conforming the patch has a different value, but maybe want you wanted is to have a uniform pressure on all the patch, right?


In that case you almost got it in the first place, I think you just needed to reconsider the integral part. I was thinking about how to deal with this term and I think it is in fact quite easy to handle. Using the exponential property you can split the weigth term inside the integral into:. Now the second part is not dependent on the integrand variable so you can put outside the integral. With this approach, what is it inside the integral is not dependent on the integral limits so you can track the sum in an external variable and update.

The resultant code is:

#{
    const scalar t = this->db().time().value(); // time reading

    const surfaceScalarField& phi = db().lookupObject<surfaceScalarField>("phi");   //flowrate reading in the outlet
    const fvsPatchField<scalar>& phip=patch().patchField<surfaceScalarField,scalar> (phi);

    const scalar dt = this ->db().time().deltaTValue(); //delta t reading for integration


    const scalar C   (4.82e-9);
    const scalar Rp  (2.79245);
    const scalar Rd  (1.7);
    const scalar p_0 (1.07);
    const scalar Q_0 (0);
    const scalar e(M_E); //nepero number

    static scalar integrand(0);
    integrand += pow(e,t/(Rd*C))*gSum(phip)*dt;

    operator==(
        Rp*gSum(phip) + pow(e,-t/(Rd*C))* ((p_0-Rp*Q_0) + (1/C)*integrand)
    );
#};

I've tried an it compiles without error in OF7, so let's hope it works for you.



PD: You may want to double check the constants. The actual value of C is of the order and Rp is of order 1, so the factor multiply the exponential is really small. I don't know where do you want to apply this formula, but to match these values you're need operation times on the order of 100 ns. Once you do times higher than those there is no point in having to take into account the integrand part, in any practical effect the value if the integral is just the value of the last instant. Moreover, the approach of split the exponential term will blow the simulation will be inmense.

[Neb]

Thanks!! This is of great help to me!
Yes, that is whant I want to do.

I have to calibrate the constants better, so they are not correct. C is probably too small.

I entered the code in my calculation but returns the following error:
"invalid operands of types 'void' and 'Foam::scalar{aka double}' to binary 'operator'

[crubio.abujas]

Originally Posted by Neb

Thanks!! This is of great help to me!
Yes, that is whant I want to do.

I have to calibrate the constants better, so they are not correct. C is probably too small.

I entered the code in my calculation but returns the following error:
"invalid operands of types 'void' and 'Foam::scalar{aka double}' to binary 'operator'

It seems that there is some king of error in the operator. I don't know which function can be returning void. What version of OF are you using? Maybe I can check the code guide.

You can try some alternative way to set the field. For example using this instead of the operator==( ... )

vectorField& field = *this;
forAll(field, facei) {
    field[facei] = (Rp*gSum(phip) + pow(e,-t/(Rd*C))* ((p_0-Rp*Q_0) + (1/C)*integrand)) 
 }

If the problem persists you can try test a dummy value i.e operator==(Rp); and check if that works. On success try adding a the next term gSum(phip) and so on. That way you can locate which point is giving you troubles.

[Neb]

I found my mistake. Simply a parenthesis problem after operator.
Now everything works! Thank you very much for your help, you have been very kind!

One last question: could you recommend any material to study to learn this programming?

[crubio.abujas]

Glad to know it worked!

Originally Posted by Neb

One last question: could you recommend any material to study to learn this programming?

You can check the following link to some tutorials to learn how to program in OF.

OpenFOAM programming tutorials for beginners


You can download it and play around. The tutorials are well documented and you can play around with them.



It may be convenient to have a good comprehension of C++, and learn how to consult doxygent to understand how the different objects are made. You can have some grasp about C++ in this link:
http://www.cplusplus.com/doc/tutorial/


You may not know what is doxygen. It is a web-based reference of all the classes and funcitions that constitute OpenFoam. You can create a local version of your own but you can also check the online version:
https://cpp.openfoam.org/v7/


It may be a little be intimidating at first, with all the templates, pointers, clasess and abstract concepts, but do not despair. I would start trying to apply some solution with the tutorial and in parallel checking some manual of C++. Checking from time to time something in doxygen to get use to it.



Good luck!

[Neb]

Thanks! It will certainly be fun!

Have a nice day