Develop

incomp_navierstokes/streamwise_periodic/pipeSlice_3d/pipeslice.geo

@@ -0,0 +1,113 @@
+//-------------------------------------------------------------------------------------//
+// T. Kattmann, 13.05.2018, 3D Butterfly mesh in a circular pipe
+// Create the mesh by calling this geo file with 'gmsh <this>.geo'.
+//-------------------------------------------------------------------------------------//
+
+// Evoque Meshing Algorithm?
+Do_Meshing= 1; // 0=false, 1=true
+// Write Mesh files in .su2 format
+Write_mesh= 1; // 0=false, 1=true
+
+//Geometric inputs, ch: channel, Pin center is origin
+Radius= 0.5e-2; // Pipe Radius
+InnerBox= Radius/2; // Distance to the inner Block of the butterfly mesh
+
+//Mesh inputs
+gridsize = 0.1; // unimportant once everything is structured
+
+//ch_box 
+Nbox = 30; // Inner Box points in x direction
+
+Ncircu = 30; // Outer ring circu. points
+Rcircu = 0.9; // Spacing towards wall
+
+sqrtTwo = Cos(45*Pi/180);
+
+//-------------------------------------------------------------------------------------//
+//Points
+// Inner Box
+Point(1) = {-InnerBox, -InnerBox, 0, gridsize};
+Point(2) = {-InnerBox, InnerBox, 0, gridsize};
+Point(3) = {InnerBox, InnerBox, 0, gridsize};
+Point(4) = {InnerBox, -InnerBox, 0, gridsize};
+
+// Outer Ring
+Point(5) = {-Radius*sqrtTwo, -Radius*sqrtTwo, 0, gridsize};
+Point(6) = {-Radius*sqrtTwo, Radius*sqrtTwo, 0, gridsize};
+Point(7) = {Radius*sqrtTwo, Radius*sqrtTwo, 0, gridsize};
+Point(8) = {Radius*sqrtTwo, -Radius*sqrtTwo, 0, gridsize};
+
+Point(9) = {0,0,0,gridsize}; // Helper Point for circles
+
+//-------------------------------------------------------------------------------------//
+//Lines
+//Inner Box (clockwise)
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+//Walls (clockwise)
+Circle(5) = {5, 9, 6};
+Circle(6) = {6, 9, 7};
+Circle(7) = {7, 9, 8};
+Circle(8) = {8, 9, 5};
+
+//Connecting lines (outward facing)
+Line(9) = {1, 5};
+Line(10) = {2, 6};
+Line(11) = {3, 7};
+Line(12) = {4, 8};
+
+//-------------------------------------------------------------------------------------//
+//Lineloops and surfaces
+// Inner Box (clockwise)
+Line Loop(1) = {1,2,3,4}; Plane Surface(1) = {1};
+
+// Ring sections (clockwise starting at 9 o'clock)
+Line Loop(2) = {5, -10, -1, 9}; Plane Surface(2) = {2};
+Line Loop(3) = {10, 6, -11, -2}; Plane Surface(3) = {3};
+Line Loop(4) = {-3, 11, 7, -12}; Plane Surface(4) = {4};
+Line Loop(5) = {12, 8, -9, -4}; Plane Surface(5) = {5};
+
+//make structured mesh with transfinite lines
+//radial
+Transfinite Line{1, 2, 3, 4, 5, 6, 7, 8} = Nbox;
+//circumferential
+Transfinite Line{9, 10, 11, 12} = Ncircu Using Progression Rcircu;
+
+Transfinite Surface{1,2,3,4,5};
+Recombine Surface{1,2,3,4,5};
+
+//Extrude 1 mesh layer
+Extrude {0, 0, 0.0005} {
+    Surface{1}; Surface{2}; Surface{3}; Surface{4}; Surface{5};
+    Layers{1};
+    Recombine; 
+}
+Coherence;
+
+//Physical groups made with GUI
+Physical Surface("inlet") = {4, 1, 5, 3, 2};
+Physical Surface("outlet") = {100, 122, 56, 78, 34};
+Physical Surface("wall") = {69, 95, 113, 43};
+Physical Volume("fluid") = {1, 2, 3, 4, 5};
+
+// ----------------------------------------------------------------------------------- //
+// Meshing
+Transfinite Surface "*";
+Recombine Surface "*";
+
+If (Do_Meshing == 1)
+    Mesh 1; Mesh 2; Mesh 3;
+EndIf
+
+// ----------------------------------------------------------------------------------- //
+// Write .su2 meshfile
+If (Write_mesh == 1)
+
+    Mesh.Format = 42; // .su2 mesh format, 
+    Save "pipe1cell3D.su2";
+
+EndIf
+

