Hot News:

Unser Angebot:

  Foren auf CAD.de (alle Foren)
  OpenFOAM
  Surce Code ForceCoeffs.so und Forces.C

Antwort erstellen  Neues Thema erstellen
CAD.de Login | Logout | Profil | Profil bearbeiten | Registrieren | Voreinstellungen | Hilfe | Suchen

Anzeige:

Darstellung des Themas zum Ausdrucken. Bitte dann die Druckfunktion des Browsers verwenden. | Suche nach Beiträgen nächster neuer Beitrag | nächster älterer Beitrag
Autor Thema:  Surce Code ForceCoeffs.so und Forces.C (1002 mal gelesen)
Ed93
Mitglied



Sehen Sie sich das Profil von Ed93 an!   Senden Sie eine Private Message an Ed93  Schreiben Sie einen Gästebucheintrag für Ed93

Beiträge: 311
Registriert: 10.10.2015

Win 7
Intel Xeon 4x3,3GHz
24GB DDR3
Catia V5

erstellt am: 16. Okt. 2017 22:16    Editieren oder löschen Sie diesen Beitrag!  <-- editieren / zitieren -->   Antwort mit Zitat in Fett Antwort mit kursivem Zitat    Unities abgeben: 1 Unity (wenig hilfreich, aber dennoch)2 Unities3 Unities4 Unities5 Unities6 Unities7 Unities8 Unities9 Unities10 Unities

Hallo,

ich setze mich gerade mit der Widerstandsbeiwertberechnung in OpenFOAM auseinander. Ich bekomme Widerstandsbeiwerte ausgegeben und möchte diese nun interpretieren. Dazu würde ich gerne wissen, wie genau die Widerstandsbeiwerte berechnet werden und habe einen Blick auf den Source Code geworfen (forcesCoeffs.C), zudem ich ein zwei Fragen habe:

Der Widerstandsbeiwert wird berechnet mit

Code:
// Calculate coefficients
    scalar CmTot = 0;
    scalar CdTot = 0;
    scalar ClTot = 0;
    forAll(liftCoeffs, i)
    {
        momentCoeffs[i] = (moment_[i] & pitchAxis_)/(Aref_*pDyn*lRef_);
        dragCoeffs[i] = (force_[i] & dragDir_)/(Aref_*pDyn);
        liftCoeffs[i] = (force_[i] & liftDir_)/(Aref_*pDyn);

        CmTot += sum(momentCoeffs[i]);
        CdTot += sum(dragCoeffs[i]);
  ClTot += sum(liftCoeffs[i]);

Der dynamische Druck ist wie folgt definiert:  pDyn = 0.5*rhoRef_*magUInf_*magUInf_ ... alles soweit verständlich, doch Frage ich mich wie die Kraft, also force_[i] berechnet wird.

Dazu habe ich einen Blick in forces.C geworfen:

Code:
forAllConstIter(labelHashSet, patchSet_, iter)
        {
            label patchi = iter.key();

            vectorField Md
            (
                mesh_.C().boundaryField()[patchi] - coordSys_.origin()
            );

            scalarField sA(mag(Sfb[patchi]));

            // Normal force = surfaceUnitNormal*(surfaceNormal & forceDensity)
            vectorField fN
            (
                Sfb[patchi]/sA
              *(
                    Sfb[patchi] & fD.boundaryField()[patchi]
                )
            );

            // Tangential force (total force minus normal fN)
            vectorField fT(sA*fD.boundaryField()[patchi] - fN);

            // Porous force
            vectorField fP(Md.size(), Zero);

            addToFields(patchi, Md, fN, fT, fP);

            applyBins(Md, fN, fT, fP, mesh_.C().boundaryField()[patchi]);
        }
    }


Ich verstehe die Kraftberechnung leider nicht. Was ist überhaupt eine porous force? Woher stammt die forceDensity? Und wie wird diese bestimmt?

