Autor
|
Thema: Surce Code ForceCoeffs.so und Forces.C (1004 mal gelesen)
|
Ed93 Mitglied
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 / zitieren --> Unities abgeben:
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
Beiträge: 2463 Registriert: 27.08.2010 OpenFOAM-dev (Foundation) OpenFOAM-xxxx (ESI)
|
erstellt am: 16. Okt. 2017 22:26 <-- editieren / zitieren --> Unities abgeben: Nur für Ed93
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
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 / zitieren --> Unities abgeben:
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
Beiträge: 2463 Registriert: 27.08.2010 OpenFOAM-dev (Foundation) OpenFOAM-xxxx (ESI)
|
erstellt am: 31. Okt. 2017 16:02 <-- editieren / zitieren --> Unities abgeben: Nur für Ed93
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
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 / zitieren --> Unities abgeben:
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
Beiträge: 2463 Registriert: 27.08.2010 OpenFOAM-dev (Foundation) OpenFOAM-xxxx (ESI)
|
erstellt am: 31. Okt. 2017 18:34 <-- editieren / zitieren --> Unities abgeben: Nur für Ed93
|
Ed93 Mitglied
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 / zitieren --> Unities abgeben:
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
Beiträge: 2463 Registriert: 27.08.2010 OpenFOAM-dev (Foundation) OpenFOAM-xxxx (ESI)
|
erstellt am: 31. Okt. 2017 19:13 <-- editieren / zitieren --> Unities abgeben: Nur für Ed93
|
Ed93 Mitglied
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 / zitieren --> Unities abgeben:
|
| Anzeige.:
Anzeige: (Infos zum Werbeplatz >>)
|