incomp_navierstokes/streamwise_periodic/pipeSlice_3d/plots.py

@@ -0,0 +1,162 @@
+#! /usr/bin/python3
+# --------------------------------------------------------------------------- #
+# Kattmann, 16.07.2019
+# This python script provides some plots to test the match between analytical
+# and simulated solution for a 3D circular laminar pipe flow, either from 
+# streamwise periodic simulation or the outlet of a suitable long pipe.
+#
+# requires: surface_flow.dat (SURFACE_TECPLOT_ASCII) in current directory
+#
+# output: plots (opened in separate window, not saved)
+#
+# optional: which plots to show
+showLineplot       = True
+show2Dsurfaceplots = False
+show3Dplots        = False
+# --------------------------------------------------------------------------- #
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+from mpl_toolkits.mplot3d import Axes3D
+from scipy.spatial import Delaunay
+from scipy.interpolate import LinearNDInterpolator
+
+# --------------------------------------------------------------------------- #
+# Import data from surface_flow.dat into pandas dataframe
+data = pd.read_csv("surface_flow.dat", nrows=4264, skiprows=3, sep='\t', header=None)
+x = data[0][:]
+y = data[1][:]
+vel_z = data[6][:]
+
+# Create Delaunay surface triangulation from scatterd dataset
+points2D = np.vstack([x,y]).T
+tri = Delaunay(points2D)
+
+# --------------------------------------------------------------------------- #
+# Create analytic solution vector on the same points as the imported data
+dynanmic_vsicosity = 1.8e-5
+pressure_drop = 1e-3
+domain_length = 5e-4
+radius = 5e-3
+
+analytic_sol = -1/(4*dynanmic_vsicosity) * (-pressure_drop/domain_length) * \
+    (radius**2 - ((x**2 + y**2)**(0.5))**2 )
+
+perc_devi_from_anal = abs(analytic_sol - vel_z) / max(analytic_sol) * 100
+maxvel = max(abs(perc_devi_from_anal)) # get absolute maximum of dataset
+
+# --------------------------------------------------------------------------- #
+# Plot velocity on line from domain midpoint to wall
+if showLineplot:
+    plt.close()
+
+    # interpolator (ip) for simulated and analytical dataset
+    ip_sim = LinearNDInterpolator(tri, vel_z)
+    ip_ana = LinearNDInterpolator(tri, analytic_sol)
+    # line (which lies on the x-axis) where values will be interpolated
+    n_sample_points = 30
+    x_line = np.linspace(0, radius-5e-6, n_sample_points)
+    y_line = np.zeros(n_sample_points)
+    ip_pos = np.vstack((x_line,y_line)).T
+
+    ax = plt.axes()
+    plt.plot(ip_sim(ip_pos), x_line, color='b', marker='', linestyle='--', linewidth=3, label='simulated')
+    plt.plot(ip_ana(ip_pos), x_line, color='r', marker='', linestyle=':' , linewidth=3, label='analytical')
+    plt.legend()
+    plt.title('Velocity profile: analytic vs simulated (interpolated values)')
+    plt.xlabel('velocity [m/s]')
+    plt.ylabel('radius [m]')
+    ax.set_aspect(aspect=max(ip_sim(ip_pos)) / max(x_line)) # make plot square
+    plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+    plt.grid(True, linestyle='--')
+    plt.show()
+
+# --------------------------------------------------------------------------- #
+# Plot various 2D surface plots of sim. and analy. data
+if show2Dsurfaceplots:
+    plt.close()
+
+    fig, ax = plt.subplots(2,2)
+
+    # 1. analytical solution
+    ax_tmp = ax[0,0]
+
+    tcf = ax_tmp.tricontourf(x, y, abs(analytic_sol))
+    ax_tmp.scatter(x,y, s=0.1, color='black', marker='.')
+
+    ax_tmp.set_title("Analytical solution")
+    ax_tmp.set_aspect('equal')
+    fig.colorbar(tcf, ax=ax_tmp)
+
+    # 2. simulated solution
+    ax_tmp = ax[1,0]
+
+    tcf = ax_tmp.tricontourf(x, y, vel_z)
+    ax_tmp.scatter(x,y, s=0.1, color='black', marker='.')
+
+    ax_tmp.set_title("Simulated solution")
+    ax_tmp.set_aspect('equal')
+    fig.colorbar(tcf, ax=ax_tmp)
+
+    # 3. absolute value deviation between analytic and simulated
+    ax_tmp = ax[0,1]
+
+    tcf = ax_tmp.tricontourf(x, y, abs(analytic_sol-vel_z), cmap=plt.cm.Greys)
+    ax_tmp.scatter(x,y, s=0.1, color='black', marker='.')
+
+    ax_tmp.set_title("abs(analytic-simulated)")
+    ax_tmp.set_aspect('equal')
+    fig.colorbar(tcf, ax=ax_tmp)
+
+    # 4. percentual deviation scaled by the maximal value
+    ax_tmp = ax[1,1]
+
+    tcf = ax_tmp.tricontourf(x, y, perc_devi_from_anal, cmap=plt.cm.Greys, vmin=0.0, vmax=maxvel)
+    ax_tmp.scatter(x,y, s=0.1, color='black', marker='.')
+
+    ax_tmp.set_title("abs(analytic-simulated) / max(analytic) * 100")
+    ax_tmp.set_aspect('equal')
+    fig.colorbar(tcf, ax=ax_tmp)
+
+    plt.show()
+
+# --------------------------------------------------------------------------- #
+if show3Dplots:
+    # Plot 3D surfaces of sim. and analy. data
+    plt.close()
+
+    # Scatter plot deviation
+    fig = plt.figure()
+    ax = fig.gca(projection='3d')
+
+    ax.scatter(x, y, perc_devi_from_anal)
+    ax.set_xlabel('x [m]')
+    ax.set_ylabel('y [m]')
+    ax.set_zlabel('z-Velocity deviation [%]')
+
+    plt.show()
+
+    # Surface plot deviation
+    fig = plt.figure()
+    ax = fig.gca(projection='3d')
+
+    surf = ax.plot_trisurf(x, y, perc_devi_from_anal, triangles=tri.simplices, cmap='jet', linewidth=0)
+    ax.set_xlabel('x [m]')
+    ax.set_ylabel('y [m]')
+    ax.set_zlabel('z-Velocity deviation [%]')
+    fig.colorbar(surf)
+
+    plt.show()
+
+    # Surface plot of velocity
+    fig = plt.figure()
+    ax = fig.gca(projection='3d')
+
+    surf = ax.plot_trisurf(x, y, vel_z, triangles=tri.simplices, cmap='jet', linewidth=0)
+    ax.set_xlabel('x [m]')
+    ax.set_ylabel('y [m]')
+    ax.set_zlabel('z-Velocity [m/s]')
+    fig.colorbar(surf)
+
+    plt.show()