Ich hoffe, dass mir jemand weiterhelfen kann. Ich bin ein Neuling und hatte bislang mit Programmiersprachen auch wenig zu tun.


Viele Grüße

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

Shor-ty
Moderator





Sehen Sie sich das Profil von Shor-ty an!   Senden Sie eine Private Message an Shor-ty  Schreiben Sie einen Gästebucheintrag für Shor-ty

Beiträge: 2463
Registriert: 27.08.2010

OpenFOAM-dev (Foundation)
OpenFOAM-xxxx (ESI)

erstellt am: 16. Okt. 2017 22:26    Editieren oder löschen Sie diesen Beitrag!  <-- editieren / zitieren -->   Antwort mit Zitat in Fett Antwort mit kursivem Zitat    Unities abgeben: 1 Unity (wenig hilfreich, aber dennoch)2 Unities3 Unities4 Unities5 Unities6 Unities7 Unities8 Unities9 Unities10 Unities Nur für Ed93 10 Unities + Antwort hilfreich

Hallo,

kurz und prägnant. Die Porous-Force wird wahrscheinlich mit Porositäten zu tun haben. Wie du siehst ist es zu Beginn nur ein Vektorfeld mit den Werten Null. Du must dir die Funktion fP anschauen um zu verstehen was da passiert. Die forceDensity ist kommentiert und wird im Programm nicht verwendet.

Jedenfalls verwendest du nicht die OpenFOAM Foundation 5.x Version. Da gibts deinen Source-Code nämlich gar nicht.

------------------
Viele Grüße,
Tobias Holzmann

OpenFOAM Tutorials | Publikationen | Für Anfänger wiki.openfoam.com

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

Ed93
Mitglied



Sehen Sie sich das Profil von Ed93 an!   Senden Sie eine Private Message an Ed93  Schreiben Sie einen Gästebucheintrag für Ed93

Beiträge: 311
Registriert: 10.10.2015

Win 7
Intel Xeon 4x3,3GHz
24GB DDR3
Catia V5

erstellt am: 31. Okt. 2017 15:05    Editieren oder löschen Sie diesen Beitrag!  <-- editieren / zitieren -->   Antwort mit Zitat in Fett Antwort mit kursivem Zitat    Unities abgeben: 1 Unity (wenig hilfreich, aber dennoch)2 Unities3 Unities4 Unities5 Unities6 Unities7 Unities8 Unities9 Unities10 Unities

Hallo,

erstmal vielen Dank für die Antwort und entschuldige bitte die späte Antwort. Ich habe die OpenFOAM Foundatio 4.1 Version.
Mir geht es weniger um die Porosität, als darum wie die Kräfte berechnet werden.

Aber ja, ich habe gar keine forceDensity, somit war ich auf dem falschen Dampfer. Also zunächst ist das gesamte Vektorfeld 0 by default. Zur Berechnung der Kräfte ist Sfb entscheident, was den Einträgen an den Boundaries entspricht, wenn ich mich nicht täusche:

Code:
void Foam::functionObjects::forces::calcForcesMoment()
{
    initialise();

    force_[0] = Zero;
    force_[1] = Zero;
    force_[2] = Zero;

    moment_[0] = Zero;
    moment_[1] = Zero;
    moment_[2] = Zero;

    if (directForceDensity_)

        ...
else
    {
        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
        const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);

        const fvMesh& mesh = U.mesh();

        const surfaceVectorField::Boundary& Sfb =
            mesh.Sf().boundaryField();

        tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
        const volSymmTensorField::Boundary& devRhoReffb
            = tdevRhoReff().boundaryField();

        // Scale pRef by density for incompressible simulations
        scalar pRef = pRef_/rho(p);

        forAllConstIter(labelHashSet, patchSet_, iter)
        {
            label patchi = iter.key();

            vectorField Md
            (
                mesh.C().boundaryField()[patchi] - coordSys_.origin()
            );

            vectorField fN
            (
                rho(p)*Sfb[patchi]*(p.boundaryField()[patchi] - pRef)
            );

            vectorField fT(Sfb[patchi] & devRhoReffb[patchi]);

            vectorField fP(Md.size(), Zero);

            applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchi]);
        }
    }




