looks promissing

ntrfc/meshquality/nondimensionals.py

@@ -9,6 +9,7 @@ from itertools import combinations
 
 import numpy as np
 import pyvista as pv
+import vtk
 
 from ntrfc.math.vectorcalc import vecAbs, vecAbs_list, vecProjection, unitvec
 
@@ -48,7 +49,7 @@ def cellSpans(solutionMesh, calcFrom):
     return cell_spans
 
 
-def get_wall_shear_stress(mesh, dynamic_viscosity, density_fieldname, velocity_fieldname):
+def get_wall_shear_stress_velocity(mesh, dynamic_viscosity, density_fieldname, velocity_fieldname):
     """
     Calculate the wall shear stress velocity at various points on a surface in a CFD simulation.
 
@@ -59,16 +60,18 @@ def get_wall_shear_stress(mesh, dynamic_viscosity, density_fieldname, velocity_f
     velocity_fieldname (str): The name of the field that contains the velocity of the fluid.
 
     Returns:
-    wall_shear_stress_velocity (np.ndarray): An array containing the velocity at which a fluid layer adjacent to the surface would need to move in order to experience the same shear stress as the actual fluid layer in contact with the surface.
+    wall_shear_stress_velocity (np.ndarray): An array containing the velocity at which a fluid layer adjacent to
+    the surface would need to move in order to experience the same shear stress as the actual fluid layer in contact
+    with the surface.
     """
+
     grad_velocity = mesh[f"grad_{velocity_fieldname}"].reshape(mesh.number_of_cells, 3, 3)
     wall_normals = mesh["wallNormal"]
     velocity_gradient_normal = vecAbs_list([np.dot(grad_velocity[i], wall_normals[i]) for i in range(mesh.number_of_cells)])
     fluid_density = mesh[density_fieldname]
-    wall_shear_stress = velocity_gradient_normal * dynamic_viscosity
-    wall_shear_stress_velocity = (wall_shear_stress / fluid_density) ** 0.5
+    tau_w = fluid_density*dynamic_viscosity*velocity_gradient_normal
 
-    return wall_shear_stress_velocity
+    return tau_w
 
 
 def constructWallMesh(surfaces):
@@ -99,7 +102,8 @@ def calc_dimensionless_yplus(volmesh, surfaces, use_velfield, use_rhofield, mu_0
     distcompute = pv.PolyData(volmesh_walladjacent["cellCenters"]).compute_implicit_distance(
         surfaceMesh.extract_surface())
     volmesh_walladjacent["cellCentersWallDistance"] = distcompute["implicit_distance"]
-    volmesh_walladjacent["uTaus"] = get_wall_shear_stress(volmesh_walladjacent, mu_0, use_rhofield, use_velfield)
+    volmesh_walladjacent["uTaus"] = get_wall_shear_stress_velocity(volmesh_walladjacent, mu_0, use_rhofield,
+                                                                   use_velfield)
     YPlus = np.abs(volmesh_walladjacent["cellCentersWallDistance"]) * volmesh_walladjacent["uTaus"] / mu_0
     volmesh_walladjacent["yPlus"] = YPlus
     return volmesh_walladjacent
@@ -121,7 +125,17 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
     surfacenormals_surface = surfaceMesh.extract_surface().compute_normals()
 
     print("finding walladjacent cells")
-    volmesh = compute_scalar_gradient(volmesh, use_velfield)
+    merged_mesh = surfaceMesh + volmesh
+
+    # Compute cell-based derivative of a vector field
+    derivatives = compute_scalar_gradient(merged_mesh,"U")
+
+    # Extract the volumetric mesh of the derivative-field
+    cell_types = merged_mesh.celltypes
+    volumetric_cells = np.where(cell_types == vtk.VTK_HEXAHEDRON)[0]
+
+    # Extract the derivative-field
+    volmesh = derivatives.extract_cells(volumetric_cells)
 
     walladjacentids = volmesh.find_containing_cell(surfacenormals_surface.cell_centers().points)
     volmesh_walladjacent = volmesh.extract_cells(walladjacentids)
@@ -138,7 +152,9 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
     volmesh_walladjacent["zSpan"] = np.array([i[2] for i in spanS])  # calculate cell span in span direction
 
     print("calculating wall-shear and friction-velocity")
-    volmesh_walladjacent["uTaus"] = get_wall_shear_stress(volmesh_walladjacent, mu_0, use_rhofield, use_velfield)
+    volmesh_walladjacent["uTaus"] = get_wall_shear_stress_velocity(volmesh_walladjacent, mu_0, use_rhofield,
+                                                                   use_velfield)
+
     uTaus = volmesh_walladjacent["uTaus"]
     volmesh_walladjacent["DeltaXPlus"] = volmesh_walladjacent[
                                              "xSpan"] * uTaus / mu_0  # calculate cell span in flow direction

tests/meshquality/test_ntrfc_nondimensionals.py

@@ -213,7 +213,7 @@ def test_calc_dimensionless_gridspacing():
 
 
 def test_compute_wallshearstress():
-    from ntrfc.meshquality.nondimensionals import get_wall_shear_stress
+    from ntrfc.meshquality.nondimensionals import get_wall_shear_stress_velocity
     # 11**3 nodes are enough
     nodes = 11
     # needs to be something simple
@@ -271,7 +271,7 @@ def test_compute_wallshearstress():
         surfacenormals_surface.point_data["Normals"][surfacenormals_surface.find_closest_point(i)]
         for i in volmesh_walladjacent.points]
 
-    utaus = get_wall_shear_stress(volmesh_walladjacent, mu_0, "rho", "U")
+    utaus = get_wall_shear_stress_velocity(volmesh_walladjacent, mu_0, "rho", "U")
 
     assert np.all(np.isclose(utaus, wallshearvelocity))