Adding example of validation of transonic potential flow scheme

README.md

@@ -30,6 +30,14 @@ Unit tests should *not* be uploaded to this repository. Please put them in the `
 - [Two-fluids dam break scenario](https://github.com/KratosMultiphysics/Examples/blob/master/fluid_dynamics/validation/two_fluid_dam_break/README.md)
 - [Two-fluids wave propagation](https://github.com/KratosMultiphysics/Examples/blob/master/fluid_dynamics/validation/two_fluid_wave/README.md)
 
+## Potential Flow
+
+**Use cases**
+- [Naca 0012 compressible - sensitivity analysis](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/use_cases/sensitivity_analysis_compressible/source)
+
+**Validation**
+- [Naca 0012 transonic scheme](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/validation/naca0012_transonic_scheme/source)
+
 ## Structural Mechanics
 
 **Use cases**

potential_flow/README.md

@@ -0,0 +1,7 @@
+## Potential Flow
+
+**Use cases**
+- [Naca 0012 compressible - sensitivity analysis](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/use_cases/sensitivity_analysis_compressible/source)
+
+**Validation**
+- [Naca 0012 transonic scheme](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/validation/naca0012_transonic_scheme/source)

potential_flow/use_cases/README.md

@@ -0,0 +1,5 @@
+# Use Cases
+
+This folder contains the use cases:
+
+- [Naca 0012 compressible - sensitivity analysis](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/use_cases/sensitivity_analysis_compressible/source)

potential_flow/validation/README.md

@@ -0,0 +1,5 @@
+# Validation
+
+This folder contains the validation cases:
+
+- [Naca 0012 transonic scheme](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/validation/naca0012_transonic_scheme/source)

potential_flow/validation/naca0012_transonic_scheme/README.md

@@ -0,0 +1,93 @@
+# Transonic Potential Flow 
+# Naca 0012 airfoil - Mach number = 0.8, Angle of attack = 0º
+
+**Author:** [Marco Antonio Zuñiga Perez](https://github.com/marco1410)
+
+**Kratos version:** 9.5
+
+**Source files:** [Naca 0012 transonic potential flow scheme](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/validation/naca0012_transonic_scheme/source)
+
+## Case Specification
+This is a 2D simulation of a NACA 0012 airfoil under inviscid transonic flow conditions (Ma = 0.8).  An adaptative upwind density scheme (1) is used as artificial compressibility method in order to stabilize the problem in supersonic flow regions.
+
+The flow around the airfoil is computed at zero angle of attack and Mach 0.8 to illustrate the effect of using variable parameters to bias the density. 
+
+The problem geometry consists in a 100[m]x100[m] domain and an airfoil with a chord length of 1[m] at zero degrees angle of attack. The computational domain was meshed with 7920 linear triangular elements (Figure 1).
+
+<p align="center">
+  <img src="data/Airfoil_geometry.png" style="width: 600px;"/>
+  <figcaption align="center"> Figure 1: NACA0012 airfoil geometry and domain. </figcaption align="center">
+</p>
+
+### Boundary Conditions
+The boundary conditions imposed in the far field are:
+
+* Free stream density = 1.225 _kg/m<sup>3</sup>_
+* Angle of attack = 0º
+* Mach infinity = 0.8
+* Wake condition is computed automatically by the app.
+
+### General Solving Settings
+The problem is solved with a finite - element transonic potential flow solver with and embedded wake approach (2). 
+
+Settings:
+
+* Maximum non-linear iterations = 30
+* Convergence criterion = "residual_criterion"
+* Relative and absolute tolerances = 1e-8
+* Non-Linear solver: Newton Raphson with a Line search strategy
+* Linear solver: "sparse_lu"
+
+### Transonic Scheme Settings
+In the current test, for the first case, the upwind factor constant and the critical mach are initialized to 2 and 0.90 to produce strong stabilization over a large portion of the flow. As the solution converges, below of a certain residual tolerance, these parameters are set to 1 and 0.95. These final values ​​were chosen because they are suitable for most cases. The second case is the same but here the parameters are not updated, and in the last case, the parameters are only initialized to 1.0 and 0.95.
+
+In summary, the _perturbation_transonic_ potential flow element is selected and the transonic scheme configurations used were:
+
+* **Case 1:**
+    * Initial critical mach = 0.90
+    * Initial upwind factor constant = 2.0
+    * Update relative residual norm = 1e-3
+    * Target critical mach = 0.95
+    * Target upwind factor constant = 1.0
+
+* **Case 2:**
+    * Initial critical mach = 0.90
+    * Initial upwind factor constant = 2.0
+    * The values are not updated
+
+* **Case 3:**
+    * Initial critical mach = 0.95
+    * Initial upwind factor constant = 1.0
+    * The values are not updated
+
+## Results
+Below three snapshots depicting the Cp distributions over the airfoil skin, Cp distributions over the whole domain and the relative residual norm convergence histogram obteined with three sets of parameters. As it can be observed, the expected symmetric solution and shocks are obtained.
+
+The Figure 2 shows the coefficient of pressure (Cp) distribution along the chord of the airfoil. It is observed how the solution is affected by the different parameters tested. A really sharp solution (blue curve), a smoothed solution (orange curve) and a solution with not-converged shock region solution (green curve). 
+
+<p align="center">
+  <img src="data/Airfoils_Cp_x.png" alt="Effect of varying the upwind parameters during the solution procedure." style="width: 600px;"/>
+  <figcaption align="center"> Figure 2: Effect of varying the upwind parameters during the solution procedure. </figcaption>
+</p>
+
+The Figure 3 shows the Cp distribution over the whole domain for different final values of the parameters. For a given grid density, changing the values of the parameters allows to fine tune the strength and the location of the shock.
+
+<p align="center">
+  <img src="data/Airfoils_Cp_distributions.png" alt="Cp distributions at transonic flow around a NACA0012 airfoil." style="width: 600px;"/>
+  <figcaption align="center"> Figure 3: Pressure distributions around the NACA 0012 airfoil at zero angle of attack and Mach 0.8 for the three cases. </figcaption>
+</p>
+
+The Figure 4 shows the evolution of the relative residual, computed as the L2 norm of the unknown potential vector, as a function of the iteration count. The Figure shows that the algorithm converges quite fast to a solution in which the shock is smeared and displaced upstream. On the other hand, setting the parameters directly to 1 and 0.95 to produce a low amount of bias (green curve) prevents the scheme from converging, and the obtained solution is oscillatory near the shock. Varying the parameters (blue curve) allows the solution procedure to converge and yields to a sharp solution.
+
+<p align="center">
+  <img src="data/residual_norm_histogram.png" alt="Relative residual norm convergence histogram." style="width: 600px;"/>
+  <figcaption align="center"> Figure 4: Transonic flow around a NACA0012 airfoil: Relative residual norm convergence histogram. </figcaption>
+</p>
+
+In practice, the optimal values of the parameters ​​are often case-dependent and sensitive to the grid density, making them difficult to choose a priori. Additionally, changing the grid density has a similar effect on the results.
+
+## References
+(1) Crovato, A. (2020). <em>Steady Transonic Aerodynamic and Aeroelastic Modeling for Preliminary Aircraft Design</em>. [Doctoral thesis, ULiège - Université de Liège]. ORBi-University of Liège. [https://hdl.handle.net/2268/251906] (https://orbi.uliege.be/handle/2268/251906)
+
+(2) López Canalejo, Iñigo Pablo. (2022). <em>A finite-element transonic potential flow solver with an embedded wake approach for aircraft conceptual design</em>. [Doctoral thesis, TUM School of Engineering and Design]. TUM-Technische Universität München. [https://nbn-resolving.de/urn/resolver.pl?urn:nbn:de:bvb:91-diss-20220505-1633175-1-3] (https://mediatum.ub.tum.de/?id=1633175)
+

potential_flow/validation/naca0012_transonic_scheme/source/MainKratos.py

@@ -0,0 +1,109 @@
+import sys
+import time
+import importlib
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import cbook
+
+import KratosMultiphysics
+
+def CreateAnalysisStageWithFlushInstance(cls, global_model, parameters):
+    class AnalysisStageWithFlush(cls):
+
+        def __init__(self, model,project_parameters, flush_frequency=10.0):
+            super().__init__(model,project_parameters)
+            self.flush_frequency = flush_frequency
+            self.last_flush = time.time()
+            sys.stdout.flush()
+
+        def Initialize(self):
+            super().Initialize()
+            sys.stdout.flush()
+
+        def FinalizeSolutionStep(self):
+            super().FinalizeSolutionStep()
+
+            if self.parallel_type == "OpenMP":
+                now = time.time()
+                if now - self.last_flush > self.flush_frequency:
+                    sys.stdout.flush()
+                    self.last_flush = now
+
+    return AnalysisStageWithFlush(global_model, parameters)
+
+if __name__ == "__main__":
+
+    critical_machs_and_upwind_factor_constants_list = []
+    critical_machs_and_upwind_factor_constants_list.append([0.9 , 2.0, 1e-3 , 0.95, 1.0])
+    critical_machs_and_upwind_factor_constants_list.append([0.9 , 2.0, 1e-30, 0.95, 1.0]) # not updated
+    critical_machs_and_upwind_factor_constants_list.append([0.95, 1.0, 1e-30, 0.95, 1.0]) # not updated
+
+    with open("ProjectParameters.json", 'r') as parameter_file:
+        parameters = KratosMultiphysics.Parameters(parameter_file.read())
+
+    analysis_stage_module_name = parameters["analysis_stage"].GetString()
+    analysis_stage_class_name = analysis_stage_module_name.split('.')[-1]
+    analysis_stage_class_name = ''.join(x.title() for x in analysis_stage_class_name.split('_'))
+
+    analysis_stage_module = importlib.import_module(analysis_stage_module_name)
+    analysis_stage_class = getattr(analysis_stage_module, analysis_stage_class_name)
+
+    for id, values in enumerate(critical_machs_and_upwind_factor_constants_list): 
+        KratosMultiphysics.kratos_utilities.DeleteFileIfExisting("residual_per_iteration.txt")
+
+        parameters["solver_settings"]["scheme_settings"]["initial_critical_mach"].SetDouble(values[0])
+        parameters["solver_settings"]["scheme_settings"]["initial_upwind_factor_constant"].SetDouble(values[1])
+
+        parameters["solver_settings"]["scheme_settings"]["update_relative_residual_norm"].SetDouble(values[2]) 
+
+        parameters["solver_settings"]["scheme_settings"]["target_critical_mach"].SetDouble(values[3])
+        parameters["solver_settings"]["scheme_settings"]["target_upwind_factor_constant"].SetDouble(values[4]) 
+
+        output_name = f'{values[0]}_{values[1]}_{values[2]}_{values[3]}_{values[4]}'
+        parameters["output_processes"]["gid_output"][0]["Parameters"]["output_name"].SetString(f'gid_output/{output_name}') 
+        parameters["output_processes"]["vtk_output"][0]["Parameters"]["output_path"].SetString(f'vtk_output/{output_name}') 
+
+        global_model = KratosMultiphysics.Model()
+        simulation = CreateAnalysisStageWithFlushInstance(analysis_stage_class, global_model, parameters)
+        simulation.Run()
+
+        modelpart = global_model["FluidModelPart.Body2D_Body"]
+        fout = open(f'{id}.dat','w')
+        for i,node in enumerate(modelpart.Nodes):
+            x  = node.X0
+            cp = node.GetValue(KratosMultiphysics.PRESSURE_COEFFICIENT)
+            fout.write("%s %s\n" %(x,cp))
+        fout.close()
+
+    x     = np.loadtxt('0.dat', usecols=(0,))
+    cp_0  = np.loadtxt('0.dat', usecols=(1,))
+    cp_1  = np.loadtxt('1.dat', usecols=(1,))
+    cp_2  = np.loadtxt('2.dat', usecols=(1,))
+
+    fig, ax = plt.subplots()
+    fig.set_figwidth(8.0)
+    fig.set_figheight(6.0)
+    # make data
+    ax.plot( x, -cp_0, ".", label='Mc = 0.90 Ufc = 2.0 updated to Mc = 0.95  Ufc = 1.0')
+    ax.plot( x, -cp_1, ".", label='Mc = 0.90 Ufc = 2.0 not updated')
+    ax.plot( x, -cp_2, ".", label='Mc = 0.95 Ufc = 1.0 not updated')
+    ax.grid()
+    # inset Axes....
+    x1, x2, y1, y2 = 0.40, 0.60, -0.1, 1.0  # subregion of the original image
+    axins = ax.inset_axes([0.65, 0.45, 0.3, 0.52])
+    axins.set_xlim(x1, x2)
+    axins.set_ylim(y1, y2)
+
+    axins.plot( x, -cp_0, ".")
+    axins.plot( x, -cp_1, ".")
+    axins.plot( x, -cp_2, ".")
+    axins.grid()
+    ax.indicate_inset_zoom(axins, edgecolor="grey")
+
+    fig = plt.ylabel('Cp')
+    fig = plt.xlabel('x')
+    fig = plt.legend()
+    fig = plt.title('Naca0012 - Alpha = 0.0º Mach = 0.8')
+    fig = plt.tight_layout()
+    fig = plt.savefig("Airfoils_Cp_x.png", dpi=400)
+    fig = plt.close('all')

potential_flow/validation/naca0012_transonic_scheme/source/ProjectParameters.json

@@ -0,0 +1,141 @@
+{
+    "analysis_stage"   : "KratosMultiphysics.CompressiblePotentialFlowApplication.potential_flow_analysis",
+    "problem_data"     : {
+        "problem_name"  : "naca0012_2D",
+        "parallel_type" : "OpenMP",
+        "echo_level"    : 0,
+        "start_time"    : 0.0,
+        "end_time"      : 1.0
+    },
+    "output_processes" : {
+        "gid_output" : [{
+            "python_module" : "gid_output_process",
+            "kratos_module" : "KratosMultiphysics",
+            "process_name"  : "GiDOutputProcess",
+            "Parameters"    : {
+                "model_part_name"        : "FluidModelPart",
+                "postprocess_parameters" : {
+                    "result_file_configuration" : {
+                        "gidpost_flags"               : {
+                            "GiDPostMode"           : "GiD_PostBinary",
+                            "WriteDeformedMeshFlag" : "WriteDeformed",
+                            "WriteConditionsFlag"   : "WriteConditions",
+                            "MultiFileFlag"         : "SingleFile"
+                        },
+                        "file_label"                  : "step",
+                        "output_control_type"         : "step",
+                        "output_interval"             : 1,
+                        "body_output"                 : true,
+                        "node_output"                 : false,
+                        "skin_output"                 : false,
+                        "plane_output"                : [],
+                        "nodal_results"               : ["VELOCITY_POTENTIAL","AUXILIARY_VELOCITY_POTENTIAL"],
+                        "nodal_nonhistorical_results" : ["PRESSURE_COEFFICIENT","VELOCITY","DENSITY"],
+                        "gauss_point_results"         : ["MACH","DENSITY","TEMPERATURE"]
+                    },
+                    "point_data_configuration"  : []
+                },
+                "output_name"            : "gid_output/naca0012_2D"
+            }
+        }]
+        ,
+        "vtk_output" : [{
+            "python_module" : "vtk_output_process",
+            "kratos_module" : "KratosMultiphysics",
+            "process_name"  : "VtkOutputProcess",
+            "help"          : "This process writes postprocessing files for Paraview",
+            "Parameters"    : {
+                "model_part_name"                    : "FluidModelPart",
+                "output_control_type"                : "step",
+                "output_interval"                    : 1,
+                "file_format"                        : "binary",
+                "output_precision"                   : 7,
+                "output_sub_model_parts"             : true,
+                "output_path"                        : "vtk_output",
+                "save_output_files_in_folder"        : true,
+                "nodal_solution_step_data_variables" : ["VELOCITY_POTENTIAL","AUXILIARY_VELOCITY_POTENTIAL"],
+                "nodal_data_value_variables"         : ["PRESSURE_COEFFICIENT","VELOCITY","DENSITY",
+                                                        "MACH","POTENTIAL_JUMP","TRAILING_EDGE","WAKE_DISTANCE","WING_TIP"],
+                "element_data_value_variables"       : ["WAKE","KUTTA"],
+                "gauss_point_variables_in_elements"  : ["PRESSURE_COEFFICIENT","VELOCITY","MACH"],
+                "condition_data_value_variables"     : []
+            }
+        }]
+    },
+    "solver_settings"  : {
+        "model_part_name"          : "FluidModelPart",
+        "domain_size"              : 2,
+        "solver_type"              : "potential_flow",
+        "model_import_settings"    : {
+            "input_type"     : "mdpa",
+            "input_filename" : "naca0012"
+        },
+        "formulation": {
+                    "element_type":"perturbation_transonic"
+        },
+        "scheme_settings"          : {
+            "is_transonic" : true,
+            "initial_critical_mach"         : 0.9,
+            "initial_upwind_factor_constant": 2.0,
+            "target_critical_mach"          : 0.95,
+            "target_upwind_factor_constant" : 1.0,
+            "update_relative_residual_norm" : 1e-3,
+            "mach_number_squared_limit"     : 3.0
+        },
+        "maximum_iterations"       : 30,
+        "echo_level"               : 1,
+        "relative_tolerance"               : 1e-8,
+        "absolute_tolerance"               : 1e-8,
+        "solving_strategy_settings"         : {
+            "type" : "line_search",
+            "advanced_settings": {
+                "first_alpha_value"         : 0.5,
+                "second_alpha_value"        : 1.0,
+                "min_alpha"                 : 0.5,
+                "max_alpha"                 : 1.0,
+                "line_search_tolerance"     : 0.5,
+                "max_line_search_iterations": 5
+            }
+        },
+        "linear_solver_settings"  : {
+                "solver_type"             : "LinearSolversApplication.sparse_lu"
+        },
+        "volume_model_part_name"   : "Parts_Parts_Auto1",
+        "skin_parts"               : ["PotentialWallCondition2D_Far_field_Auto1","Body2D_Body"],
+        "no_skin_parts"            : [],
+        "reform_dofs_at_each_step" : false,
+        "auxiliary_variables_list" : []
+    },
+    "processes"        : {
+        "boundary_conditions_process_list" : [{
+            "python_module" : "apply_far_field_process",
+            "kratos_module" : "KratosMultiphysics.CompressiblePotentialFlowApplication",
+            "process_name"  : "FarFieldProcess",
+            "Parameters"    : {
+                "model_part_name" : "FluidModelPart.PotentialWallCondition2D_Far_field_Auto1",
+                "free_stream_density"   : 1.225,
+                "angle_of_attack"       : 0.0,
+                "mach_infinity"         : 0.8,
+                "perturbation_field"    : true 
+            }
+        }
+        ,{
+            "python_module" : "define_wake_process_2d",
+            "kratos_module" : "KratosMultiphysics.CompressiblePotentialFlowApplication",
+            "process_name"  : "DefineWakeProcess2D",
+            "Parameters"    : {
+                "model_part_name" : "FluidModelPart.Body2D_Body"
+            }
+        }
+        ,{
+            "python_module" : "compute_nodal_value_process",
+            "kratos_module" : "KratosMultiphysics.CompressiblePotentialFlowApplication",
+            "process_name"  : "ComputeNodalValueProcess",
+            "Parameters"    : {
+                "model_part_name" : "FluidModelPart",
+                "elemental_variables_list_to_project": ["VELOCITY", "PRESSURE_COEFFICIENT"]
+            }
+        }],
+        "auxiliar_process_list"            : []
+    }
+}