Hier ist ersichtlich, wie fN bestimmt wird. Ich muss dazu aber wissen, wie Sfb definiert wird:

Code:
Foam::tmp<Foam::volSymmTensorField>
Foam::functionObjects::forces::DevRhoReff() const
{
    typedef compressible::turbulenceModel cmpTurbModel;
    typedef incompressible::turbulenceModel icoTurbModel;

    if (obr_.foundObject<cmpTurbModel>(cmpTurbModel: ropertiesName))
    {
        const cmpTurbModel& turb =
            obr_.lookupObject<cmpTurbModel>(cmpTurbModel: ropertiesName);

        return turb.devRhoReff();
    }
    else if (obr_.foundObject<icoTurbModel>(icoTurbModel: ropertiesName))
    {
        const incompressible::turbulenceModel& turb =
            obr_.lookupObject<icoTurbModel>(icoTurbModel: ropertiesName);

        return rho()*turb.devReff();
    }
    else if (obr_.foundObject<fluidThermo>(fluidThermo::DictName))
    {
        const fluidThermo& thermo =
            obr_.lookupObject<fluidThermo>(fluidThermo::DictName);

        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);

        return -thermo.mu()*dev(twoSymm(fvc::grad(U)));
    }
    else if
    (
        obr_.foundObject<transportModel>("transportProperties")
    )
    {
        const transportModel& laminarT =
            obr_.lookupObject<transportModel>("transportProperties");

        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);

        return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
    }
    else if (obr_.foundObject<dictionary>("transportProperties"))
    {
        const dictionary& transportProperties =
            obr_.lookupObject<dictionary>("transportProperties");

        dimensionedScalar nu(transportProperties.lookup("nu"));

        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);

        return -rho()*nu*dev(twoSymm(fvc::grad(U)));
    }
    else
    {
        FatalErrorInFunction
            << "No valid model for viscous stress calculation"
            << exit(FatalError);

        return volSymmTensorField::null();
    }
}


Der für mich zutreffende Abschnitt ist meiner Meinung nach:

Code:
else if (obr_.foundObject<dictionary>("transportProperties"))
    {
        const dictionary& transportProperties =
            obr_.lookupObject<dictionary>("transportProperties");

        dimensionedScalar nu(transportProperties.lookup("nu"));

        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);

        return -rho()*nu*dev(twoSymm(fvc::grad(U)));
    }


Die letzte Zeile gibt mir somit an, wie mein Vektorfeld berechnet wird.
Ich bitte um Verbesserung, wenn ich falsch liege! Ich Frage mich noch was twoSymm und fvc:: bedeutet.

Außerdem habe ich evtl einen Denkfehler, ich komme bei der Berechnung der Kraft: rho(p)*Sfb[patchi]*(p.boundaryField()[patchi] - pRef) nicht auf die Einheit Newton, vielmehr auf einen Ausdruck mit kg³, da in rho,p und Sfb jeweils kg im Nenner stehen. Wie kann das sein? Wo liegt mein Fehler?

Viele Grüße und einen schönen Feiertag!

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

Shor-ty
Moderator





Sehen Sie sich das Profil von Shor-ty an!   Senden Sie eine Private Message an Shor-ty  Schreiben Sie einen Gästebucheintrag für Shor-ty

Beiträge: 2463
Registriert: 27.08.2010

OpenFOAM-dev (Foundation)
OpenFOAM-xxxx (ESI)

erstellt am: 31. Okt. 2017 16:02    Editieren oder löschen Sie diesen Beitrag!  <-- editieren / zitieren -->   Antwort mit Zitat in Fett Antwort mit kursivem Zitat    Unities abgeben: 1 Unity (wenig hilfreich, aber dennoch)2 Unities3 Unities4 Unities5 Unities6 Unities7 Unities8 Unities9 Unities10 Unities Nur für Ed93 10 Unities + Antwort hilfreich

