diff --git a/ntrfc/math/vectorcalc.py b/ntrfc/math/vectorcalc.py
index 1042176ac273d5f7cc2385f3992edbde96db0300..0daef1b5bed4671897f7c5e894d34846cb3aa2b5 100644
--- a/ntrfc/math/vectorcalc.py
+++ b/ntrfc/math/vectorcalc.py
@@ -209,6 +209,16 @@ def ellipsoidVol(sig):
def vecProjection(direction, vector):
+ """
+ Calculate the projection of a vector onto a direction vector.
+
+ Parameters:
+ direction (list or numpy array): The direction vector onto which the projection will be calculated.
+ vector (list or numpy array): The vector to be projected onto the direction vector.
+
+ Returns:
+ projection (numpy array): A vector representing the projection of the input vector onto the direction vector.
+ """
unitDir = vecDir(direction)
return np.dot(vector, unitDir) * unitDir
diff --git a/ntrfc/meshquality/nondimensionals.py b/ntrfc/meshquality/nondimensionals.py
index e9be7ab83af38b5a6cb1ab3a454b8f6e41cb57de..638ff60c59b93e6eb4b3ace635fb6843208ef734 100644
--- a/ntrfc/meshquality/nondimensionals.py
+++ b/ntrfc/meshquality/nondimensionals.py
@@ -24,8 +24,8 @@ from ntrfc.math.vectorcalc import vecAbs, vecAbs_list, vecProjection, vecAngle,
def cellDirections(cellUMean, wallNorm):
x = unitvec(cellUMean) # mainDirection
- z = unitvec(wallNorm) # tangential
- y = unitvec(np.cross(x, z)) # spanwise
+ y = unitvec(wallNorm) # tangential
+ z = unitvec(np.cross(x, y)) # spanwise
return np.array([x, y, z])
@@ -67,20 +67,28 @@ def cellSpans(solutionMesh, calcFrom):
return cell_spans
-def getWalluTaus(solutionMesh, surfaceMesh, mu_0, rhofieldname, velfieldname):
- surfaceNormals = surfaceMesh.extract_surface().compute_normals()
- cell_centers = solutionMesh["cellCenters"]
- relatedpointonwallids = [surfaceNormals.find_closest_point(i) for i in cell_centers]
- relatedpointonwall = [surfaceNormals.points[i] for i in relatedpointonwallids]
- relatedpointonwall_samples = pv.PolyData(relatedpointonwall).sample(surfaceMesh)
- graduwall = relatedpointonwall_samples["gradient"].reshape(relatedpointonwall_samples.number_of_points, 3, 3)
- flowdirections = unitvec_list(solutionMesh[velfieldname])
- gradUyWs = vecAbs_list(
- [np.dot(flowdirections[i], graduwall[i]) for i in range(relatedpointonwall_samples.number_of_points)])
- rhoWs = relatedpointonwall_samples[rhofieldname]
+def getWalluTaus(solutionMesh, mu_0, rhofieldname, velfieldname):
+ """
+ Calculate the wall shear stress at various points on a surface in a CFD simulation.
+
+ Parameters:
+ solutionMesh (pv.PolyData): A mesh that represents the solution of the CFD simulation.
+ surfaceMesh (pv.PolyData): A mesh that represents the surface on which the wall shear stress is to be calculated.
+ mu_0 (float): The dynamic viscosity of the fluid.
+ rhofieldname (str): The name of the field that contains the density of the fluid.
+ velfieldname (str): The name of the field that contains the velocity of the fluid.
+
+ Returns:
+ u_taus (np.ndarray): An array containing the wall shear stress at each point on the surface.
+ """
- tauWs = gradUyWs * mu_0 * rhoWs
+
+ graduwall = solutionMesh[f"grad_{velfieldname}"].reshape(solutionMesh.number_of_cells,3,3)
+ wallnormals = solutionMesh["wallNormal"]
+ dudy = vecAbs_list([np.dot(graduwall[i], wallnormals[i]) for i in range(solutionMesh.number_of_cells)])
+ rhoWs = solutionMesh[rhofieldname]
+ tauWs = dudy * mu_0
u_taus = (tauWs / rhoWs) ** 0.5
return u_taus
@@ -115,6 +123,13 @@ def constructWallMesh(surfaces):
# #
#################################################################################
+
+def compute_scalar_gradient(mesh,arrayname):
+ mesh = mesh.compute_derivative(scalars=arrayname)
+ mesh[f"grad_{arrayname}"] = mesh["gradient"].reshape((mesh.number_of_cells,3,3))
+ return mesh
+
+
def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield, mu_0):
"""
:param volmesh: pyvista-vtk object
@@ -127,31 +142,27 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
print("constructing surfacemesh from wall meshes ...")
surfaceMesh = constructWallMesh(surfaces)
- surfaceMeshcopy = surfaceMesh.copy()
- surfaceMeshcopy.clear_data()
-
- print("preparing processData from meshes")
-
- volmesh = volmesh.compute_derivative(scalars=use_velfield)
- volmesh_walladjacent = volmesh.extract_cells(volmesh.find_containing_cell(volmesh.extract_surface().points))
- surfaceMeshcopy = surfaceMeshcopy.sample(volmesh_walladjacent)
- surfaceMesh["gradient"] = surfaceMeshcopy["gradient"]
- volmesh_walladjacent["cellCenters"] = volmesh_walladjacent.cell_centers().points
-
print("calculating wall-normal vectors...")
surfacenormals_surface = surfaceMesh.extract_surface().compute_normals()
+ print("finding walladjacent cells")
+ volmesh = compute_scalar_gradient(volmesh,use_velfield)
+ walladjacentids = volmesh.find_containing_cell(surfacenormals_surface.points)
+ volmesh_walladjacent = volmesh.extract_cells(walladjacentids)
+
+ volmesh_walladjacent["cellCenters"] = volmesh_walladjacent.cell_centers().points
volmesh_walladjacent["wallNormal"] = [
surfacenormals_surface.point_data["Normals"][surfacenormals_surface.find_closest_point(i)]
for i in volmesh_walladjacent.points]
- print("calculating cell spans from WallNormals and CellEdges...")
+
+ print("calculating cell spans from WallNormals...")
spanS = cellSpans(volmesh_walladjacent, use_velfield)
volmesh_walladjacent["xSpan"] = np.array([i[0] for i in spanS]) # calculate cell span in flow direction
volmesh_walladjacent["ySpan"] = np.array([i[1] for i in spanS]) # calculate cell span in wall normal direction
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"] = getWalluTaus(volmesh_walladjacent, surfaceMesh, mu_0, use_rhofield, use_velfield)
+ volmesh_walladjacent["uTaus"] = getWalluTaus(volmesh_walladjacent, mu_0, use_rhofield, use_velfield)
print("calculating grid spacing")
gridSpacings = gridSpacing(mu_0, volmesh_walladjacent)
@@ -161,3 +172,4 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
volmesh_walladjacent["DeltaZPlus"] = gridSpacings[2]
return volmesh_walladjacent
+
diff --git a/tests/meshquality/test_ntrfc_nondimensionals.py b/tests/meshquality/test_ntrfc_nondimensionals.py
index 8ad6acb89b56122c3853216b505d72f27ac459bc..5fd477a3e8f3657ed41710704b2c262e7af11ce5 100644
--- a/tests/meshquality/test_ntrfc_nondimensionals.py
+++ b/tests/meshquality/test_ntrfc_nondimensionals.py
@@ -20,9 +20,9 @@ def test_calc_dimensionless_gridspacing():
velocity = 2
gradient = velocity / height
# analytical solution
- span_x = length / (nodes - 1)
- span_y = width / (nodes - 1)
- span_z = height / (nodes - 1)
+ span_x = width / (nodes - 1)
+ span_y = height / (nodes - 1)
+ span_z = length / (nodes - 1)
wallshearstress = mu_0 * gradient
wallshearvelocity = np.sqrt(wallshearstress / rho)
deltaxplus = wallshearvelocity * span_x / mu_0
@@ -36,7 +36,13 @@ def test_calc_dimensionless_gridspacing():
grid = pv.StructuredGrid(x, y, z)
# scale the mesh
grid.points /= nodes - 1
- grid.points *= np.array([length, height, width])
+ grid.points *= np.array([width, height, length])
+
+ # define velocityfield
+ bounds = grid.bounds
+ min_z = bounds[4]
+ grid["U"] = [gradient * (grid.cell_centers().points[::,1][i]-min_z) * np.array([0,0,1]) for i in range(grid.number_of_cells)]
+ grid["rho"] = np.ones(grid.number_of_cells)
# extract surface
surface = grid.extract_surface()
facecellids = surface.surface_indices()
@@ -50,19 +56,16 @@ def test_calc_dimensionless_gridspacing():
upperwallids.append(faceid)
lowerwall = surface.extract_cells(lowerwallids)
upperwall = surface.extract_cells(upperwallids)
- # define velocityfield
- grid["U"] = np.array([gradient, 0, 0]) * grid.cell_centers().points
- lowerwall["U"] = np.array([gradient, 0, 0]) * lowerwall.cell_centers().points
- upperwall["U"] = np.array([gradient, 0, 0]) * upperwall.cell_centers().points
- grid["rho"] = np.ones(grid.number_of_cells)
+ lowerwall["U"] = [gradient * (lowerwall.cell_centers().points[::,1][i]-min_z) * np.array([0,0,1]) for i in range(lowerwall.number_of_cells)]
+ upperwall["U"] = [gradient * (upperwall.cell_centers().points[::,1][i]-min_z) * np.array([0,0,1]) for i in range(upperwall.number_of_cells)]
lowerwall["rho"] = np.ones(lowerwall.number_of_cells) * rho
upperwall["rho"] = np.ones(upperwall.number_of_cells) * rho
dimless_gridspacings = calc_dimensionless_gridspacing(grid, [lowerwall, upperwall], "U", "rho", mu_0)
- assert all(np.isclose(deltaxplus, dimless_gridspacings[
+ assert all(np.isclose(deltazplus, dimless_gridspacings[
"DeltaXPlus"], rtol=0.05)), "calc_dimensionelss_gridspcing in x direction not accurate"
assert all(np.isclose(deltayplus, dimless_gridspacings[
"DeltaYPlus"], rtol=0.05)), "calc_dimensionelss_gridspcing in y direction not accurate"
- assert all(np.isclose(deltazplus, dimless_gridspacings[
+ assert all(np.isclose(deltaxplus, dimless_gridspacings[
"DeltaZPlus"], rtol=0.05)), "calc_dimensionelss_gridspcing in z direction not accurate"
runtest(height=1, length=1, width=1)