Merge pull request #154 from KratosMultiphysics/potentail_fllow/add_incompresible_example

[POTENTIAL] Add example for incompressible potential flow

Author: juancamarotti

potential_flow/README.md

@@ -4,4 +4,5 @@
 - [Naca 0012 compressible - sensitivity analysis](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/use_cases/sensitivity_analysis_compressible/README.md)
 
 **Validation**
+- [Naca 0012 incompressible potential flow](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/validation/naca0012_potential_incompressible_flow/README.md)
 - [Naca 0012 transonic scheme](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/validation/naca0012_transonic_scheme/README.md)

potential_flow/validation/naca0012_potential_incompressible_flow/README.md

@@ -0,0 +1,60 @@
+# Incompressible Potential Flow 
+# Naca 0012 airfoil - Mach number = 0.0588, Angle of attack = 5º
+
+**Author:** [Juan Ignacio Camarotti](https://github.com/juancamarotti)
+
+**Kratos version:** 10.4
+
+**Source files:** [Naca 0012 incompressible potential flow](https://github.com/KratosMultiphysics/Examples/tree/master/potential_flow/validation/naca0012_potential_incompressible_flow/source)
+
+## Case Specification
+This example presents a two-dimensional simulation of inviscid incompressible flow around a NACA 0012 airfoil using the potential flow solver available in the Kratos CompressiblePotentialFlowApplication.
+
+The flow around the airfoil is computed at AOA=5° and Mach 0.0588. Although the formulation is incompressible, the Mach number is provided to define the free-stream velocity.
+
+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 = 5º
+* Mach infinity = 0.0588
+* Wake condition is computed automatically by the app
+* Slip boundary condition on the airfoil surface, enforcing zero normal velocity
+
+## Results
+
+Two snapshots of the pressure coefficient ($C_p$) are shown below: one depicting the distribution along the airfoil surface and another showing the pressure field over the entire computational domain.
+
+Figure 2 presents the pressure coefficient ($C_p$) distribution along the chord of the airfoil.
+
+<p align="center">
+  <img src="data/Cp_distribution_x.png" alt="Pressure coefficient distribution along the airfoil chord." style="width: 600px;"/>
+  <figcaption align="center"> Figure 2: Pressure coefficient (Cp) distribution along the chord of the NACA 0012 airfoil. </figcaption>
+</p>
+
+Figure 3 shows the spatial distribution of the pressure coefficient over the entire computational domain.
+
+<p align="center">
+  <img src="data/Cp_distribution_airfoil.png" alt="Pressure coefficient distribution around the NACA 0012 airfoil." style="width: 600px;"/>
+  <figcaption align="center"> Figure 3: Pressure coefficient (Cp) distribution in the computational domain around the NACA 0012 airfoil. </figcaption>
+</p>
+
+The aerodynamic coefficients obtained from the simulation can be compared with reference data using the lift curve shown in Figure 4.
+
+<p align="center">
+  <img src="data/naca0012_lift_curve.png" alt="NACA 0012 Lift vs AOA Curve" style="width: 600px;"/>
+  <figcaption align="center"> Figure 4: Reference lift coefficient (Cl) as a function of the angle of attack for the NACA 0012 airfoil. </figcaption>
+</p>
+
+For an angle of attack of $\alpha = 5^\circ$, the computed lift coefficient is $C_l = 0.569$. This result shows very good agreement with the reference data presented in Figure 4, which is used here for validation purposes.
+
+## References
+(1) https://aerospaceweb.org/question/airfoils/q0259c.shtml
+

potential_flow/validation/naca0012_potential_incompressible_flow/data/Airfoil_geometry.png

@@ -0,0 +1,107 @@
+import sys
+import time
+import importlib
+import numpy as np
+import matplotlib.pyplot as plt
+
+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__":
+
+    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)
+
+    # Create model and run simulation
+    global_model = KratosMultiphysics.Model()
+    simulation = CreateAnalysisStageWithFlushInstance(analysis_stage_class, global_model, parameters)
+    simulation.Run()
+
+    # Extract Cp
+    modelpart = global_model["FluidModelPart.Body2D_Body"]
+
+    upper_x = []
+    upper_cp = []
+    lower_x = []
+    lower_cp = []
+
+    for node in modelpart.Nodes:
+        x = node.X0
+        y = node.Y0
+        cp = node.GetValue(KratosMultiphysics.PRESSURE_COEFFICIENT)
+
+        if y >= 0.0:
+            upper_x.append(x)
+            upper_cp.append(cp)
+        else:
+            lower_x.append(x)
+            lower_cp.append(cp)
+
+    # convert to numpy
+    upper_x = np.array(upper_x)
+    upper_cp = np.array(upper_cp)
+    lower_x = np.array(lower_x)
+    lower_cp = np.array(lower_cp)
+
+    # sort along chord
+    upper_order = np.argsort(upper_x)
+    lower_order = np.argsort(lower_x)
+
+    upper_x = upper_x[upper_order]
+    upper_cp = upper_cp[upper_order]
+
+    lower_x = lower_x[lower_order]
+    lower_cp = lower_cp[lower_order]
+
+    # Plot
+    fig, ax = plt.subplots()
+    fig.set_figwidth(8)
+    fig.set_figheight(6)
+
+    ax.plot(upper_x, -upper_cp, "o", label="Upper surface")
+    ax.plot(lower_x, -lower_cp, "o", label="Lower surface")
+
+    ax.set_xlabel("x/c")
+    ax.set_ylabel("-Cp")
+    ax.set_title("NACA0012 Potential Flow")
+    ax.grid()
+
+    # aerodynamic convention
+    ax.invert_yaxis()
+
+    ax.legend()
+    plt.tight_layout()
+    plt.savefig("Cp_distribution.png", dpi=400)
+    plt.show()

potential_flow/validation/naca0012_potential_incompressible_flow/data/Cp_distribution_airfoil.png

@@ -0,0 +1,155 @@
+{
+    "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": "incompressible"
+        },
+        "scheme_settings"          : {
+            "is_transonic" : false
+            ,
+            "initial_critical_mach"         : 0.0,
+            "initial_upwind_factor_constant": 0.0,
+            "target_critical_mach"          : 0.0,
+            "target_upwind_factor_constant" : 0.0,
+            "update_relative_residual_norm" : 0.0,
+            "mach_number_squared_limit"     : 0.0
+        },
+        "maximum_iterations"       : 30,
+        "compute_reactions": true,
+        "echo_level"               : 2,
+        "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.pardiso_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_and_wake_process",
+            "kratos_module" : "KratosMultiphysics.CompressiblePotentialFlowApplication",
+            "process_name"  : "ApplyFarFieldAndWakeProcess",
+            "Parameters"    : {
+                "model_part_name" : "FluidModelPart.PotentialWallCondition2D_Far_field_Auto1",
+                "angle_of_attack_units": "degrees",
+                "free_stream_density"   : 1.225,
+                "mach_infinity": 0.0588235,
+                "angle_of_attack"       : 5.0,
+                "perturbation_field"    : false,
+                "define_wake": true,
+                "wake_type": "Operations.KratosMultiphysics.CompressiblePotentialFlowApplication.Define2DWakeOperation", 
+                "wake_parameters"    : {
+                    "body_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"]
+            }
+        }
+        ,
+        {
+            "python_module" : "compute_lift_process",
+            "kratos_module" : "KratosMultiphysics.CompressiblePotentialFlowApplication",
+            "process_name"  : "ComputeLiftProcess",
+            "Parameters"    : {
+                "model_part_name": "FluidModelPart.Body2D_Body",
+                "moment_reference_point" : [0.0,0.0,0.0],
+                "mean_aerodynamic_chord" : 1.0,
+                "far_field_model_part_name": "FluidModelPart.PotentialWallCondition2D_Far_field_Auto1",
+                "trailing_edge_model_part_name": "FluidModelPart.Body2D_Body",
+                "is_infinite_wing": false
+            }
+        }
+    ],
+        "auxiliar_process_list"            : []
+    }
+}