Ich habs nur überflogen:

fvc:= finite volume calculus -> explizit mit alten Werten. Implizit geht hier auch gar nicht
twoSymm() ist der symmetrische Anteil der Matrix von grad(U) multipliziert mit 2. Das kommt vom Spannungstensor Sigma und dem Deformationstensor.

Wie kommst du den eigentlich drauf das es genau die letzte if else Anweisung ist. Kommt halt drauf an was du für Settings hast. Entsprechend wird eben der Druckteil vom Spannungstensor zurückgegeben.

------------------
Viele Grüße,
Tobias Holzmann

OpenFOAM Tutorials | Publikationen | Für Anfänger wiki.openfoam.com

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

Ed93
Mitglied



Sehen Sie sich das Profil von Ed93 an!   Senden Sie eine Private Message an Ed93  Schreiben Sie einen Gästebucheintrag für Ed93

Beiträge: 311
Registriert: 10.10.2015

Win 7
Intel Xeon 4x3,3GHz
24GB DDR3
Catia V5

erstellt am: 31. Okt. 2017 18:21    Editieren oder löschen Sie diesen Beitrag!  <-- editieren / zitieren -->   Antwort mit Zitat in Fett Antwort mit kursivem Zitat    Unities abgeben: 1 Unity (wenig hilfreich, aber dennoch)2 Unities3 Unities4 Unities5 Unities6 Unities7 Unities8 Unities9 Unities10 Unities

Die letze if else, weil ich eine dictionary namens "TransportProperties" habe, in dem ich mein RAS Modell defniniert habe und zudem inkompressibel rechne. Deswegen ist es für mich die naheliegendste if else Bedingung.

Super danke, das hilft mir schon einmal weiter! Ich frage mich noch, was dev() bedeutet, die Funktion wird ebenfalls auf grad(U), bzw. den symmetrischen Anteil davon angewandt. Ich dachte zunächst an die Divergenz, aber es heißt ja nicht div().

Dev taucht auch in der tangential Kraft wieder auf, aber nicht als Funktion, sondem im devRhoReffb:

vectorField fT(Sfb[patchi] & devRhoReffb[patchi]);

Ich weiß auch noch nicht genau, was mir dieser Ausdruck sagt. Hier müsste ja eigentlich die Reibung für die Tangentialkraft berücksichtig werden. In fT sollten alle Einträge stehen, bei denen Sfb und devRHoReffb je patch übereinstimmen, jedoch weiß ich nicht genau, was devRhoReffb eigentlich ist.


Viele Grüße,
ED93

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

Shor-ty
Moderator





Sehen Sie sich das Profil von Shor-ty an!   Senden Sie eine Private Message an Shor-ty  Schreiben Sie einen Gästebucheintrag für Shor-ty

Beiträge: 2463
Registriert: 27.08.2010

OpenFOAM-dev (Foundation)
OpenFOAM-xxxx (ESI)

erstellt am: 31. Okt. 2017 18:34    Editieren oder löschen Sie diesen Beitrag!  <-- editieren / zitieren -->   Antwort mit Zitat in Fett Antwort mit kursivem Zitat    Unities abgeben: 1 Unity (wenig hilfreich, aber dennoch)2 Unities3 Unities4 Unities5 Unities6 Unities7 Unities8 Unities9 Unities10 Unities Nur für Ed93 10 Unities + Antwort hilfreich

Wenn du ein RANS Modell hast dann gehst du wohl in die zweite ifelse. Du solltest schon wissen was du machst. dev ist der Deviatorische Anteil der Matrix die du Berechnest. Heißt also der hydrostatische Anteil (1/3 tr(Matrix) I) wird abgezogen.

------------------
Viele Grüße,
Tobias Holzmann

