diff --git a/PotentialFlow/generateOptFile.py b/PotentialFlow/generateOptFile.py deleted file mode 100644 index 7f121d051cbe7be23549e32571fc38bc8a57d24a..0000000000000000000000000000000000000000 --- a/PotentialFlow/generateOptFile.py +++ /dev/null @@ -1,26 +0,0 @@ -#!/usr/bin/env python -# This Python script generates a formay file for the "phi.pos" -# map generated by the model. - -#print input[0] -phiTrail = input[0] - -fout = open('KJ.pos.opt', 'w') -fout.write('l=PostProcessing.NbViews-1;\n') -fout.write('View[l].IntervalsType = 3;\n') -fout.write('Mesh.SurfaceEdges = 0;\n') -fout.write('View[l].RangeType = 2;\n') -fout.write('View[l].CustomMin = ' + str(phiTrail) + ';\n') -fout.write('View[l].CustomMax = ' + str(1.001*phiTrail) + ';\n') -fout.write('View[l].NbIso = 10;\n\n') - -fout.close() - - - - - - - - - diff --git a/PotentialFlow/magnus.pro b/PotentialFlow/magnus.pro index ab45cee046c9fedff1f4b85122b3638876a14209..c2597cf817c48419aa7eb5bd7de2cf9dc58786a9 100644 --- a/PotentialFlow/magnus.pro +++ b/PotentialFlow/magnus.pro @@ -2,14 +2,13 @@ Tutorial 6 : Potential flow, lift and Magnus effect Features: - - Discontinuous scalar field + - Potential flow, irrotational flow, lift and Magnus effect - Multivalued scalar field - - Potential flow - Non-linear iteration on an Associated global quantity - - Using a run-time variable in a Post-Operation (via Python) + - Use a run-time variable in a Post-Operation To compute the solution in a terminal: - getdp magnus.pro -solve LPE -pos LPE + getdp magnus.pro -solve PotentialFlow -pos PotentialFlow gmsh magnus.geo velocity.pos To compute the solution interactively from the Gmsh GUI: @@ -77,19 +76,21 @@ This condition states that the actual value of "Dmdt" is the one for which the stagnation point (singularity of the velocity field) - is located at the trailing edge of the airfoil, - i.e., when the function + is located at the trailing edge of the airfoil,. + This is true when the function argVTrail[] = ( Atan2[CompY[$1], CompX[$1] ] - Incidence ) / deg ; returns zero when evaluated for the velocity at a point "P_edge" - close to the trailing edge. + close to the trailing edge of the airfoil. A non linear iteration is thus done by means of a pseudo-newton scheme, updating the value of the imposed "Dmdt" until the norm of the residual f(V) = argVTrail[ {d phi } ] OnPoint {1.0001,0,0} - comes below a fixed tolerance. + comes below a fixed tolerance. + + */ Include "magnus_common.pro"; @@ -196,6 +197,9 @@ FunctionSpace{ NameOfConstraint DeltaPhi;} Else If( Flag_Object == 0 ) // if Cylinder only + // In case of the Airfoil, the contraint on Dmdt is omitted + // and replaced by an equation "Dmdt=$newDmdt" + // See the "Resolution" section. {NameOfCoef Dmdt; EntityType GroupsOfNodesOf; NameOfConstraint MassFlowRate;} EndIf @@ -205,7 +209,7 @@ FunctionSpace{ } Formulation{ - {Name LPE; Type FemEquation; + {Name PotentialFlow; Type FemEquation; Quantity{ {Name phi; Type Local; NameOfSpace Vh;} { Name phiCont ; Type Local ; NameOfSpace Vh[phiCont] ; } @@ -228,9 +232,9 @@ Formulation{ } Resolution{ - {Name LPE; + {Name PotentialFlow; System{ - {Name A; NameOfFormulation LPE;} + {Name A; NameOfFormulation PotentialFlow;} } Operation{ InitSolution[A]; @@ -241,9 +245,11 @@ Resolution{ Else // Pseudo-Newton iteration to determine Dmdt in the Airfoil case. - // The run-time variable $newDmdt used in Generate[A] - // whereas $circ, $dmdt and argV are evaluated - // by PostOperation[Trailing] + // The run-time variable $newDmdt is used in Generate[A] + // whereas $circ, $dmdt, $argV and $phiTrailing are evaluated + // by the PostOperation[Trailing]. + // A file containing "PhiTrailing = $phiTrailing;" is after all + // written on disk to be used for formatting the DeleteFile["KJiter.txt"]; @@ -256,7 +262,7 @@ Resolution{ PostOperation[Trailing]; Print[{$syscount, $circ, $dmdt, $argV}, - Format "iter = %3g Circ = %7.5e Dmdt = %7.5e argV = %7.5e"]; + Format "iter = %3g Circ = %5.2f Dmdt = %5.2f argV = %5.2e"]; While[ Norm[ $argV ] > 1e-3 && $syscount < 50] { @@ -273,20 +279,17 @@ Resolution{ Evaluate[$argVp = $argV]; PostOperation[Trailing]; - Print[{$syscount, $circ, $dmdt, $argV}, - Format "iter = %3g Circ = %7.5e Dmdt = %7.5e argV = %7.5e"]; + Print[{$syscount, $circ, $dmdt, $jac, $argV}, + Format "iter = %3g Circ = %5.2f Dmdt = %5.2f jac=%5.2f argV = %5.3e"]; } - // A python script is used to generate a format file "KJ.pos.opt" - // for the color map "KJ.pos" that depends - // on the run-time variable "phiTrailing". - Evaluate[ Python[$phiTrailing]{"generateOptFile.py"} ]; + Print[{$phiTrailing}, Format "PhiTrailing=%g;", File "PhiTrailing.txt"]; EndIf } } } PostProcessing{ - {Name LPE; NameOfFormulation LPE; + {Name PotentialFlow; NameOfFormulation PotentialFlow; Quantity{ {Name phi; Value { Local{ [ {phi} ] ; In Dom_Vh; Jacobian Vol; } } } @@ -296,8 +299,10 @@ PostProcessing{ Local { [ { phiDisc } ] ; In Dom_Vh ; Jacobian Vol ; } } } {Name velocity; Value { Local { [ {d phi} ]; In Dom_Vh; Jacobian Vol; } } } + {Name normVelocity; Value { + Local { [ Norm[{d phi}] ]; In Dom_Vh; Jacobian Vol; } } } {Name pressure; Value { - Local { [ -0.5*rho[]*SquNorm[ {d phi} ]]; + Local { [-0.5*rho[]*SquNorm[ {d phi} ]]; In Dom_Vh; Jacobian Vol; } } } {Name Angle; Value { Local{ [ argVTrail[{d phi}] ]; @@ -331,7 +336,7 @@ PostProcessing{ } PostOperation{ - {Name LPE; NameOfPostProcessing LPE; + {Name PotentialFlow; NameOfPostProcessing PotentialFlow; Operation{ Print[Circ, OnRegion Sur_Cut, File > "output.txt", Color "Ivory", @@ -354,35 +359,41 @@ PostOperation{ //Print[pressure, OnElementsOf Vol_rho, File "p.pos"]; Print[ velocity, OnElementsOf Vol_rho, File "velocity.pos"]; - Echo[ StrCat["l=PostProcessing.NbViews-1;", - "View[l].VectorType = 1;", - "View[l].LineWidth = 2;", - "View[l].ArrowSizeMax = 100;", - "View[l].CenterGlyphs = 1;"], + Echo[ Str["l=PostProcessing.NbViews-1;", + "View[l].VectorType = 1;", + "View[l].LineWidth = 2;", + "View[l].ArrowSizeMax = 100;", + "View[l].CenterGlyphs = 1;"], File "tmp.geo", LastTimeStepOnly] ; - If( Flag_Circ || ( Flag_Object == 0 ) ) - //same condition as in the Resolution section - - Print[phi, OnElementsOf Vol_rho, File "phi.pos"]; - Echo[ StrCat["l=PostProcessing.NbViews-1;", - "View[l].IntervalsType = 3;", - "View[l].NbIso = 40;"], - File "tmp.geo", LastTimeStepOnly] ; + // Stagnation points are points where Norm[{d phi}] is zero + // This is better seen in log scale. + Print[normVelocity, OnElementsOf Vol_rho, File "p.pos"]; + Echo[ Str["l=PostProcessing.NbViews-1;", + "View[l].Name = 'stagnation points';", + "View[l].ScaleType = 2;"], + File "tmp.geo", LastTimeStepOnly] ; - Else - Print[phi, OnElementsOf Vol_rho, File "KJ.pos"]; - Echo[ StrCat["l=PostProcessing.NbViews-1;", - "View[l].IntervalsType = 3;", - "View[l].NbIso = 40;"], - File "tmp.geo", LastTimeStepOnly] ; - EndIf + // Show the isovalue phi=$PhiTrailing + // which should be perpendicular to the airfoil + // if the Kutta condition is fulfilled + Print[phi, OnElementsOf Vol_rho, File "KJ.pos"]; + Echo[Str["Merge 'PhiTrailing.txt';", + "l=PostProcessing.NbViews-1;", + "View[l].Name = 'isovalue phiTrailing';", + "View[l].IntervalsType = 3;", + "Mesh.SurfaceEdges = 0;", + "View[l].RangeType = 2;", + "View[l].CustomMin = PhiTrailing;", + "View[l].CustomMax = 1.001*PhiTrailing;", + "View[l].NbIso = 10;"], + File "tmp.geo", LastTimeStepOnly] ; } } } PostOperation{ // for the Airfoil model - {Name Trailing; NameOfPostProcessing LPE; + {Name Trailing; NameOfPostProcessing PotentialFlow; Operation{ Print[Circ, OnRegion Sur_Cut, File > "KJiter.txt", Format Table, StoreInVariable $circ, @@ -403,8 +414,8 @@ PostOperation{ // for the Airfoil model DefineConstant[ - R_ = {"LPE", Name "GetDP/1ResolutionChoices", Visible 0}, + R_ = {"PotentialFlow", Name "GetDP/1ResolutionChoices", Visible 0}, C_ = {"-solve -pos", Name "GetDP/9ComputeCommand", Visible 0}, - P_ = {"LPE", Name "GetDP/2PostOperationChoices", Visible 0} + P_ = {"PotentialFlow", Name "GetDP/2PostOperationChoices", Visible 0} ];