CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

verification of a new BC whit codedfixedValue

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 4, 2020, 05:46
Default help for new BC with codedfixedValue
  #1
Neb
Member
 
Join Date: Mar 2020
Posts: 66
Rep Power: 6
Neb is on a distinguished road
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!!
Attached Images
File Type: png equazione.png (7.1 KB, 30 views)

Last edited by Neb; June 4, 2020 at 11:56.
Neb is offline   Reply With Quote

Old   June 4, 2020, 13:25
Default Some comments
  #2
Senior Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9
crubio.abujas is on a distinguished road
Hi Neb,

Some things to consider:
  • I think there is a small typo in the formula. Just a minus sign in the pressure corrector.
Quote:
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:
Code:
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.
crubio.abujas is offline   Reply With Quote

Old   June 5, 2020, 03:57
Default
  #3
Neb
Member
 
Join Date: Mar 2020
Posts: 66
Rep Power: 6
Neb is on a distinguished road
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?
Neb is offline   Reply With Quote

Old   June 5, 2020, 06:01
Default A solution for uniform application
  #4
Senior Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9
crubio.abujas is on a distinguished road
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 e^{x+y}=e^xe^y you can split the weigth term inside the integral into:e^{\frac{\tau-t}{R_2C}} = e^{\frac{\tau}{R_2C}}e^{\frac{-t}{R_2C}}. Now the second part is not dependent on the integrand variable \tau 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:
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);
    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 10^-9 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 e^{\frac{\tau}{R_2C}} will be inmense.
crubio.abujas is offline   Reply With Quote

Old   June 5, 2020, 07:42
Default
  #5
Neb
Member
 
Join Date: Mar 2020
Posts: 66
Rep Power: 6
Neb is on a distinguished road
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'
Neb is offline   Reply With Quote

Old   June 5, 2020, 08:08
Thumbs up
  #6
Senior Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9
crubio.abujas is on a distinguished road
Quote:
Originally Posted by Neb View Post
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==( ... )

Code:
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.
crubio.abujas is offline   Reply With Quote

Old   June 5, 2020, 09:19
Default
  #7
Neb
Member
 
Join Date: Mar 2020
Posts: 66
Rep Power: 6
Neb is on a distinguished road
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?
Neb is offline   Reply With Quote

Old   June 5, 2020, 10:58
Default
  #8
Senior Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9
crubio.abujas is on a distinguished road
Glad to know it worked!


Quote:
Originally Posted by Neb View Post
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!
crubio.abujas is offline   Reply With Quote

Old   June 5, 2020, 11:02
Default
  #9
Neb
Member
 
Join Date: Mar 2020
Posts: 66
Rep Power: 6
Neb is on a distinguished road
Thanks! It will certainly be fun!

Have a nice day
Neb is offline   Reply With Quote

Reply

Tags
codedfixedvalue, integral computing, new bc

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
including parameter file in codedFixedValue Loekatoni OpenFOAM Running, Solving & CFD 4 November 9, 2023 16:56
Why does the result not the same as the verification manuals? Heqian Zhao FLUENT 1 April 20, 2020 15:29
codedFixedValue: accessing other patch causes crash in parallel RL-S OpenFOAM Running, Solving & CFD 2 December 24, 2019 21:20
Problem to run simpleFoam using qsub? be_inspired OpenFOAM 1 December 22, 2015 12:53
Using codedFixedvalue to apply totalPressure Boundary Condition cdm OpenFOAM Running, Solving & CFD 2 June 22, 2013 14:10


All times are GMT -4. The time now is 15:06.