OpenFOAM Tutorials | Publikationen | Für Anfänger wiki.openfoam.com

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

Ed93
Mitglied



Sehen Sie sich das Profil von Ed93 an!   Senden Sie eine Private Message an Ed93  Schreiben Sie einen Gästebucheintrag für Ed93

Beiträge: 311
Registriert: 10.10.2015

Win 7
Intel Xeon 4x3,3GHz
24GB DDR3
Catia V5

erstellt am: 31. Okt. 2017 19:01    Editieren oder löschen Sie diesen Beitrag!  <-- editieren / zitieren -->   Antwort mit Zitat in Fett Antwort mit kursivem Zitat    Unities abgeben: 1 Unity (wenig hilfreich, aber dennoch)2 Unities3 Unities4 Unities5 Unities6 Unities7 Unities8 Unities9 Unities10 Unities

Ah super, danke!

Ist meine erste Rechnung und ich versuche nachzuvollziehen wie das Programm arbeitet. Habe da wenig Erfahrung und tue mich etwas schwer, sry, aber irgendwo muss ich ja anfangen.

Wieso gehe ich da in die das zweite if else?

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

Shor-ty
Moderator





Sehen Sie sich das Profil von Shor-ty an!   Senden Sie eine Private Message an Shor-ty  Schreiben Sie einen Gästebucheintrag für Shor-ty

Beiträge: 2463
Registriert: 27.08.2010

OpenFOAM-dev (Foundation)
OpenFOAM-xxxx (ESI)

erstellt am: 31. Okt. 2017 19:13    Editieren oder löschen Sie diesen Beitrag!  <-- editieren / zitieren -->   Antwort mit Zitat in Fett Antwort mit kursivem Zitat    Unities abgeben: 1 Unity (wenig hilfreich, aber dennoch)2 Unities3 Unities4 Unities5 Unities6 Unities7 Unities8 Unities9 Unities10 Unities Nur für Ed93 10 Unities + Antwort hilfreich

Wenn du turbulent rechnest, gibt es ein turbulence Objekt in der FOAM db() ... die Berechnung muss ja anders sein als laminar aufgrund des Spannungstensors Sigma. Ich verweise diesbezüglich mal auf mein Buch.

------------------
Viele Grüße,
Tobias Holzmann

OpenFOAM Tutorials | Publikationen | Für Anfänger wiki.openfoam.com

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

Ed93
Mitglied



Sehen Sie sich das Profil von Ed93 an!   Senden Sie eine Private Message an Ed93  Schreiben Sie einen Gästebucheintrag für Ed93

Beiträge: 311
Registriert: 10.10.2015

Win 7
Intel Xeon 4x3,3GHz
24GB DDR3
Catia V5

erstellt am: 01. Nov. 2017 19:23    Editieren oder löschen Sie diesen Beitrag!  <-- editieren / zitieren -->   Antwort mit Zitat in Fett Antwort mit kursivem Zitat    Unities abgeben: 1 Unity (wenig hilfreich, aber dennoch)2 Unities3 Unities4 Unities5 Unities6 Unities7 Unities8 Unities9 Unities10 Unities

Ok, danke für die Hilfe! Etwas schlauer bin ich, hab aber noch einen weiten Weg vormir.

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

Anzeige.:

Anzeige: (Infos zum Werbeplatz >>)

Darstellung des Themas zum Ausdrucken. Bitte dann die Druckfunktion des Browsers verwenden. | Suche nach Beiträgen

nächster neuerer Beitrag | nächster älterer Beitrag
Antwort erstellen


Diesen Beitrag mit Lesezeichen versehen ... | Nach anderen Beiträgen suchen | CAD.de-Newsletter

Administrative Optionen: Beitrag schliessen | Archivieren/Bewegen | Beitrag melden!

Fragen und Anregungen: Kritik-Forum | Neues aus der Community: Community-Forum

(c)2023 CAD.de | Impressum | Datenschutz