NUMERICAL AND EXPERIMENTAL INVESTIGATION OF WALL SHEAR STRESS
GENERATED BY A BERNOULLI PAD

By

Anshul Singh Tomar

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

Mechanical Engineering—Doctor of Philosophy

2024

ABSTRACT

Bernoulli pads can create a significant normal force on an object without contact, which allows

them to be traditionally used for non-contact pick-and-place operations in industry. In addition to

the normal force, the pad produces shear forces, which can be utilized in cleaning a workpiece

without contact. The motivation for the present work has been to understand the flow physics of the

Bernoulli pad such that they can be employed for non-contact biofouling mitigation of ship hulls.

Numerical investigations have shown that the shear stress distribution generated by the action of the

Bernoulli pad on the workpiece is concentrated and results in maximum shear stress very close to

the neck of the pad. The maximum value of wall shear stress is an important metric for determining

the cleaning efficacy of the Bernoulli pad. We use numerical simulations over a range of parameter

space to develop a relationship between the inlet fluid power and the maximum shear stress obtained

on the workpiece. To increase the shear force distribution, we explore the possibility of adding

mechanical power to the system in addition to the fluid power. The flow field between the Bernoulli

pad and the workpiece typically involves a recirculation region and transition between laminar and

turbulent flow. The maximum shear stress occurs in the vicinity of the recirculation region and to

gain confidence in the numerical solver’s ability to estimate these stresses accurately, experiments

were conducted with a hot-film sensor.

The main contributions of this work are as follows: a direct relationship is obtained between

the maximum shear stress on the workpiece and inlet fluid power using dimensional analysis. A

relationship between the maximum shear stress and the inlet Reynolds number is also obtained,

and implications of these scaling relationships are studied. A direct relationship between the inlet

fluid power and the shear losses motivates us to explore other methods of providing power to the

system with the objective of increasing shear forces and thereby improving cleaning efficacy. We

numerically investigate a Bernoulli pad in which additional mechanical power is added by rotating

the pad. This additional power increases both the normal and shear forces on the workpiece for the

same inlet fluid power. In the context of the rotating Bernoulli pad, it is found that for a given normal

attractive force, a stable equilibrium configuration can exist for two different mass flow rates, with

the higher mass flow rate resulting in a higher stiffness of the flow field. This phenomenon has

not been reported in the literature. The wall shear stress distribution, obtained using numerical

simulations, is validated using experiments for the first time. A constant temperature anemometer

is used with a hot-film sensor and water as the working fluid; the sensor is calibrated using a

fully developed channel flow. An experimental setup is designed to calibrate and later measure

the wall shear stress in a Bernoulli pad assembly. The maximum wall shear stress is observed

very close to the neck of the pad due to flow constriction and separation; the hot-film experiments

accurately capture the magnitude of the maximum shear stress and its location. This provides us

with confidence in the numerical solver, which can be used to optimize the Bernoulli pad design to

improve its cleaning efficacy.

Copyright by
ANSHUL SINGH TOMAR
2024

Dedicated to my wife Nikita Talukdar

v

ACKNOWLEDGEMENTS

I want to express my sincere gratitude to my thesis advisor, Dr. Ranjan Mukherjee, for his invaluable

guidance, support, and encouragement throughout the process of completing this thesis. His

expertise, patience, and unwavering commitment have been instrumental in shaping the direction

of my research and helping me navigate the challenges along the way. I am deeply thankful for his

insightful feedback, constructive criticism, and dedication to excellence, which have significantly

enriched the quality of this work. His mentorship has enhanced my academic skills and inspired me

to strive for excellence in all my endeavors. Furthermore, I am grateful to him for his mentorship

beyond academia, his willingness to share their knowledge and wisdom, and for serving as a role

model of integrity and professionalism. I would also like to thank Dr. Aren Hellum for his guidance

and contribution to this research work. Aren was always available whenever I needed suggestions

and feedback. I will always be grateful for his patience; it has greatly helped me.

I would also like to thank Dr. Ricardo M. Alvarez for his continuous support throughout this

journey. I have learned a lot from him, and his insights into my research work have helped me

greatly. I would also like to take this opportunity to thank my thesis committee member, Dr. Andre

Benard, for providing me with valuable suggestions for shaping this thesis.

My journey as a Ph.D. student was a roller coaster ride, but it had always had a constant catalyst

for pacifying me - my wife, Nikita. I can not thank you enough for trusting me and making me

what I am today. I can have a good or bad day, but when I am home, you make me forget about

everything with that smile. Most days, we will go out for walks, and you will listen to my work and

boost my morale without even saying a word. I know you don’t need this, but thank you for being

my best friend and a superb life partner.

I would also like to thank my parents and sisters. Your constant support and well-being ensured

that I am always motivated to perform better. Your sacrifices have inspired me in more ways than

one can imagine. You guys are the best. I am blessed to have wonderful in-laws who cheered

for me throughout my journey. I am thankful for my wonderful friends here: Sai, Nikhil, Vipesh,

Arjun, Revanth, and Mudita. Thank you all for being wonderful listeners and kind to me. This

vi

journey was impossible without the support I got from my lab mates - Aakash, Shaede, Sanders,

and Robert.

At last, I would like to acknowledge the Office of Naval Research, ONR Grant N00014220-1-

2170 for funding this research work.

vii

CHAPTER 1

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

1

TABLE OF CONTENTS

CHAPTER 2

Introduction .

RELATIONSHIP BETWEEN INLET FLUID POWER AND WALL
SHEAR STRESS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.

6
6
.
6
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
. 16
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. 18
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

.
.
2.1
.
. .
2.2 Methodology .
.
2.3 Scaling Analysis .
.
.
2.4 Results .
.
2.5 Discussion .
.
.
2.6 Concluding Remarks

.
.

.
.

.

CHAPTER 3

.
.
. .

.
. .

Introduction .

FLOW PHYSICS OF A ROTATING BERNOULLI PAD . . . . . . . .

.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 24
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.1
3.2 Approach .
. 26
3.3 Flow physics of a rotating pad . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.4 Flow physics of a rotating Bernoulli Pad . . . . . . . . . . . . . . . . . . . . . 33
. 39
. . . . .
3.5 Rotating Bernoulli Pad: Variation of Normal Force with Gap Height
3.6 Rotating Bernoulli Pad: Energy Balance and Advantages . . . . . . . . . . .
. 45
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.7 Conclusion .

.

.

.

.

CHAPTER 4

Introduction .

WALL SHEAR MEASUREMENT: HOT-FILM ANEMOMETRY AND
MODEL TESTING . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 55
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.1
.
.
. 56
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Background . .
. 60
4.3 Calibration of Hot-Film Sensor . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Bernoulli Pad Experimental Setup . . . . . . . . . . . . . . . . . . . . . . .
. 64
4.5 Wall Shear Stress Measurement . . . . . . . . . . . . . . . . . . . . . . . . . . 69
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.6 Conclusion .

.
.
. .

.

.

.

.

CHAPTER 5

5.1 Summary . .
5.2 Future Scope

SUMMARY AND FUTURE WORK . . . . . . . . . . . . . . . . . . . 77
. 77
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

. .
.
.

. .
.
.

BIBLIOGRAPHY .

.

.

.

.

. .

. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

APPENDIX

. .

. .

. .

. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

viii

CHAPTER 1

INTRODUCTION

Bernoulli pads or Bernoulli grippers are used widely to grip objects or surfaces without physically

contacting them. In its simplest implementation, this device consists of an axial jet surrounded by

a surface that is parallel to the exit of the jet - see Figure 1.1. Other implementations exist, which

use a simple (Brun and Melkote (2009)) or complex (Wagner et al. (2008)) center body to redirect

the axial flow as it exits the jet. Because of their ability to manipulate delicate parts in a clean

environment, Bernoulli pads and similar technologies (Paivanas and Hassan (1981)) see widespread

use in the semiconductor industry, and a rich patent literature focused on this application has also

been developed (Logue (1970), Frey (1999), and McIlwraith and Christie (2003)). Because these

devices are able to exert force to grip soft and pliable workpieces, they have also been explored

in medical (Erturk and Erzincanli (2020)), apparel (Fantoni et al. (2014)), and meat processing

(Misimi et al. (2016)) applications.

Figure 1.1 A typical Bernoulli gripper used in pick and place industries (P.C.:
https://youtu.be/phSlc8B4kFo)

The direction and magnitude of the normal force produced by the pad depend on the distance

between the pad and the opposite surface.

In the limit where the pad contacts the surface, the

repulsive force is equal to the product of the jet’s area and feed pressure, whereas at very large

1

Figure 1.2 A Bernoulli pad showing the inlet and outlet flow parameters.

distances, the repulsive force scales with the jet momentum. An attractive force is produced at

some intermediate distance, which can be understood, with an important caveat, in the context of

the Bernoulli equation. By combining the Bernoulli equation with mass conservation, we find an

expression for the local pressure for the geometry in Figure 1.2. For small values of the product

ℎ𝑟, a net negative gage pressure results, which can produce a net attractive force integrated over

the pad area. Because the force is repulsive in the limits of large and small ℎ and attractive for

intermediate values, the system has two equilibrium points, one stable and the other unstable (Li

and Kagawa (2014), Kamensky et al. (2019)).

𝑝(𝑟) = 𝑝(0) −

𝜌

2

(cid:19)

(cid:18) 𝜋𝑑2𝑢𝑖𝑛
8𝜋ℎ𝑟

(1.1)

Assuming a sufficient flow rate for the applied force, the system will operate at its stable

equilibrium height ℎ𝑒𝑞, which dictates the distance between the pad and the opposing surface. This

height is important in a variety of applications (Wagner et al. (2008), Olsson and Williams (1969))

and is also linked to the amount of fluid power required to operate the device (Kamensky et al.

(2019)).

Kamensky et al. (2019) produced a scaling relationship between operating power and equilib-

rium height, the derivation of which indicates that work added to the system is primarily required

to overcome shear at the wall. This wall shear is the reason that the inviscid analysis in (1.1) is not

fully correct (Wark and Foss (1984)). For most applications, the wall shear is a by-product of the

2

Figure 1.3 The colonization of microorganisms, plants, algae, or animals affecting the ships life
(Pic Credit: https://safety4sea.com/cm-understanding-marine-biofouling-how-anti-fouling-
systems-prevent-growth/)

normal force that supports the workpiece. Although the effect of pressure distribution on delicate

(Brun and Melkote (2009)) and flexible (Dini et al. (2009)) workpieces has been examined, shear is

less well-studied. Typically, deformation or destruction of the opposing surface due to shear forces

is not a substantial concern because the maximum normal stresses are much larger.

The present work is motivated by a grooming application in which a Bernoulli pad is moved

around the submerged hull of a marine vessel to remove biofouling. In this application, the pad

keeps itself close to the surface of the vessel using the attractive normal forces generated and

simultaneously uses the wall shear forces for biofouling mitigation (Kamensky et al. (2020)).

This method of cleaning is consistent with the operation of commercially available fouling-release

coatings (Hu et al. (2020)), which use boundary-layer shear while the vessel is underway to slough

fouling organisms. Biofouling colonization on a ship’s hull is most problematic in port, where

these forces are not present- see Figure 1.3. For the grooming device to keep the vessel clean, it

is necessary for the local shear stress imposed at the hull to be large enough to dislodge fouling

organisms, a value which depends on the organism and settling time (Menesses et al. (2017)).

Therefore, the maximum shear stress created by the flow field is a parameter of interest for this

application.

This work aims to understand the flow physics involving a Bernoulli pad and employ them as

3

a non-contact shear-based grooming tool. A power law relation between the power available at

the inlet and the maximum shear stress produced on the workpiece is obtained. This relation will

enable us to explore the other ways of providing inlet power to the Bernoulli pad. The study aims

to increase the wall shear stress value for better grooming efficacy of biofouling mitigation and

increase normal suction force for better adhering capacity. The CFD simulations in the present

work are carried out using commercial solvers such as ANSYS Fluent and COMSOL. The flow

domain is relatively simple, but the flow field involves a circulation region and a transition from a

laminar to a turbulence regime. Hence, to accurately resolve the flow field, the numerical solver

used for the simulations requires validation from experiments. We have designed an experimental

setup to calibrate the hot-film sensor in the present work. Further, we have designed a setup that

uses a Bernoulli pad and proximally kept workpiece to measure the wall shear stress using the

calibrated sensor.

In chapter 2, we simulate a large number of Bernoulli pad geometries using water as the working

fluid. The methods used in these simulations are validated by a Particle Tracking Velocimetry (PTV)

experiment performed on one of the geometries. These simulations indicate that the maximum shear

stress increases as the equilibrium height decreases and the inlet fluid power increases. We present a

scaling analysis to obtain a power-law relation for the inlet fluid power based on the maximum shear

stress and equilibrium height. The CFD data also indicate a power-law relationship between the

maximum shear stress and the inlet Reynolds number. The empirical form found for the relationship

between maximum shear stress and inlet power in Appendix also indicates a direct relationship

between the inlet Reynolds number and the length scale which characterizes the dissipation of

power by wall shear.

In chapter 3, we provide rotational motion to the Bernoulli pad. The solver used for the

numerical simulation is validated against the Particle Tracking Velocimetry (PTV) experiment

performed for one of the geometries and the angular speed of zero rad/s. The working fluid is water

for all the simulations. The effect of rotation on the computational domain is presented with and

without inlet fluid power. The effect of rotation on the pressure profile on the workpiece and force

4

variation with gap height is discussed. Later, work analysis is carried out, and the work available

at the inlet, outlet, and mechanical work is tabulated. The effect of angular speed on the force and

shear stress on the workpiece is also presented.

In chapter 4, hot-film experiments are discussed to measure the wall shear stress obtained on

the workpiece. We have used a fully developed flow through a channel to calibrate the hot-film

probe. The design of the channel for achieving an early, fully developed flow is presented. Later,

the sensor is mounted on the channel based on the location of the fully developed regime. The

readings from the pressure transducers connected with the channel and hot-film probes are used to

obtain the calibration coefficient.

In chapter 5, we have summarized the findings of the present work, and a brief discussion is

presented about the future work to be carried out. The future scope involves studies with cavitation

effects, rigid curved Bernoulli pads and experiments for rotating Bernoulli pad.

5

CHAPTER 2

RELATIONSHIP BETWEEN INLET FLUID POWER AND WALL SHEAR STRESS

2.1

Introduction

A large part of this chapter was used to publish an article in the Flow journal Tomar et al.

(2022). The motivation for the present work is to employ the Bernoulli pad for non-contact

biofouling mitigation on the ship hulls. The axial flow from the inlet of the Bernoulli pad strikes

the workpiece and passes radially in the outward direction through a constricted passage. The gap

between the pad and the workpiece is very small, and due to the sharp curvature of the pad, a

recirculation region is formed very close to the neck of the pad. The velocity gradients under the

recirculation regions are very high, and therefore, below the belly of the recirculation region, the

maximum value of wall shear stress is observed.

The maximum value of wall shear stress is an essential parameter for the application of non-

contact biofouling mitigation of ship hulls. Hence, the focus of this work revolves around the

quantification of wall shear stress. The wall shear stress generated on the workpiece depends upon

the gap height, pad dimensions, and the inlet fluid power provided to the system. The inlet fluid

power is easily measurable for real-life problems. Therefore, in this chapter, we have used CFD

simulations to obtain a power-law relationship between the inlet fluid power and the maximum

shear stress generated on the workpiece.

2.2 Methodology

2.2.1 Computational domain

The Bernoulli pad is axially symmetric, and therefore, the fluid dynamics is investigated using a

two-dimensional axisymmetric computational model; this significantly reduces the computational

time. The computational domain is defined by the inner diameter 𝑑, outer diameter 𝐷, and gap

height ℎ - see Figure 1.2. An axial flow is provided at the inlet; assuming incompressibility, the flow

parameters include the inlet velocity 𝑢𝑖𝑛, the inlet pressure 𝑝𝑖𝑛, and fluid density 𝜌. After impinging

the wall, the fluid flows radially outward through the gap; the pressure at the outlet is atmospheric,

6

𝑝𝑎𝑡𝑚. The surface roughness of the pad and the wall are 𝜀 𝑝 and 𝜀𝑤, respectively, and no-slip

boundary conditions are imposed on both surfaces. Following the work by Kamensky Kamensky

et al. (2019), a four equation-based transition-SST model is used for turbulence modeling. As

discussed therein, there exists a stable equilibrium gap height ℎ𝑒𝑞, where the net normal force

on the wall is zero. In this study, the same model is used to determine the inlet velocities that

correspond to different configurations of stable equilibrium gap, ℎ𝑒𝑞.

It should be noted that

although ℎ𝑒𝑞 is an imposed condition for the CFD model, in a physical system, ℎ𝑒𝑞 is the result

of the imposed flow conditions and pad geometry. In this study, we investigate the relationship

between inlet fluid power and maximum wall shear stress at ℎ𝑒𝑞.

2.2.2 Equilibrium height

The constricted flow of the fluid in the gap introduces a net normal force 𝐹 on the wall; the

magnitude and direction of this force depend on the gap height ℎ. As discussed by Kamensky et al.

Kamensky et al. (2019), the pad and the wall repel each other (𝐹 is positive) when the gap is small

Figure 2.1 Variation of force 𝐹 with gap height ℎ for a Bernoulli pad with 𝑑 = 31.75 mm and
𝐷 = 300 mm for six different inlet velocities 𝑢𝑖𝑛 = {1.0, 1.5, 2.0, 2.5, 2.847, 3.187} m/s.

7

Figure 2.2 (a) Experimental setup for PTV and (b) CAD model of Bernoulli pad (Kamensky,
Kamensky (2020))

and attract (𝐹 is negative) when the gap is increased beyond a certain value. There exists a critical

gap height for which the net normal force on the wall is zero (𝐹 = 0); this gap height is referred

to as the equilibrium height ℎ𝑒𝑞. The plot in Figure 2.1 shows the variation of 𝐹 with ℎ for six

different values of inlet velocity 𝑢𝑖𝑛 for a specific Bernoulli pad. It is evident that 𝐹 is positive for

small values of ℎ and becomes negative when ℎ is increased; for each inlet velocity, it is possible

to find a value of ℎ = ℎ𝑒𝑞 where the net normal force is equal to zero.

In this chapter, we investigate the relationship between inlet fluid power and maximum wall

shear stress; this relationship will be investigated at the stable equilibrium point where ℎ = ℎ𝑒𝑞.

For a given Bernoulli pad, the value of ℎ𝑒𝑞 depends completely on the inlet velocity. To reduce

the computational time associated with continuous re-meshing of the computational domain due to

change in ℎ, we keep ℎ constant and vary 𝑢𝑖𝑛 till we reach the equilibrium point, i.e., 𝐹 = 0. This

allows us to determine the inlet velocity for a desired ℎ𝑒𝑞.

2.2.3 Validation with PTV

The CFD simulations were carried out using an element size of 6 × 10−5 m; this resulted in

61, 302 elements and a maximum 𝑦+ value of 7.5. A smaller element size of 5 × 10−5 m resulted

in 75, 315 elements but the percentage change in the maximum wall normal and maximum wall

shear stresses (0.11% and 0.87%, respectively) did not justify the additional computational cost.

8

Figure 2.3 Comparison of PTV and CFD data: plot of nondimensional radial velocity at 𝑧 = ℎ/2
with respect to nondimensional radial position

To gain confidence in our CFD results, we compare them with experimental results obtained using

3D Particle Tracking Velocimetry (PTV). The experimental setup consists of four ultra high-speed

cameras (Phantom v2512) placed in an arc-like pattern; these cameras focus on a volume of

rectangular cross-section between the pad and tank wall. A dual-head laser is placed between the

four cameras; it is used to illuminate the particles in the flow - see Figure 2.2 (a).

A CAD model of the Bernoulli pad showing the volume in which the particles are imaged is

shown in Figure 2.2 (b). The camera view is masked to this region to reduce the disk space and

memory requirement. The region includes a portion of axial flow through the stem and radial flow

extending beyond half of the outer radius of the Bernoulli pad. Further details of the experimental

setup can be found in Kamensky Kamensky (2020)).

The PTV data is used to determine the ensemble-averaged velocity on an 𝑟-𝑧 plane of the flow

domain (refer to Figure 1.2) and is compared with CFD results in Figure 2.3 for 𝑧 = ℎ/2: The radial

velocity was nondimensionalized with respect to its maximum value 𝑢𝑚𝑎𝑥, and the radial position

9

was nondimensionalized with respect to the outer radius of the pad. There is a good agreement

between the velocity profiles obtained using CFD simulations and PTV experiments; the location

of the maximum peak is also observed at the same value of 𝑟 (𝑟 = 16.62 mm), which is very near

to the entrance of the fluid in the radial direction (𝑟 = 𝑑/2 = 15.875 mm). This provides us with

confidence to proceed with the analysis of CFD data.

2.2.4 CFD data

The nominal parameter values used for the CFD simulations are

𝑑 = 31.75 𝑚𝑚, 𝐷 = 200 𝑚𝑚,

ℎ𝑒𝑞 = 1.5 𝑚𝑚,

𝜀 𝑝 = 𝜀𝑤 = 0.04 𝑚𝑚

CFD simulations are carried out for a Bernoulli pad with the nominal dimensions using ANSYS

Fluent. The radial velocity and total pressure contours obtained are shown in Figure 2.4 and 2.5.

The gap between the pad and the workpiece is very small; therefore, to provide a better view of the

contours, the contour plots in Figure 2.4 and 2.5 are presented on a magnified and smaller portion

of the complete computational domain. The radial velocity very close to the neck of the pad is

Figure 2.4 Radial velocity contour for the nominal pad

10

negative, which represents the recirculation region, and the maximum velocity is observed below

the recirculation bubble.

The total pressure is maximum at the inlet due to the inlet fluid power provided to the system,

which, after striking the workpiece, reduces in the impinging region. After the flow enters the

gap between the pad and the workpiece, the flow area increases, and therefore, the flow velocity

reduces. This, in turn, implies that the pressure magnitude will increase in the gap. The pressure

boundary condition at the outlet is atmospheric, which implies that the pressure at the neck of the

workpiece is negative. The negative pressure at the neck of the workpiece is responsible for the

suction force generation on the workpiece.

A total of 45 design points were generated using three values of 𝑑 (0.5, 0.707, and 1 times the

nominal value), five different values of 𝐷 (0.5, 0.707, 1, 1.414 and 2 times the nominal value) and

three different values of ℎ𝑒𝑞 = {1.4, 1.5, 1.6} mm. The value of 𝜀 𝑝 and 𝜀𝑤 used is consistent with

Material. A brief discussion of the hydraulic roughness can be found in the Discussion section.

Figure 2.5 Total pressure contour for the nominal pad

11

Figure 2.6 Plot of inlet fluid power versus maximum wall shear stress

The inlet fluid power (cid:164)𝑊𝑖𝑛, defined in Kamensky et al. Kamensky et al. (2019) as

(cid:164)𝑊𝑖𝑛 = 2𝜋

∫ 𝑑/2

0

(cid:20) 1
2

𝜌 [𝑈 (𝑟)]2 + 𝑝(𝑟)

(cid:21)

𝑈 (𝑟) 𝑟𝑑𝑟

and the maximum shear stress on the wall 𝜏𝑤,𝑚𝑎𝑥 are obtained for all these design points; they are

plotted in Figure 2.6 in log-log scale. The 45 data points are grouped into nine sets of five points.

Each set of five points is joined by straight line segments and corresponds to the five different values

of 𝐷 for specific values of 𝑑 and ℎ𝑒𝑞; the direction in which 𝐷 increases is shown in the figure.

The data points shown with ∗, ◦, and △ symbols correspond to specific values of ℎ𝑒𝑞, shown in the

figure.

For clarity, a magnified image of a portion of the plot is shown inset; the top three lines (dashed)

correspond to 𝑑 = 31.75 mm, the middle three lines (solid) correspond to 𝑑 = 31.75×0.707 = 22.44

mm, and the bottom three lines (dotted) correspond to 𝑑 = 31.75 × 0.5 = 15.875 mm. It is evident

that maximum wall shear stress varies linearly with inlet fluid power for a given value of 𝑑 and ℎ𝑒𝑞;

furthermore, as expected, the maximum wall shear stress increases as ℎ𝑒𝑞 decreases for the same

12

level of inlet power.

2.3 Scaling Analysis

In a physical system like the one considered in this study, the gap between the pad and wall

is not imposed. Rather, the system reaches the equilibrium gap ℎ𝑒𝑞 as the result of the imposed

pad geometry (inner and outer diameters, 𝑑 and 𝐷 respectively), and inlet flow conditions (inlet

velocity 𝑢𝑖𝑛, characteristic length 𝑑, and fluid physical properties). Also as a result of those

imposed conditions, the flow will exhibit a maximum wall-shear 𝜏𝑤,𝑚𝑎𝑥 somewhere along the

radius of the gap, and consequently will demand a certain amount of power (cid:164)𝑊𝑖𝑛 to sustain the inlet

flow conditions.

With this in mind, the independent variables governing the flow through a Bernoulli pad in

equilibrium conditions are the geometric parameters 𝑑 and 𝐷 and the inlet flow parameters 𝜌, 𝜇

and 𝑢𝑖𝑛, because those are the only variables controlled by the experimentalist; and in consequence

(cid:164)𝑊𝑖𝑛, ℎ𝑒𝑞, and 𝜏𝑤,𝑚𝑎𝑥 are the dependent variables because they result from the imposed geometric

conditions and inlet flow conditions. Therefore, the dependent variables can be expressed as a

function of the independent variables as follows:

(cid:164)𝑊𝑖𝑛 = 𝑓 (𝜇, 𝑢𝑖𝑛, 𝐷, 𝑑, 𝜌) ,

ℎ𝑒𝑞 = 𝑔 (𝜇, 𝑢𝑖𝑛, 𝐷, 𝑑, 𝜌) ,

𝜏𝑤,𝑚𝑎𝑥 = ℎ (𝜇, 𝑢𝑖𝑛, 𝐷, 𝑑, 𝜌)

(2.1)

After choosing the inlet diameter 𝑑, dynamic viscosity 𝜇, and inlet velocity 𝑢𝑖𝑛 as repeating

variables, a dimensional analysis over each one of the above functional forms provide the following

set of dimensionless parameters:

(cid:101)(cid:164)𝑊 𝑖𝑛 =

(cid:164)𝑊𝑖𝑛
𝑖𝑛 𝑑
𝜇 𝑢2

,

ℎ∗ =

ℎ𝑒𝑞
𝑑

,

𝜏𝑤,𝑚𝑎𝑥 =
(cid:101)

𝜏𝑤,𝑚𝑎𝑥 𝑑
𝜇 𝑢𝑖𝑛

,

𝐷∗ =

𝐷
𝑑

,

𝑅𝑒𝑖𝑛 =

𝜌 𝑢𝑖𝑛 𝑑
𝜇

(2.2)

The nondimensional maximum shear stress (cid:101)

𝜏𝑤,𝑚𝑎𝑥 is the ratio between the maximum wall-shear

stress and a combination of parameters that together also have the dimension of stress: 𝜇 𝑢𝑖𝑛/𝑑.

This combination of parameters is however representative of a shear scale in the inlet pipe rather

13

than in the pad gap. A more appropriate shear scaling would be of the form 𝜇 𝑢𝑐/ℎ𝑒𝑞, where 𝑢𝑐

is the characteristic velocity at the inlet of the pad gap (since the maximum shear occurs in this

region), and can be determined from the volumetric flow rate (cid:164)𝑉 and the inlet cross-section of the

pad gap as follows:

(cid:164)𝑉 ≜ (cid:16) 𝜋
4

𝑑2(cid:17)

𝑢𝑖𝑛

⇒

𝑢𝑐 =

(cid:164)𝑉
𝜋𝑑ℎ𝑒𝑞

=

1
4

·

𝑢𝑖𝑛 𝑑
ℎ𝑒𝑞

Hence, the shear scaling for the inlet region of the pad gap is:

𝜇

𝑢𝑐
ℎ𝑒𝑞

=

1
4

·

𝜇𝑢𝑖𝑛 𝑑
ℎ2
𝑒𝑞

∼

𝜇𝑢𝑖𝑛 𝑑
ℎ2
𝑒𝑞

(2.3)

(2.4)

The factor 1/4 can be ignored because it will be ultimately absorbed by a fitting coefficient. The

nondimensional shear stress (cid:101)

𝜏𝑤,𝑚𝑎𝑥 can now be manipulated based on this shear scaling as follows:

𝜏𝑤,𝑚𝑎𝑥 =
(cid:101)

=

𝜏𝑤,𝑚𝑎𝑥 𝑑
𝜇 𝑢𝑖𝑛

𝜏𝑤,𝑚𝑎𝑥 𝑑
𝜇 𝑢𝑖𝑛 (𝑑/ℎ2

𝜏𝑤,𝑚𝑎𝑥
𝜇 𝑢𝑖𝑛 (𝑑/ℎ2
𝜏𝑤,𝑚𝑎𝑥 is the combination of a new expression for dimensionless

· [ℎ∗]−2

𝑑
ℎ2
𝑒𝑞

(2.5)

𝑒𝑞)

𝑒𝑞)

=

·

Equation 2.5 shows that (cid:101)

maximum shear and ℎ∗, but ℎ∗ is one of the dependent nondimensional variables. To eliminate this

redundancy, (cid:101)

𝜏𝑤,𝑚𝑎𝑥 is replaced by 𝜏∗

𝑤,𝑚𝑎𝑥, which is defined as follows:

𝜏∗
𝑤,𝑚𝑎𝑥 =

𝜏𝑤,𝑚𝑎𝑥
𝜇 𝑢𝑖𝑛 (𝑑/ℎ2

𝑒𝑞)

(2.6)

The set of dimensionless parameters in (2.2) can now be modified to:

(cid:101)(cid:164)𝑊 𝑖𝑛 =

(cid:164)𝑊𝑖𝑛
𝑖𝑛 𝑑
𝜇 𝑢2

,

ℎ∗ =

ℎ𝑒𝑞
𝑑

,

𝜏∗
𝑤,𝑚𝑎𝑥 =

𝜏𝑤,𝑚𝑎𝑥
𝜇 𝑢𝑖𝑛 (𝑑/ℎ2

𝑒𝑞)

,

𝐷∗ =

𝐷
𝑑

,

𝑅𝑒𝑖𝑛 =

𝜌 𝑢𝑖𝑛 𝑑
𝜇

(2.7)

This set of parameters reduces the physical system to three dimensionless functions of the form:

(cid:101)(cid:164)𝑊 𝑖𝑛 = Φ (𝐷∗, 𝑅𝑒𝑖𝑛) ,

ℎ∗ = Γ (𝐷∗, 𝑅𝑒𝑖𝑛) ,

𝑤,𝑚𝑎𝑥 = Ψ (𝐷∗, 𝑅𝑒𝑖𝑛)
𝜏∗

(2.8)

14

The dimensionless system described by (2.8) has two independent variables (𝐷∗ and 𝑅𝑒𝑖𝑛) as

opposed to the five independent variables (𝜇, 𝑢𝑖𝑛, 𝐷, 𝑑, 𝜌) of the dimensional system in (2.1). While

there are no analytical solutions available to determine the exact form of the relationships in (2.8),

it is possible to use regression analysis to determine how the independent variables predict the

dependent variables. Figure 2.6 suggests that the inlet power (cid:164)𝑊𝑖𝑛 is linearly related to the maximum

shear stress 𝜏𝑤,𝑚𝑎𝑥 in a log-log space. This suggests that the relationships in (2.8) are also linear in

a log-log space:

ln (cid:101)(cid:164)𝑊 𝑖𝑛 = 𝐶11 ln 𝐷∗ + 𝐶12 ln 𝑅𝑒𝑖𝑛 + 𝐶13

ln ℎ∗ = 𝐶21 ln 𝐷∗ + 𝐶22 ln 𝑅𝑒𝑖𝑛 + 𝐶23

ln 𝜏∗

𝑤,𝑚𝑎𝑥 = 𝐶31 ln 𝐷∗ + 𝐶32 ln 𝑅𝑒𝑖𝑛 + 𝐶33

(2.9a)

(2.9b)

(2.9c)

Equations (2.9a), (2.9b) and (2.9c) can be used to eliminate 𝐷∗ and 𝑅𝑒𝑖𝑛 and obtain the

relationship between the dependent variables:

ln (cid:101)(cid:164)𝑊 𝑖𝑛 = 𝜅1 ln ℎ∗ + 𝜅2 ln 𝜏∗

𝑤,𝑚𝑎𝑥 + 𝜅3

(2.10)

where the coefficients 𝜅1, 𝜅2 and 𝜅3 can be obtained using the relations

𝜅1 =

𝜅2 =

𝜅3 =

𝐶12𝐶31 − 𝐶11𝐶32
𝐶22𝐶31 − 𝐶21𝐶32
𝐶11𝐶22 − 𝐶12𝐶21
𝐶22𝐶31 − 𝐶21𝐶32
𝐶13𝐶22𝐶31 − 𝐶12𝐶23𝐶31 − 𝐶13𝐶21𝐶32 + 𝐶11𝐶23𝐶32 + 𝐶12𝐶21𝐶33 − 𝐶11𝐶22𝐶33
𝐶22𝐶31 − 𝐶21𝐶32

(2.11a)

(2.11b)

(2.11c)

If the assumption of linearity in log-log space is true, the linear model in (2.10) can also be

used for a direct fit to the data, and the resulting values of the coefficients 𝜅𝑖, 𝑖 = 1, 2, 3, should

be consistent with the values determined from equations (2.11a), (2.11b) and (2.11c). This latter

approach can be viewed as the indirect fit.

15

2.4 Results

A direct fit of the linear model in (2.10) to the CFD data results in the following values of the

coefficients:

(cid:104)

𝜅1

𝜅2

𝜅3

(cid:104)

(cid:105)

=

−1.9683 2.3416 0.1115

(cid:105)

(2.12)

These coefficient values can be substituted in (2.10) to obtain the following power-law expression:

(cid:101)(cid:164)𝑊 𝑖𝑛 = 1.1180 [𝜏∗

𝑤,𝑚𝑎𝑥 ]2.3416 [ℎ∗]−1.9683

(2.13)

The power-law obtained in (2.13) collapses the data on a straight line in a log-log plane - this is

shown in Figure 2.8.

By fitting equations (2.9a), (2.9b) and (2.9c) to the CFD data, we get the following results for the

matrix of coefficients C = [𝐶𝑖 𝑗 ], 𝑖, 𝑗 = 1, 2, 3; the coefficients can then be computed indirectly as

follows:

C =












−1.4703

2.1924 −5.9854

0.8564 −0.7021

1.9268

0.0867

0.3477 −0.9874












(cid:104)

⇒

𝜅1

𝜅2

𝜅3

(cid:105)

(cid:104)

=

−1.9554 2.3568 0.1094

(cid:105)

(2.14)

A comparison of the results in (2.12) and (2.13), obtained through direct and indirect fit of the data,

indicates that 𝜅1 and 𝜅2 differ by 0.65%, the results for 𝜅3 differ by 1.9%. A goodness-of-fit analysis

of the direct fit was conducted via calculation of the relation coefficient and an Anderson-Darling

normality test (Stephens, Stephens (1974); D’Agostino et al., D’Agostino and Stephens (1986)).

The direct fit exhibits a relation coefficient of 𝑅2 = 0.9951, which suggests an adequate

goodness-of-fit. But the Anderson-Darling normality test (Stephens, Stephens (1974); D’Agostino

et al., D’Agostino and Stephens (1986)) is a stronger descriptor because it determines if the error

between the model and the data is random and normally distributed. If the model fails the test,

it indicates that the model is likely excluding higher-order terms. The null hypothesis, 𝐻0, states

that if the data comes from a normal distribution. To test 𝐻0, one needs the standardized residuals

of the experimental measurements 𝜀𝑖. These are the differences between the computed (using

16

the coefficients obtained with the direct fit) and experimental values of ln (cid:101)(cid:164)𝑊, normalized by their

unbiased standard deviation:

𝜀𝑖 =

(cid:17)

(cid:16)
(cid:101)(cid:164)𝑊 𝑖

ln

− (cid:2)𝜅1 ln (ℎ∗

𝑖 ) + 𝜅2 ln (𝜏∗

𝑖 ) + 𝜅3

(cid:3)

(cid:118)(cid:117)(cid:117)(cid:117)(cid:116)

𝑁
(cid:205)
𝑗=1

(cid:110)

(cid:17)

(cid:16)
(cid:101)(cid:164)𝑊 𝑗

ln

−

(cid:104)

𝜅1 ln (ℎ∗

𝑗 ) + 𝜅2 ln (𝜏∗

𝑗 ) + 𝜅3

(cid:105) (cid:111)2

(2.15)

where 𝑁 = 45 is the number of CFD data points. The following test statistic is used to assess 𝐻0:

𝑁 − 1

𝐴2 = −𝑁 −

𝑁
∑︁

𝑖=1

2𝑖 − 1
𝑁

(cid:8)ln 𝐹 (𝜀𝑖) + ln (cid:2)1 − 𝐹 (cid:0)𝜀 𝑁 +1−𝑖

(cid:1)(cid:3) (cid:9)

where 𝐹 (𝜀𝑖) is the Empirical Cumulative Distribution Function (ECDF) and is defined as:

𝐹 (𝜀𝑖) =

𝑝(𝜀𝑖)

∑︁

𝜀𝑖 ≤𝜖

(2.16)

(2.17)

Equation (2.17) will be used in equation (2.16) to calculate 𝐴2, which will provide a quantitative test

of 𝐻0. Considering a significance level of 𝛼 = 0.05, 𝐴2 has a critical value of 0.7385 (D’Agostino

et al., 1986), which means that 𝐻0 should be rejected if 𝐴2 > 0.7385 when calculated with the

experimental data. For the current data set, 𝐴2 = 0.3352 < 0.7385, showing that the standardized

error follows a normal distribution. Hence, the linear model in (2.10) is adequate to describe the

data. A graphical verification of this result is provided in Figure 2.7 where the ECDF is plotted

together with the Normal Cumulative Distribution Function (NCDF), which is defined as:

𝐹 (𝜀) =

1
√
2𝜋

𝜎

𝜀

∫

−∞

𝑒− 1

2

𝜂2 𝑑𝜂

(2.18)

In the study by Kamensky et al. Kamensky et al. (2019),

(cid:164)𝑊𝑖𝑛 was non-dimensionalized as

follows:

(cid:164)𝑊 ∗

𝑖𝑛 =

(cid:164)𝑊𝑖𝑛
𝑖𝑛𝐷 (𝐷/𝑑)

𝜇 𝑢2

(2.19)

To put the current results in context with that study, (cid:101)(cid:164)𝑊 𝑖𝑛 can be shown to be equal to (cid:164)𝑊 ∗

𝑖𝑛𝐷∗2.

Therefore, we get from (2.13)

(cid:164)𝑊 ∗

𝑖𝑛 = 1.1180 [𝜏∗

𝑤,𝑚𝑎𝑥 ]2.3416 [ℎ∗]−1.9683 [𝐷∗]−2

(2.20)

17

Figure 2.7 Plot showing the ECDF 𝜀𝑖 and the NCDF 𝜀

Based on theoretical grounds, Kamensky et al. Kamensky et al. (2019) proposed a proportionality

of the form:

(cid:164)𝑊 ∗

𝑖𝑛 ∝ [ℎ∗]−3, which was near to the scaling observed in that work. But equation

(2.20) suggests a proportionality closer to:

(cid:164)𝑊 ∗

𝑖𝑛 ∝ [ℎ∗]−2. These results are different because they

describe different phenomena. While the former is based on the global effect of inlet fluid power

on equilibrium height, the latter focuses on the local effects (where the shear stress is maximum)

that relate fluid power with maximum shear stress and equilibrium height.

2.5 Discussion

2.5.1 Relationship between maximum shear stress and inlet Reynolds number

The location of the maximum shear stress is just beyond the location of the corner, 𝑟 ≈ 𝑑/2;

this short development length indicates a laminar boundary layer over the range of flow speeds

examined here. The shear stress for a large number of laminar flows takes the form 𝜏 ∝ 𝜌𝑈2𝑅𝑒𝑛Λ,

where 𝑈 is the flow velocity and Λ is a flow-dependent length scale. The value 𝑛 = −1/2 has been

found for pure boundary layers with zero and non-zero pressure gradients (White White (1974)).

This value of 𝑛 was also found for the developed portion of the flow under a Bernoulli pad by (Guo

18

Figure 2.8 Plot showing CFD data collapsing on a straight line corresponding to the power-law in
(2.13)

et al. Guo et al. (2017)). The same is true for the maximum wall shear produced by both the planar

and axisymmetric wall jet, which is an unconfined analog of the developing region of the Bernoulli

pad flow (Phares et al. Phares et al. (2000)).

The latter work uses laminar boundary layer theory to yield (in our nomenclature) the expression

𝜏𝑤,𝑚𝑎𝑥 ∝ 𝜌𝑢2

𝑖𝑛𝑅𝑒𝑛

𝑖𝑛 [ℎ∗]−2 for the maximum shear stress. The use of the gap height arises in the

wall jet because of the development of the shear layer at the edge of the jet. Using our scaling

relationships, this expression predicts that

𝜏∗
𝑤,𝑚𝑎𝑥 =

𝜏𝑤,𝑚𝑎𝑥
𝜇𝑢𝑖𝑛𝑑/ℎ2
𝑒𝑞

∝

𝜌𝑢2
𝑖𝑛
𝜇𝑢𝑖𝑛𝑑/ℎ2
𝑒𝑞

𝑅𝑒𝑛

𝑖𝑛 [ℎ∗]−2 =

(cid:18) 𝜌𝑢𝑖𝑛𝑑
𝜇

(cid:33)

(cid:19) (cid:32) ℎ2
𝑒𝑞
𝑑2

𝑅𝑒𝑛

𝑖𝑛 [ℎ∗]−2 = 𝑅𝑒𝑛+1
𝑖𝑛

⇒ 𝜏∗

𝑤,𝑚𝑎𝑥 ∝ 𝑅𝑒𝑛+1
𝑖𝑛

that is, the non-dimensional shear stress is not dependent on any length scale other than the one

appearing in the Reynolds number. We have found this to be the case (see Figure 2.9), where the

(2.21)

19

value 𝑛 + 1 = 0.42, or 𝑛 = −0.58, is broadly in line with the findings for other flow fields. We can

write equation (2.13) as a proportionality in terms of 𝑅𝑒𝑖𝑛 and ℎ∗

(cid:101)(cid:164)𝑊 𝑖𝑛 ∝ (cid:2)𝑅𝑒0.42

𝑖𝑛

(cid:3) 2.3416 [ℎ∗]−1.9683 = [𝑅𝑒𝑖𝑛]0.98 [ℎ∗]−1.9683

(2.22)

Figure 2.9 Best fit of CFD data showing power law relationship between 𝜏∗

𝑤,𝑚𝑎𝑥 and 𝑅𝑒𝑖𝑛

The exponents which appear in Eq.(2.22) are enticingly near to round numbers but we have not been

able to predict them using analytical methods. This statement also holds for expressions presented

in the following section. We therefore present these relationships as implications of the data rather

than derivations.

2.5.2

Implications of the observed scaling relationship

We can write the power (cid:164)𝑊𝑖𝑛 required to operate the pad in terms of the shear stress

(cid:164)𝑊𝑖𝑛 ∝

∫ 2𝜋

∫ 𝐷/2

0

0

𝜏𝑤 (𝑟)𝑈 (𝑟) 𝑟𝑑𝑟𝑑𝜃 ∝

∫ 𝑟𝑤

0

𝜏𝑤,𝑚𝑎𝑥 𝑈 (𝑟) 𝑟𝑑𝑟

(2.23)

20

where 𝑟𝑤 is a radial location which characterizes an area over which we have to integrate the shear

stress, and 𝑈 (𝑟) is the mean radial component of the velocity in the gap. There are two critical

assumptions being made. The first assumption is that the power required to operate the system

is well-characterized by the integral of the wall shear over the entire pad area; this is consistent

with a low fraction of the inlet power being used to provide kinetic energy to the fluid (Kamensky

Kamensky et al. (2019)). The second assumption is that the integral of the wall shear over the

entire pad area can be characterized using the maximum shear stress, integrated over a much

smaller area. We justify this partially on the basis of the mean radial velocity 𝑈 (𝑟) ∝ 𝑟 −1. Because

𝜏𝑤 (𝑟) ∝ [𝑈 (𝑟)]2, 𝜏𝑤 (𝑟) 𝑈 (𝑟)𝑟 ∝ 𝑟 −2, and therefore most of the power used to overcome wall shear

is expended at small 𝑟. This assumption is further justified by the observed empirical form, in

which the outer diameter 𝐷 does not explicitly appear. We now manipulate the integrand, using

the result from mass conservation:

𝑈 (𝑟) ∝

𝑢𝑖𝑛𝑑
ℎ∗𝑟

⇒ 𝜏𝑤,𝑚𝑎𝑥 𝑈 (𝑟)𝑟 ∝ 𝜏∗

𝑤,𝑚𝑎𝑥

(cid:18) 𝜇𝑢𝑖𝑛
𝑑 [ℎ∗]2

(cid:19) (cid:18) 𝑢𝑖𝑛𝑑
ℎ∗𝑟

(cid:19)

𝑟 ∝ 𝜇𝑢2

𝑖𝑛𝜏∗

𝑤,𝑚𝑎𝑥 [ℎ∗]−3

Placing this form into equation (2.23) and making the result non-dimensional, yields

(cid:164)𝑊𝑖𝑛 ∝

∫ 𝑟𝑤

0

𝜇𝑢2

𝑖𝑛𝜏∗

𝑤,𝑚𝑎𝑥 [ℎ∗]−3𝑑𝑟 ⇒ (cid:101)(cid:164)𝑊 𝑖𝑛 ∝

1
𝜇𝑢2
𝑖𝑛𝑑

∫ 𝑟𝑤

0

𝜇𝑢2

𝑖𝑛𝜏∗

𝑤,𝑚𝑎𝑥 [ℎ∗]−3𝑑𝑟 = 𝜏∗

𝑤,𝑚𝑎𝑥 [ℎ∗]−3 (cid:16) 𝑟𝑤
𝑑

(cid:17)

(2.24)

The empirical relationship in equation (2.13) has the form [𝜏∗

𝑤,𝑚𝑎𝑥]𝑚 [ℎ∗]−2, where 𝑚 = 2.34. A

comparison between this form and the above relationship, and substitution of equation (2.21) along

with the value of the observed constant 𝑛 = −0.58, yields

𝑟𝑤
𝑑

∝ [𝜏∗

𝑤,𝑚𝑎𝑥]𝑚−1ℎ∗ ⇒

𝑟𝑤
ℎ𝑒𝑞

∝ 𝑅𝑒(𝑚−1)(𝑛+1)
𝑖𝑛

≈ 𝑅𝑒0.56

𝑖𝑛

(2.25)

This indicates that the length scale in the problem which best represents the power dissipated by

wall shear is the equilibrium gap height ℎ𝑒𝑞, but with a non-trivial dependence on the inlet Reynolds

number. The reason for this Reynolds number dependence is likely related to both the growth of

21

the boundary layer on the wall, and the scale of the recirculation bubble on the pad (Nakabayashi

Nakabayashi et al. (2002)).

As a final comment, we note that the “entrance length” of a pipe or channel is defined as a region

with wall shear substantially above the value in the fully developed flow (White White (1974)).

This is analogous to our definition of 𝑟𝑤 within the gap between the pad and the wall. For planar

channels, the dimensionless entrance length scales with 𝑅𝑒1.0. Equation (2.25) indicates that the

entrance length to a radial outflow scales instead with 𝑅𝑒0.56

𝑖𝑛

, notably reduced from the planar case.

2.5.3 Hydraulic roughness

The roughness used in our simulations follows the model of monodisperse sand grain (Niku-

radse Nikuradse (1950)). That source uses the roughness Reynolds number to define the wall as

hydraulically smooth for 𝜀+ ≲ 5, transitionally rough for 5 < 𝜀+ < 70, and rough for 𝜀+ ≳ 70,

where

𝜀+ ≜ 𝜀𝑢𝜏 𝜌/𝜇,

𝑢𝜏 ≜ √︁𝜏𝑤/𝜌

Although the present study was not undertaken in a zero-pressure gradient boundary layer, we

can use these values to provide context for the material roughness employed in the simulations by

defining a local value for 𝜀+ based on the local value of 𝑢𝜏.

The flow was hydraulically smooth predominantly over most of the of the pad’s surface in all

simulations. The worst case was observed for the nominal pad: 𝑑 = 31.75 mm, 𝐷 = 200 mm,

ℎ𝑒𝑞 = 1.5 mm, 𝜀 𝑝 = 𝜀𝑤 = 0.04 mm, which showed 𝜀+ ≈ 8.096 at the point of maximum shear

stress. For this worst case, there was a mild transitionally rough behavior around the region of

maximum shear.

A full examination of the roughness in the context of the proposed biofouling removal application

would probably require a more detailed description of the roughness, to include material properties

and inhomogeneities. For the present work, we have restricted ourselves to a real material for which

we have experimental data, and with a roughness which is small with respect to the gap height.

22

2.6 Concluding Remarks

A computational study of Bernoulli pads has been performed in order to determine a relationship

between the input power required to maintain stable equilibrium and the maximum shear stress.

The scaling law presented collapses the power and shear stress results of a large number of pad

geometries and operating conditions onto a single curve. This collapse is extremely robust over

the range of parameters examined. Because the inlet power is computed using a combination of

the pressure and velocity at the inlet, it is not obvious ita priori that the inlet power would be so

well-described by the maximum value of shear stress. Kamensky et al. Kamensky et al. (2019),

indicated that the fluid power scaled with “worst’ case” choices of characteristic parameters: largest

velocity, largest gradient, and largest area. Although the largest velocity and gradient occur near

the location of maximum shear stress, this position is much nearer to the jet exit 𝑑 than the outer

diameter 𝐷. These findings are not inconsistent, since the outer diameter is being used in Eq.(2.8)

to scale the power, equilibrium height, and shear stress. It appears that (per the arguments made in

the Discussion) the power expended to overcome wall shear over the entire pad is proportional to

the maximum value.

Unusually, this work has used water as the working fluid in its simulations of the Bernoulli

pad flow field. This is consistent with our proposed application, which is to use a Bernoulli pad

to groom underwater surfaces without contact. This choice of working fluid introduces a great

deal of potential richness - particularly cavitation - which has been ignored in the present work.

Because cavitation occurs when the local pressure is below the vapor pressure, bubble formation

tends to happen where the local velocity is large - precisely where we observe the highest wall

shear stress. The design of alternate flow geometries which create sufficiently high levels of shear

stress to groom surfaces while controlling cavitation will be examined in future work. Although it

is unlikely that the scaling laws proposed here will work ”off-the-shelf” on these new geometries,

we anticipate that the techniques used here will be employed.

23

CHAPTER 3

FLOW PHYSICS OF A ROTATING BERNOULLI PAD

3.1

Introduction

In the last chapter, the power-law relationship obtained for a Bernoulli pad suggests that there

is a direct relationship between the fluid power delivered to the system and the losses due to shear

in the radial direction. It is important to maximize shear forces in applications such as non-contact

biofouling mitigation Kamensky et al. (2020); Tomar et al. (2022). To this end, it is important to

consider other methods of delivering power to the system. As shown in Figure 3.1, we can envision

rotating the pad about its axis, whereby mechanical power is provided in addition to fluid power.

The additional power delivered to the system through rotation will increase losses due to shear,

but the shear forces will now act in both radial and tangential directions. It is interesting to note

that shear forces in the radial direction decrease with an increase in radial distance, whereas shear

forces in the tangential direction increase with the increase in radial distance.

In the absence of a prime mover, a rotating pad produces attractive forces for all gap heights and

does not produce a stable equilibrium configuration with the workpiece. The generation of attractive

forces (or subatmospheric pressure) through the action of rotation (mechanical power) was first

demonstrated experimentally by Demour in 1730 Houghtalen et al. (1996) using a T-shaped pipe;

here it is demonstrated using simulations in Section 3.3. Interestingly, the addition of fluid power

Figure 3.1 A schematic diagram of a Bernoulli pad provided with the additional capability of
rotation.

24

to the rotating pad changes the nature of equilibria from unstable to stable and results in a rotating

Bernoulli pad; this is presented in Section 3.4. A rotating Bernoulli pad has a unique property

that is not observed with the standard non-rotating Bernoulli pad, namely the ability to produce a

stable equilibrium configuration for two different mass flow rates; the attractive force is the same

for both mass flow rates. The ability to produce the same force with two different flow rates can

also be observed with a standard non-rotating Bernoulli pad but only in the repulsive regime; this

has not been discussed in the literature, potentially because of its limited range of applications.

Another important observation is that for the same level of inlet fluid power, a rotating Bernoulli pad

generates a higher level of maximum shear stress compared to the standard non-rotating Bernoulli

pad; this is discussed in Section 3.6.

The additional shear forces produced by a rotating Bernoulli pad increase the efficacy of

grooming surfaces subjected to biofouling. "Grooming” Tribou and Swain (2010) is the gentle and

frequent mechanical maintenance of a ship’s hull to keep it clear of fouling organisms. This is an

alternative to a less-frequent high-intensity “cleaning” procedure to avoid damaging coatings on

the hull. These coatings are the first line of defense against biofouling and fall largely into two

categories: antifouling (AF) and fouling-release (FR). AF coatings are biocidal, while FR coatings

are slippery to fouling organisms. This allows shear forces to clean the ship while in motion

Donnelly et al. (2022), but some level of preventative maintenance is also required Tribou and

Swain (2010); Yeginbayeva and Atlar (2018), depending on coating type, local fouling pressure,

time of year, and length of time in port. The in-water mechanical grooming tools are shear-based,

whether by direct contact like brushesDahlgren et al. (2022) or using fluid forces like water jets

Oliveira and Granhag (2020), cavitating jets Song and Cui (2020), ultrasonics Trickey et al. (2022),

bubbles Menesses et al. (2017), or a standard Bernoulli pad Kamensky et al. (2020); the standard

Bernoulli pad can be augmented by rotation, which we describe in the present work.

In this work, we use a computational approach to examine the features of rotating pads and

rotating Bernoulli pads. The computational approach is described and validated by our published

data Tomar et al. (2022) in Section 3.2. Rotating pads act as centrifugal pumps and produce only

25

unstable equilibria; in contrast, rotating Bernoulli pads produce stable and unstable equilibria. We

provide findings on the stiffness of the rotating systems, energy balance, and improvements to the

normal force and maximum shear stress resulting from the addition of rotation.

3.2 Approach

3.2.1 Computational Domain

A rotating Bernoulli pad is shown in Figure 3.1; 𝑑 and 𝐷 denote the inner and outer diameter,

ℎ denotes the gap between the pad and the workpiece, 𝜀p and 𝜀w denote the height of the roughness

of the pad and the workpiece surfaces, 𝜌 denotes the density of the fluid, 𝑢in represents the fluid

velocity at the inlet, 𝑝in represents the pressure at the inlet, and the angular speed of rotation of the
pad is denoted by 𝜔 = (cid:164)𝜃. The pressure at the outlet of the flow is atmospheric, i.e., 𝑝out = 𝑝atm.

The values of the parameters used for the Computational Fluid Dynamics (CFD) simulations

are provided in Table 3.1; these values match those used in our previous investigation Tomar et al.

(2022) with the exception of the value of 𝑑. A larger value of 𝑑 is used here to ensure that a stable

equilibrium can be achieved with a zero net normal force on the workpiece Kamensky et al. (2019)

for a reasonable mass flow rate. The working fluid is water, and the surface roughness values used

here are consistent with the material of the pad and workpiece used in the experiments Tomar et al.

(2022); Kamensky (2020) i.e., structural steel.

The investigation method used in this work is CFD simulations carried out using ANSYS

Fluent. A 2-D steady-state incompressible Navier-Stokes equation is solved using a coupled

scheme, which provides a robust and efficient steady-state solution compared to other segregated

Table 3.1 Parameter values used in CFD simulations

𝜔 [rad/s]
𝑑
𝐷
ℎ
𝜀 𝑝 = 𝜀𝑤
𝜌
𝜇

(cid:164)𝑚in [kg/s]
63.5 mm
200 mm
1.5 mm
0.04 mm
998.2 kg/m3
0.001 Ns/m2

26

Figure 3.2 Axisymmetric flow domain of the rotating Bernoulli pad with hybrid mesh.

schemes ANSYS (2022). Because the geometry and boundary conditions are axisymmetric, we

used a 2-D axisymmetric solver to reduce run time. The pad rotation 𝜔 about the 𝑧-axis is imposed

using an axisymmetric swirl module, which assumes zero gradient but non-zero velocity in the

circumferential direction 𝜃. At the inlet of the computational domain, we impose a mass flux

boundary condition (cid:164)𝑚𝑖𝑛. At the outlet, the gauge pressure is zero, i.e., 𝑝atm = 0. The workpiece is

kept stationary with a no-slip boundary condition at a gap height ℎ from the rotating pad.

The pressure is discretized using the PRESTO (PREssure STaggering Option) since it is well

suited for steep pressure gradients generated by the flow impinging on the workpiece Islam et al.

(2020). The Reynolds Stress Model (RSM) is used to model turbulence, which is widely used for

studies involving swirling motion Gibson et al. (2017) and vortices Neto et al. (2011). The RSM

is well suited to the severe streamline curvature caused by the swirling motion of the pad Jawarneh

and Vatistas (2006). The spatial terms for momentum, swirl velocity, and Reynolds stresses are

discretized using a second-order upwind scheme. The mesh of the computational domain is shown

in Figure 3.2; a magnified view of the mesh at the neck of the Bernoulli pad is shown inset. The

mesh is hybrid and is comprised of quadrilateral (dominant) and triangular elements.

3.2.2 Experimental Validation

The numerical methods are compared with the results from a Particle Tracking Velocimetry

(PTV) experiment Kamensky (2020) to gain confidence in the numerical solver and discretization

27

Figure 3.3 Comparison of PTV (data taken from Tomar et al. (2022)) and CFD data: plot of the
nondimensional radial velocity at 𝑧 = ℎ/2 with respect to the nondimensional radial position.

schemes. The experimental work Kamensky (2020) used a different set of parameters (𝑑 = 15.875

mm, 𝐷 = 100 mm, ℎ = 1.275 mm); these values were used in the numerical simulation for

comparison. The pad was not rotating (𝜔 = 0), and the value of the mass flow rate was (cid:164)𝑚in = 0.59

kg/s.

If 𝑣 denotes the radial velocity at the middle of the gap (𝑧 = ℎ/2), the non-dimensional

radial velocity is defined as 𝑣∗ = 𝑣/𝑣max, where 𝑣max is the maximum value of the radial velocity.

Similarly, the non-dimensional radial position is defined as 𝑟 ∗ = 2𝑟/𝐷. The non-dimensional

radial velocity for the numerical solution and the experiment are plotted as a function of the non-

dimensional radial position in Figure 3.3. The velocity profile shows good agreement between

the numerical simulation and PTV experiments. In particular, both cases show 𝑣max at nearly the

same location (𝑟 ≈ 16.62 mm), very close to the neck, where the fluid enters the gap in the radial

direction. This indicates that the numerical model captures the recirculation just downstream of

the neck in a satisfactory manner.

28

Figure 3.4 Variation of y+ values over the workpiece in the radial direction.

3.2.3 Grid Independence Study

The CFD simulations are carried out using an element size of 4 × 10−5 m, which resulted in

94, 547 elements and a maximum 𝑦+ value of 0.27 over the entire range of the workpiece as shown

in Figure 3.4. A smaller element size of 3×10−5 m resulted in 127, 189 elements, but the percentage

change in the values of the variables of interest, namely, the normal force and the maximum shear

stress on the workpiece, were found to be 1.09 % and 0.54 %, respectively, and did not justify the

additional computational cost.

3.2.4 Comparison of axisymmetric swirl module with three-dimensional domain

We have carried out a numerical simulation of a three-dimensional computational domain with

identical dimensions used in Table3.1. The domain is provided with a mass flow rate of 0.5 kg/s

and a rotational speed of 10 rad/s. The results obtained with the three-dimensional simulations

are used to compare the results obtained with the axisymmetric swirl module. A comparison of

wall shear stress generated on the workpiece and radial velocity at the center of the gap height is

shown in Figure 3.5 and Figure 3.6. The qualitative nature of the variation of wall shear stress and

29

Figure 3.5 Variation of wall shear stress over the workpiece in the radial direction for 3D and
axisymmetric swirl module.

radial velocity is captured well in both cases. The error between the maximum value of wall shear

stress and radial velocity is also small, which provides us with confidence in the axisymmetric swirl

module of the solver. The difference in the magnitude of the maximum velocity can be attributed

to the coarser mesh used in the axisymmetric domain and it can be avoided by further refining the

mesh but it is out of the scope for this kind of work.

3.3 Flow physics of a rotating pad

A Bernoulli pad involves a prime mover that pumps fluid through the inlet - see Figure 3.1.

To understand the effect of rotation of the Bernoulli pad on the flow physics, we first consider a

rotating pad with the inlet exposed to atmospheric pressure, i.e., 𝑝in = 𝑝atm. Numerical simulations

are carried out to simulate the flow physics of a rotating pad for the parameter values in Table 3.1

and the set of angular speeds 𝜔 ∈ {10, 20, 30, 40, 50, 60} rad/s. The rotating pad can be considered

to be a centrifugal pump.

For the rotating pad, numerical computation of the flow is challenging since both the inlet and

30

Figure 3.6 Variation of radial velocity at ℎ = ℎ/2 over the workpiece in the radial direction for 3D
and axisymmetric swirl module.

Table 3.2 (cid:164)𝑚in required to achieve 𝑝in = 𝑝atm for different 𝜔

𝜔 [rad/s]
10
20
30
40
50
60

(cid:164)𝑚in [kg/s]
0.0520
0.1100
0.1700
0.2280
0.2830
0.3389

the outlet have the same pressure boundary conditions, namely, 𝑝in = 𝑝out = 𝑝atm. To address this

challenge, the pressure boundary condition is imposed at the outlet, and the inlet mass flow rate

(cid:164)𝑚in is varied iteratively till the condition 𝑝in = 𝑝atm is satisfied; this condition corresponds to zero

static pressure at the inlet. The (cid:164)𝑚in required to satisfy the pressure boundary conditions for the

chosen set of angular velocities is provided in Table 3.2.

The pressure 𝑝 on the workpiece surface is presented in Figure 3.7 for the set 𝜔 ∈ {10, 20, 30, 40, 50, 60}

31

rad/s. It is equal to the stagnation pressure in the region below the inlet; it is small but different

for different rotational speeds. The difference in the stagnation pressure can be attributed to the

difference in the mass flow rates at the inlet (cid:164)𝑚in - see Table 3.2. A higher speed of rotation results

in a higher value of (cid:164)𝑚in and a higher value of stagnation pressure at the workpiece. The pressure

on the workpiece at the neck, where the axial flow becomes radial, drops rapidly due to the large

increase in the radial velocity owing to the constriction of the flow - see Figure 3.8.

For a given pad and rotational speed, the gap between the pad and the workpiece dictates the

normal force generated on the workpiece. The normal force 𝐹, which is obtained by integrating the

pressure on the workpiece (see Figure 3.7), is negative and hence attractive (attempts to pull the pad

and the workpiece closer). An external force, equal and opposite to the normal force, can maintain

equilibrium, but the nature of this equilibrium (stable or unstable) will depend on the stiffness of

the flow field in the 𝑧 direction - see Figure 3.1. The stiffness can be computed from the change in

Figure 3.7 Rotating pad with parameters in Table 3.1: Pressure on the workpiece surface as a
function of the radial position for different angular speeds.

32

Figure 3.8 Contours of (a) pressure and (b) nondimensional radial velocity (𝑣/𝑢in) for a rotating
pad. A portion of the stagnation region 𝑟 ∈ [0, 8] mm and the neck region of the pad 𝑟 ∈ [24, 40]
mm are provided along with a magnified view of the neck region.

the normal force 𝐹 due to a change in the gap, as follows:

𝑘 = −

Δ𝐹
Δℎ

(3.1)

where Δ𝐹 is the change in the normal force and Δℎ is the change in the gap. The stiffness of the

flow field is shown in Figure 3.9 for different angular speeds about the nominal gap of ℎ = 1.5 mm.

The negative value of the stiffness indicates that a bias force can be used to generate an equilibrium,

but the nature of the equilibrium is unstable for all 𝜔.

3.4 Flow physics of a rotating Bernoulli Pad

We now consider a rotating Bernoulli pad, where a prime mover does non-negative work to

pump the fluid through the inlet, and the pad has an angular speed 𝜔. Similar to the rotating pad,

the inlet boundary condition is the mass flux (cid:164)𝑚in, but the mass flux is not iterated to find an inlet

33

Figure 3.9 Stiffness of the flow field generated by the rotating pad at ℎ = 1.5 mm for
𝜔 ∈ {10, 20, 30, 40, 50, 60} rad/s.

pressure of zero; i.e., 𝑝in ≠ 𝑝atm. The force on the workpiece produced by varying the mass flow

rate is provided in Figure 3.10 for 𝜔 ∈ {0, 20, 40, 60}. It can be seen that the lowest value of (cid:164)𝑚in is

different for each curve corresponding to 𝜔; these different (cid:164)𝑚in values correspond to the boundary

condition 𝑝in = 𝑝atm, provided in Table 3.2. The motivation to use these (cid:164)𝑚in values as the starting

point is to ensure that the work done by the prime mover is not negative.

In Figure 3.10, the 𝜔 = 0 curve corresponds to the standard non-rotating Bernoulli pad. The

variation in the normal force 𝐹 with increasing mass flow rate (cid:164)𝑚in is in agreement with results in

the literature Li and Kagawa (2014): as (cid:164)𝑚in is increased, 𝐹 initially increases and becomes positive

(repulsive), reaches a maximum value, then decreases and becomes negative (attractive). When

the pad is provided an angular speed, the concave-downward trend in the variation of 𝐹 with (cid:164)𝑚in

remains unchanged, but the values of the normal force uniformly decrease. The decrease in the

normal force is proportional to the increase in the angular speed; this is in agreement with our

discussion in Section 3.3.

The rotating Bernoulli pad can be viewed as a modification of the standard non-rotating Bernoulli

34

Figure 3.10 Bernoulli pad with parameters in Table 3.1: Normal force on the workpiece as a
function of the mass flow rate for (cid:164)𝑚in

Figure 3.11 Stiffness of the flow field generated by the rotating Bernoulli pad at ℎ = 1.5 mm and
(cid:164)𝑚in = 0.5 kg/s, for 𝜔 ∈ {10, 20, 30, 40, 50, 60} rad/s.

35

pad where the radial velocity in the gap is primarily due to the fluid power provided at the inlet,

and the tangential velocity in the gap is due to pad rotation. Therefore, as (cid:164)𝑚in is increased, the

effect of the radial velocity becomes prominent compared to that of the tangential velocity, and

the curves for nonzero 𝜔 in Figure 3.10 sequentially converge to the curve corresponding to the

standard non-rotating Bernoulli pad with 𝜔 = 0. The convergence of the curves corresponding to

𝜔 = 0 and 𝜔 = 20 is observed for a lower mass flow rate compared to that for convergence of the

curves corresponding to 𝜔 = 0, 𝜔 = 20 and 𝜔 = 40. All the curves converge together for the mass

flow rate of (cid:164)𝑚in ≈ 1.6 kg/s.

Figure 3.11 shows the stiffness of the flow field generated by the rotating Bernoulli pad as a

function of 𝜔 at ℎ = 1.5 mm and for (cid:164)𝑚in = 0.5 kg/s. Unlike the rotating pad, the stiffness is positive

for all 𝜔, which indicates that an equilibrium configuration of the pad will be stable. This may be

reconciled with the flow physics of the rotating pad discussed in Section 3.3. As the mass flow rate

increases due to work done by a prime mover, the pressure at the inlet increases, and this provides

a repulsive normal force on the workpiece for small values of 𝑟.

The repulsive force balances the attractive force generated by the flow in the gap, which now

includes a significant radial component compared to the rotating pad and yields a stable equilibrium.

Thus, the mass flow rate provided by a prime mover converts a rotating pad into a rotating Bernoulli

pad and an unstable equilibrium configuration into a stable equilibrium configuration. Although the

results in Figure 3.11 pertain to (cid:164)𝑚in = 0.5 kg/s, all points on the four curves in Figure 3.10, except

the starting point on each curve which do not require a prime mover, results in stable equilibrium

configurations.

Another important observation that can be made from Figure 3.10 is that for a specific angular

velocity, a horizontal line will provide the same force values at two different mass flow rates. For

example, mass flow rates of (cid:164)𝑚in = 0.8 kg/s and (cid:164)𝑚in = 1.305 kg/s generate the same normal force of

𝐹 = −3.027 N when the Bernoulli pad is rotating with 𝜔 = 60 rad/s. Therefore, with a bias force

of +3.027 N, it is possible to have an equilibrium configuration with two different mass flow rates.

The variation of the pressure on the workpiece for these two mass flow rates is presented in Figure

36

3.12. The integral of the area for both these curves represents the normal force 𝐹 generated on the

workpiece and turns out to be equal in this case.

The stiffness values corresponding to these mass flow rates are, however, different, equal to

𝑘 = 16, 412 N/m for (cid:164)𝑚in = 0.8 kg/s and 𝑘 = 31, 982 N/m for (cid:164)𝑚in = 1.305 kg/s. The positive values

of the stiffness indicate that the equilibrium configuration will be stable, but the flow field will

have a higher stiffness for a higher mass flow rate, and vice versa. This unique phenomenon occurs

due to the concave-downward nature of the variation of 𝐹 with (cid:164)𝑚in (see Figure 3.10); therefore,

the phenomenon is observed for all angular speeds, including zero angular speed or the standard

non-rotating Bernoulli pad. Although a rotating Bernoulli pad has not been investigated in the

literature, the phenomenon of a stable equilibrium where the flow field can have two different

stiffnesses corresponding to two different mass flow rates has not been discussed in the context of

the standard non-rotating Bernoulli pad. Such equilibria will henceforth be referred to as bi-stiffness

stable equilibria.

Figure 3.12 Rotating Bernoulli pad with parameters in Table 3.1 and 𝜔 = 60 rad/s: Pressure on
the workpiece as a function of radial position for (cid:164)𝑚in = 0.8 kg/s and (cid:164)𝑚in = 1.305 kg/s.

37

Figure 3.13 A workpiece of weight 𝑊 being supported by a Bernoulli pad from (a) below and (b)
above.

A horizontal line that intersects the 𝜔 = 0 curve at two points (corresponding to two mass flow

rates) in Figure 3.10 results in a positive value of 𝐹, whereas a horizontal line that intersects the

𝜔 = 60 curve results in a negative value of 𝐹. This implies that a standard non-rotating Bernoulli

pad can support a workpiece in a bi-stiffness stable equilibrium configuration only from below - see

Figure 3.13 (a). In contrast, the Bernoulli pad rotating with 𝜔 = 60 rad/s can support a workpiece

in a bi-stiffness stable equilibrium configuration only from above - see Figure 3.13 (b). In industry

applications, a standard non-rotating Bernoulli pad is typically used to support a workpiece from

above ( (cid:164)𝑚in is chosen such that 𝐹 is negative in Figure 3.10; this produces a stable equilibrium

configuration); this may explain why bi-stiffness stable equilibria have not been discussed in the

literature.

The concept of bi-stiffness equilibria is not limited to a specific pair of mass flow rate curves.

In fact, the production of bi-stiffness equilibria is observed for any arbitrary choice of mass flow

rate values in the CFD simulations. The effect of rotational speed on the normal suction force is

presented in Figure 3.14 for a range of mass flow rate values. The value of normal suction force

increases rapidly with the rotational speed for lower mass flow rate values. Hence, every single

curve in Figure 3.14 will intersect with another curve at a unique location. This unique location

represents a pad with fixed gap height at which the system has the same normal suction force for a

given rotational speed.

A magnified view of the Figure 3.14 is shown in Figure 3.15. The legends in this figure for

38

Figure 3.14 Variation of force with angular speed for different (cid:164)𝑚in

different mass flow rates are similar to those in Figure 3.14. In Figure 3.15, data on the y-axis is

clipped to the positive value of the forces. The (cid:164)𝑚in = 0.5 kg/s curve intersect with the (cid:164)𝑚in = 0.75

kg/s curve at 𝜔 ≈ 15 rad/s and (cid:164)𝑚in = 1 kg/s at 𝜔 ≈ 32 rad/s. The (cid:164)𝑚in = 0.5 kg/s curve will not

intersect the (cid:164)𝑚in = 0.25 kg/s curve because the effect of rotation is higher on the smaller mass flow

rate curves and, therefore, the suction force generation is rapid.

3.5 Rotating Bernoulli Pad: Variation of Normal Force with Gap Height

The normal force generated by a standard non-rotating Bernoulli pad strongly depends upon the

gap height. The variation of the force with the gap height is well explained in the literature, Li and

Kagawa (2014). The sudden contraction of the fluid at the neck results in a pressure less than the

atmospheric pressure, thereby producing attractive forces. Here, we present CFD simulation results

to compare the pressure profile on the workpiece for rotating and non-rotating Bernoulli pads with

parameters in Table 3.1; the results are shown in Figure 3.16 for 𝜔 = 0 rad/s and 𝜔 = 60 rad/s.

For both cases, the inlet boundary condition was (cid:164)𝑚in = 0.5 kg/s. The normal force 𝐹, obtained by

39

Figure 3.15 A magnified view of the variation of force with angular speed for different (cid:164)𝑚in with
the legends similar to Figure 3.14

integrating the pressure curves in Figure 3.16, is equal to 1.4037 N for 𝜔 = 0 rad/s and −8.1153 N

for 𝜔 = 60 rad/s. This is expected as the pressure curve for 𝜔 = 60 rad/s uniformly lies below the

pressure curve for 𝜔 = 0 rad/s. The decrease in the value of 𝐹 from a positive to a negative number

(or from a repulsive to attractive force) is due to the effect of rotation. These forces are shown in

Figure 3.17, which will be discussed next.

The effect of rotation on the variation of the normal force on the workpiece with gap height is

shown in Figure 3.17 for a mass flow rate of (cid:164)𝑚in = 0.5 kg/s. Every design point in the figure, with

the addition of a suitable bias force, corresponds to a stable equilibrium configuration; this can be

verified from the uniformly negative slope of the curves corresponding to each angular speed 𝜔. It

can be seen that an increase in the angular speed uniformly lowers the value of the normal force

𝐹. This observation, made earlier for ℎ = 1.5 mm, is found to be true for all gap heights, implying

that a rotating Bernoulli pad can lift or support more weight than a non-rotating Bernoulli pad for

a given mass flow rate.

40

Figure 3.16 Rotating Bernoulli pad with parameters in Table 3.1 and (cid:164)𝑚in = 0.5 kg/s: Pressure on
the workpiece as a function of radial position for 𝜔 = 0 rad/s and 𝜔 = 60 rad/s.

A similar trend is observed for smaller values of mass flow rate and rotational rate. The effect of

lower rotational speed for a mass flow rate of (cid:164)𝑚in = 0.1 kg/s and (cid:164)𝑚in = 1.045 kg/s is shown in Figure

3.18 and Figure 3.19 respectively. All the design points on these plots correspond to the stable

equilibrium; this can be verified from the uniformly negative slope of the curves corresponding to

each angular speed 𝜔 ∈ {0, 5, 10, 15} rad/s. For a particular value of gap height, the addition of

rotation speed results in an increase of normal suction force (or a decrease in the normal repulsive

force). The effect of rotation on the Bernoulli pad is negligible for the design points with higher

values of mass flow rate, as shown in Figure 3.19. In this case, the suction flow physics of the

Bernoulli pad is dominated by the power available at the inlet.

The 𝜔 = 10 rad/s curve from the Figure 3.18 and Figure 3.19 is taken out and plotted against

each other in Figure 3.20. The curves corresponding to two different mass flow rates intersect each

other at ℎ = 1.5 mm, and the normal suction force corresponding to this orientation is zero. In fact,

these mass flow rate values were initially chosen to make sure that they correspond to the stable

41

equilibrium configuration with zero force. The figure also shows the difference in the stiffness value

for the same configuration at stable equilibrium. For the smaller value of stiffness, the stiffness

value is also small, which in turn means that a very small value of external force is required to pull

this configuration out of the stable equilibrium.

On the other hand, the external force required to disturb the equilibrium for (cid:164)𝑚in = 1.045 kg/s

is higher. This bi-stiffness equilibria is not limited to one rotational speed. Numerical simulations

are extended to obtain the zero normal force stable equilibrium configuration for 𝜔 = 5 rad/s and

𝜔 = 20 rad/s. For 𝜔 = 5 rad/s, the mass flow rate values (cid:164)𝑚in = 0.037 kg/s and (cid:164)𝑚in = 1.05 kg/s

results in zero normal force stable equilibrium configuration as shown in the Figure 3.21 and for

𝜔 = 20 rad/s, the mass flow rate values (cid:164)𝑚in = 0.257 kg/s and (cid:164)𝑚in = 1.036 kg/s results in zero

normal force stable equilibrium configuration as shown in the Figure 3.22. A similar behavior with

bi-stiffness equilibria is observed in both figures.

The variation of total pressure for the mass flow rate values plotted in the Figure 3.20 are shown

Figure 3.17 Bernoulli pad with parameters in Table 3.1, except ℎ: Normal force on the workpiece
as a function of the gap height for (cid:164)𝑚in = 0.5 kg/s and 𝜔 ∈ {0, 20, 40, 60} rad/s.

42

Figure 3.18 Normal force on the workpiece as a function of the gap height for (cid:164)𝑚in = 0.1 kg/s.

Figure 3.19 Normal force on the workpiece as a function of the gap height for (cid:164)𝑚in = 1.045 kg/s.

43

Figure 3.20 Force curve to obtain bi-stiffness eqilibria for 𝜔 = 10 rad/s.

Figure 3.21 Force curve to obtain bi-stiffness eqilibria for 𝜔 = 5 rad/s.

44

Figure 3.22 Force curve to obtain bi-stiffness eqilibria for 𝜔 = 20 rad/s.

in the Figure 3.23, Figure 3.21 in Figure 3.24 and the variation of total pressure for the mass flow

rate values plotted in the Figure 3.22 are shown in the Figure 3.25. The integral of the area for both

the curves in Figure 3.23, Figure 3.24 and Figure 3.25 represents the normal force 𝐹 generated on

the workpiece and turns out to be zero for all the cases. The total pressure scales for the plots are

drastically different because of the difference in mass flow rate values, and therefore, the plots have

different 𝑦 axis ranges for the different mass flow rate curves.

3.6 Rotating Bernoulli Pad: Energy Balance and Advantages

3.6.1 Energy Balance

Consider the fluid control volume shown by the dotted lines in Figure 3.1. For this control

volume, the energy balance equation will be comprised of the following terms:

45

Figure 3.23 Variation of total pressure for (cid:164)𝑚in in Figure 3.20.

Figure 3.24 Variation of total pressure for (cid:164)𝑚in in Figure 3.21.

46

Figure 3.25 Variation of total pressure for (cid:164)𝑚in in Figure 3.22.

1. Inlet Fluid Power: The power associated with the flow at the inlet is:

(cid:164)𝑊in =

𝜋𝑑2

4

𝑝total 𝑢in

(3.2)

where 𝑝total is the total pressure, comprised of the static and dynamic pressures. For the

rotating pad and the rotating Bernoulli pad, 𝑝total is given by the relations:

𝑝total =





1
2

𝜌𝑢2
in

: rotating pad

1
2

𝜌𝑢2
in

+ 𝑝in

: rotating Bernoulli pad

2. Outlet Fluid Power: The power associated with the flow at the outlet is:

(cid:164)𝑊out =

𝜋𝑑2

8

𝜌𝑢3

out

,

𝑢out ≜

𝑑2
8𝐷ℎ

𝑢in

where 𝑢out is the velocity of the fluid at the outlet.

(3.3)

(3.4)

47

3. Mechanical Power: The mechanical power required to rotate the pad is given by the expres-

sion:

(cid:164)𝑊mech = 2𝜋𝜔

∫ 𝐷/2

𝑑/2

𝜇

𝜕𝑢𝜃
𝜕𝑧

𝑟 2𝑑𝑟

4. Losses: The losses will be given by the expression:

(cid:164)𝑊loss =





(cid:164)𝑊mech − (cid:164)𝑊in − (cid:164)𝑊out : rotating pad

(cid:164)𝑊mech + (cid:164)𝑊in − (cid:164)𝑊out

: rotating Bernoulli pad

(3.5)

(3.6)

and can be explained as follows. For a rotating pad, the inlet fluid power (cid:164)𝑊in is the result
of the mechanical power expended in rotating the pad (cid:164)𝑊mech. The losses are, therefore,

computed by subtracting the inlet and outlet fluid power from the mechanical power added to

the system. In contrast to the rotating pad, the rotating Bernoulli pad uses a prime mover to

provide the inlet fluid power, and hence, the losses for the system are computed by subtracting

the outlet fluid power from the sum of the mechanical power and inlet fluid power.

The losses are a measure of the shear stresses on the pad and the workpiece. In the present work,

we focus on the shear stress on the workpiece - this is motivated by the application of non-contact

shear-based grooming of biofouled surfaces Kamensky et al. (2020); Tomar et al. (2022).

3.6.1.1 Work analysis for zero static pressure at the inlet

The flow physics involving the rotating pad without an inlet pump (Static pressure = 0) is

explained in detail in section 3 of this chapter.

In the Figure 3.26, the corresponding power

available at the outlet and the mechanical power required to rotate the pad are presented. The mass

flow rate at the inlet getting sucked in increases with an increase in the angular speed of the pad.

Since the pressure available at the inlet consists of the dynamic pressure, the power available at the

inlet is very small in magnitude.

In the present control volume, the energy given to the system is in the form of mechanical

energy to rotate the Bernoulli pad. The pad acts as a centrifugal pump and sucks the water from the

inlet. Therefore, the work available at the inlet (dynamic pressure) is a byproduct of the mechanical

48

Figure 3.26 Outlet fluid power and mechanical power required for rotation

work. The losses in the system account for the rest of the energy provided to the system. These

losses can be expressed as:

𝑊𝑙 = 𝑊𝑚𝑒𝑐ℎ − (𝑊𝑖 + 𝑊𝑜)

3.6.1.2 Work analysis for mass flow rate = 0.5 kg/s

The addition of the inlet pump provides the static pressure head at the inlet. The combined

effect of rotation and inlet fluid power provides a stable equilibrium to the configuration. Numerical

simulations are carried out for mass flow rate, m = 0.5 kg/s, and the work available at the inlet,

mechanical work for the rotation, and the work available at the outlet are presented in the Figure

3.27.

In addition to the mechanical work for rotation, inlet fluid power is also added to the system.

The contribution of mechanical power to the inlet fluid is negligible (refer to Table 3.2). Therefore,

the losses in the system in this configuration are estimated as:

𝑊𝑙 = (𝑊𝑚𝑒𝑐ℎ + 𝑊𝑖) − 𝑊𝑜

49

Figure 3.27 Inlet fluid power, outlet fluid power, and mechanical power required for rotation

𝜔 (rad/s)
10
20
30
40
50
60

(cid:164)𝑚 (kg/s) 𝑊𝑖𝑛 (W) 𝑊𝑜 (W) 𝑊𝑚𝑒𝑐ℎ (W) Losses (W)

1
1
1
1
1
1

6.1394
6.1237
6.0882
6.0335
5.9558
5.8515

0.7101
0.7500
0.8400
0.9968
1.2343
1.5701

0.1571
0.6780
1.6800
3.2840
5.6000
8.7456

5.5864
6.0517
6.9282
8.3207
10.3215
13.027

Table 3.3 Work analysis for a mass flow rate of 1 kg/s

3.6.1.3 Work analysis for mass flow rate = 1 kg/s

Numerical simulations are carried out for mass flow rate, 𝑚 = 1 kg/s, and the work available

at the inlet, mechanical work for the rotation, and the work available at the outlet are tabulated in

Table 3.3. The inlet fluid power added to the system reduces with an increase in rotation. The

contribution of mechanical power to the inlet fluid is again neglected. The losses in the system in

this configuration are estimated as:

𝑊𝑙 = (𝑊𝑚𝑒𝑐ℎ + 𝑊𝑖) − 𝑊𝑜

50

(cid:164)𝑚 (kg/s) 𝑊𝑖𝑛 (W) 𝑊𝑜 (W) 𝑊𝑚𝑒𝑐ℎ (W)
0.5
1

0.0956
0.7003

0.8304
6.1447

0
0

Table 3.4 Work analysis without rotation

3.6.1.4 Work analysis without rotation

To understand the computational domain without any rotation, numerical simulations are carried

out with inlet fluid power alone. The work available at the inlet and the outlet is tabulated in the

Table 3.4. This data from the table also gives a perspective about the magnitude of inlet fluid power

added to the system when compared to the cases with a rotating pad.

3.6.2

Improvement in Normal Force and Maximum Shear Stress for a Given Inlet Fluid
Power

The normal force 𝐹 and the maximum shear stress 𝜏max are important for the application of

biofouling mitigation through non-contact shear-based grooming. The normal force (attractive)

will be used to adhere to the surface, and the maximum shear stress will define the grooming

efficacy. The normal force 𝐹 has been discussed at length in Sections 3.4 and 3.5; therefore, we

present some discussion on the shear stress 𝜏 on the workpiece before presenting the main results.

For the same inlet fluid power (cid:164)𝑊in, the radial variation of the shear stress on the workpiece is shown

in Figure 3.28 for a Bernoulli pad with 𝜔 = 0 and 𝜔 = 60 rad/s.

The maximum shear stress on the workpiece 𝜏max is observed at a location close to the entrance

of the gap (neck) of the Bernoulli pad - see the figure inset, which shows a magnified view of the

region where the flow turns from the axial direction to the radial direction. The axial flow from the

inlet strikes the workpiece and passes radially through the gap between the pad and the workpiece.

At the neck, the sharp change in the flow direction does not allow the fluid to immediately attach

to the pad surface; this produces a recirculating flow in the separation and reattachment regions

Shi and Li (2016). This, in turn, reduces the effective area of the flow between the pad and the

workpiece.

A sudden rise in velocity in this region is observed, which results in higher velocity gradients on

the workpiece. Since the recirculating flow is located close to the neck, the location of the maximum

51

Figure 3.28 Radial variation of the shear stress 𝜏 for the Bernoulli pad with parameters in Table
3.1, for 𝜔 = 0 and 𝜔 = 60 rad/s with the same (cid:164)𝑊in: the maximum value 𝜏max occurs below the
recirculating flow, at the boundary of the separation and reattachment regions - see the figure
inset, taken fromShi and Li (2016).

shear stress is independent of the speed of rotation. The velocity gradients below the recirculating

flow are also proportional to the stagnation pressure associated with the fluid entering the gap. The

addition of rotation affects the stagnation pressure and, therefore, changes the magnitude of shear

stress on the workpiece.

We present simulation results to illustrate the effect of rotation on the normal force 𝐹 and

maximum shear stress 𝜏max on the workpiece for the same level of inlet fluid power, (cid:164)𝑊in ≈ 0.83 W

- see Figure 3.29. For a given angular speed, the mass flow rate at the inlet is changed to make sure

that the inlet fluid power remains the same: the mass flow rate is equal to (cid:164)𝑚in = 0.5 kg/s at 𝜔 = 0

rad/s and (cid:164)𝑚in = 0.555 kg/s at 𝜔 = 60 rad/s. The increase in rotational speed from 𝜔 = 0 rad/s to

𝜔 = 60 rad/s results in a 14.93% rise in the value of the maximum shear stress; the normal force

changes from +1.4037 N (repulsive) to −6.6827 N (attractive). The increase in the maximum shear

stress value will result in better grooming of the surface, and an increase in the attractive normal

52

force will enable the Bernoulli pad to adhere strongly to the workpiece.

3.7 Conclusion

A standard non-rotating Bernoulli pad uses a prime mover, such as a blower or pump, to produce

its radial outflow. The axial jet provided by the prime mover at 𝑟 ≤ 𝑑/2 is turned by the workpiece

to move the flow in the radial direction. For a given mass flow rate, the Bernoulli pad produces a

stable equilibrium configuration within a certain range of gap heights. In this configuration, which

is characterized by positive stiffness of the flow field, the repulsive force between the workpiece and

the pad due to the high pressure of the axial jet and the attractive force due to the subatmospheric

pressure in the region 𝑑/2 ≤ 𝑟 ≤ 𝐷/2, balance each other. The equilibrium configuration is

associated with a specific mass flow rate, and a change in the mass flow rate causes the equilibrium

gap height to change.

A rotating pad, in contrast, does not use a separate prime mover. The shear imposed by the

Figure 3.29 Variation in the normal force and maximum shear stress on the workpiece with
angular speed for a fixed inlet fluid power (cid:164)𝑊in ≈ 0.83 W; this power results in a mass flow rate of
(cid:164)𝑚in = 0.5 kg/s at 𝜔 = 0 rad/s.

53

rotating pad drives a tangential flow within the gap, which is accelerated in the radial direction.

The mass flux produced by this radial flow enters the flow field through the inlet; the rotating

pad is its own prime mover, essentially a bladeless centrifugal pump. Because the inlet velocity

is small compared to the tangential velocity, the pressure of the axial jet is small, and this flow

field produces a strong, attractive force between the pad and the workpiece. Unlike the standard

non-rotating Bernoulli pad, the flow field produced by the rotating pad has negative stiffness and

hence no stable equilibrium configuration, meaning that the attractive force between the pad and

the workpiece becomes stronger as they get closer.

A rotating Bernoulli pad uses a combination of fluid power provided by a prime-mover and

mechanical power for pad rotation. Similar to the standard non-rotating Bernoulli pad, this system

also exhibits stable equilibrium configurations. A unique characteristic of this system is that a

given stable equilibrium configuration is associated with two different stiffnesses corresponding to

two different mass flow rates for a given speed of rotation. The standard non-rotating Bernoulli

pad is a special case of the rotating Bernoulli pad and also exhibits this property, not reported in

the literature heretofore.

The rotating flow in the gap adds additional normal force between the pad and workpiece and

drives a small amount of additional mass flux at the same inlet pressure. This allows the balance

between the impinging pressure of the axial jet and low pressure within the gap to balance at an

additional value of mass flow rate. The pad rotation also increases the shear experienced by the

workpiece at constant inlet power, both by increasing the mass flux and by providing a tangential

velocity component. This additional shear, which comes at the cost of additional mechanical power,

can be used to improve grooming efficacy in biofouling mitigation.

For a specific pad with the same inlet fluid power, an increase in the rotational speed by 60 rad/s

resulted in the maximum wall shear stress increasing by ≈ 15% and the normal force changing

from +1.4 N (repulsive) to −6.6 N (attractive). This work is intended as an overview of some of

the many rich features of this flow field; we expect that many of these results can be extended and

generalized in future work.

54

CHAPTER 4

WALL SHEAR MEASUREMENT: HOT-FILM ANEMOMETRY AND MODEL TESTING

4.1

Introduction

A large part of this chapter is submitted for publication as an article in the Experiments in Fluid

journal. In Chapter 2 and 3, we have obtained an estimation of the wall shear stress generated by

the Bernoulli pad on the opposing workpiece.

The abrupt change in the flow direction introduces separation and recirculation near the neck of

the pad Shi and Li (2016). The separation and constriction of the flow results in a large magnitude

of shear stresses on the workpiece, which is an important parameter for the application of bio-

fouling mitigation. In our previous work Tomar et al. (2022), we have used computational fluid

dynamics simulations to develop a better understanding of the flow physics, including the location

and magnitude of the maximum shear stress. It was found that the magnitude of wall shear stress

is maximum at the belly of the recirculation region. Resolving the region of flow separation is

computationally expensive, which makes it challenging to predict the wall shear stress accurately.

The numerical results need validation, and for the first time in this work, a constant temperature

hot-film anemometer is used with water as the working fluid to measure the wall shear stress.

Over the last few decades, various methods for wall shear stress measurements have been

developed and reviewed in the literature Winter (1979); Fernholz et al. (1996); Sheplak et al.

(2006). The popular methods among researchers include the use of laser Doppler anemometer,

floating element sensor, Preston or Stanton tubes, hot wire or hot-film, etc. Hot-film sensors have

been used to detect flow separation and detection of transition from laminar to turbulent flow - see

Bellhouse and Schultz (1966); Owen (1970); Jiang et al. (2000), for example. In the present work,

we use a flush-mounted hot-film sensor for the measurement of wall shear stress generated by a

Bernoulli pad.

The main principle behind this technique is to correlate heat transfer from the sensor with wall

shear stress. To this end, the sensor needs to be calibrated under known wall shear conditions.

There are two main challenges associated with these experiments. First, velocity gradients are

55

measured at the wall using typical wall shear measurement techniques. The measurement of these

velocity gradients is not possible since the sensor is a flush-mounted sensor, and hence, the sensor

has to be calibrated against the wall shear stress itself. Second, sensor measurements are highly

dependent on ambient conditions, and therefore, calibration must be performed in situ.

Experimental investigations with Bernoulli pads have been reported in the literature. Li and

Kagawa Li and Kagawa (2014), for example, conducted experiments to understand the various

factors that affect the wall normal forces. A majority of the investigations in the literature have

focused on wall-normal forces and used air as the working fluid. To the best of our knowledge,

the only work with water as the working fluid was reported by Kamensky Kamensky (2020);

in this work, Particle Tracking Velocimetry (PTV) experiments were conducted to measure the

velocity components of the flow field. PTV experiments do not provide reliable measurements

close to the wall and cannot be used to accurately measure wall shear stress. Hence, the wall shear

measurements presented in this work fill an important gap in the literature.

This chapter is organized as follows. A brief review of the working principle of hot-film sensors

is provided in Section 4.2. The calibration of the hot-film sensor is carried out using fully-developed

channel flow; an analytical solution for wall shear in channel flow is presented in Section 4.2 and

the procedure for calibration of the hot-film sensor is presented in Section 4.3. The experimental

setup for wall shear measurements is described in Section 4.4. Section 4.5 provides experimental

results and compares them with results obtained from numerical simulations. Concluding remarks

are provided in Section 4.6.

4.2 Background

4.2.1 Analytical solution for wall shear in channel flow

For two-dimensional steady-state fully developed channel flow - see Figure 4.1, the Navier-

Stokes equation for the streamwise velocity becomes:

𝜇

𝑑2𝑢
𝑑𝑦2

=

𝑑𝑝
𝑑𝑥

where 𝑝 = 𝑝(𝑥) is the pressure, 𝑢 = 𝑢(𝑦) is the streamwise velocity, and 𝜇 denotes the dynamic

56

Figure 4.1 Channel setup used for calibration of the hot-film sensor; figures are not drawn to scale: (a)
top-view of the channel showing the eight pressure ports, (b) exploded view of the inlet flow conditioner,
(c) exploded view of section c-c of channel with gasket, polycarbonate sheet, aluminum plate, and clamps
(d) outlet flow conditioner showing internal vanes.

viscosity. On integrating the above equation with respect to 𝑦 once, and substituting the symmetry

condition 𝑑𝑢/𝑑𝑦 = 0 at 𝑦 = ℎ/2, we get:

𝜇

𝑑𝑢
𝑑𝑦

(cid:18)

𝑦 −

=

(cid:19) 𝑑𝑝
𝑑𝑥

ℎ

2

(4.1)

Note that ℎ is the channel’s height. This result does not assume any particular flow regime, laminar

or turbulent. From the above result, the shear stress 𝜏 can be written as:

𝜏 ≜ 𝜇

(cid:18)

𝑦 −

=

𝑑𝑢
𝑑𝑦

(cid:19) 𝑑𝑝
𝑑𝑥

ℎ

2

At the wall, where 𝑦 = 0, the wall shear stress 𝜏w is then given by:

𝜏w = −

ℎ

2

𝑑𝑝
𝑑𝑥

(4.2)

Therefore, the wall shear stress in channel flow can be determined by measuring the pressure

gradient.

57

4.2.2 Working principle of hot-film anemometry

Hot-film anemometry is used to measure the velocity and turbulence properties of fluid flows.

The hot-film sensor used in our experiments is shown in Figure 4.2. The H-shaped film at the tip

of the sensor is a very thin electrical resistor. The main principle behind a hot-film anemometer

is that the heat generated by electric current through the film is dissipated via thermal convection.

The heat transfer rate from the film into the fluid varies with flow velocity according to King’s law

(King, 1914):

and the heat generated by electric current is:

𝑞conv = 𝑎 + 𝑏𝑢𝑛,

𝑞gen =

𝐸 2
a
𝑅f

(4.3)

(4.4)

In (4.3) and (4.4), 𝑎, 𝑏, and 𝑛 are constants, 𝑢 is the characteristic velocity around the hot-film, 𝐸a

is the voltage applied across the film, and 𝑅f is the film’s electrical resistance. To satisfy energy

conservation, 𝑞conv = 𝑞gen. Hence, from equations (4.3) and (4.4):

a = 𝑎𝑅f + (𝑏𝑅f)𝑢𝑛
𝐸 2

(4.5)

Figure 4.2 A schematic of the hot-film sensor 55R46 by Dantec Dynamics Dantec Dynamics (2023).

58

Since the hot-film anemometer will be calibrated with fully developed channel flow, (4.1) can be

integrated with respect to 𝑦, using the non-slip boundary condition 𝑢 = 0 at 𝑦 = 0, to obtain an

expression for 𝑢:

𝜇𝑢 = (𝑦2 − ℎ𝑦)

𝑑𝑝
𝑑𝑥

1
2

Very near the wall, for 𝑦 = 𝛿nw << 1, the linear term dominates and the second order term becomes

negligible.

In turbulent boundary layers, 𝛿nw is equal or less than the thickness of the laminar

sublayer. Under these conditions, the last result simplifies to:

From (4.2), we get:

𝜇𝑢 ≈ −

𝑑𝑝
𝑑𝑥

ℎ

2

𝛿𝑙

𝑢 ≈

𝜏w
𝜇

𝛿nw

Using (4.6), (4.5) can now be expressed in terms of the wall shear stress:

𝐸 2
a = 𝑎𝑅f + 𝑏𝑅f

(cid:19) 𝑛

(cid:18) 𝛿nw
𝜇

𝜏𝑛
w

(4.6)

(4.7)

𝑅f is maintained constant when the sensor is operated with a Constant Temperature Anemometer

(CTA). To this end, the hot-film is implemented within a Wheatstone bridge to constantly monitor

and correct changes in its resistance due to cooling. The main compensation is achieved by changing

the electric current through the film in order to maintaining its temperature constant. Therefore,

the coefficients of equation (4.7) can be absorbed into two calibration constants as follows:

a = 𝐴 + 𝐵𝜏𝑛
𝐸 2

w

,

(4.8)

These constants are determined by correlating the voltage output of the sensor with known wall

shear stress in fully-developed channel flow.

Though (4.8) was obtained without assuming any specific flow regime, heat transfer charac-

teristics in laminar flows differ significantly from those of turbulent flows. In laminar flows, heat

59

Figure 4.3 A CAD describing the channel used for hot-film probe calibration

transfer is diffusion-dominated, whereas it is advection-dominated in turbulent flows because of

transversal motions outside of the laminar sublayer. Therefore, the calibration constants of (4.8)

should be different for laminar and turbulent flows.

4.3 Calibration of Hot-Film Sensor

4.3.1 Design of channel

The channel used for calibration of the hot-film sensor is motivated by prior work reported

in the literature Sun et al. (2018). A schematic of the channel is shown in Figure 4.1(a) and a

sectional view through a pressure port is shown in Figure 4.1(c). The channel was fabricated using

Aluminum and its cross-section (𝑤 = 50 mm, ℎ = 5 mm) was chosen to yield a Reynolds number

of Re = 10,000 with the available pump. A CAD model is presented in the Figure 4.3.

The length of the channel was chosen to be 𝐿 = 1,200 mm to ensure that the flow reaches a

state of fully-developed flow under smooth-wall conditions, and remains so for a significant portion

of the channel. The top of the channel is covered with a clear polycarbonate sheet and leakage is

prevented through the use of a gasket along the length of the channel. An aluminum plate is placed

over the polycarbonate sheet and clamps are used to apply uniform pressure on the gasket along the

length of the channel to make it leak proof - see Figure 4.1(c) and Figure 4.5(a).

A submersible utility pump was used to circulate water through the channel. The flow rate is

controlled with a gate valve installed upstream of the channel’s inlet. The flow is conditioned at

60

Figure 4.4 Pressure measurement using a DP15 differential pressure transducer

the inlet of the channel with a series of baffles, followed by a honeycomb, and then a series of three

meshes of decreasing hole size - see Figure 4.1(b) and Figure 4.5(b). To reduce the developing

length, the boundary layer is tripped with a coarse-grit sand paper strip at the inlet of the channel.

The outlet of the channel is connected to a flow conditioner with four internal vanes to minimize

Figure 4.5 Assembled view of the channel setup in Figure 4.1: (a) top-view of the channel showing sensor
mount (without sensor) and eight pressure ports (see Figure 4.1) connected via three-way valves, (b) inlet
flow conditioner, (c) hot-film sensor, (d) magnified view of sensor mount with sensor, (e) sectional view of
channel setup through the sensor mount and sensor.

61

Figure 4.6 Pressure at the eight different ports of the channel (see Figure 4.1), computed based on pressure
differential measurements relative to port 1 and assignment of an arbitrary pressure to port 8 . Note that
the straight line fit was obtained by using the data from ports 4 though 8 .

end-effects - see Figure 4.1(d).

To measure pressure along the channel, eight pressure taps are placed on the side of the channel

along its length. The ports are drilled with a size of 0.8 mm. The ports are connected by a three-

way valve - see Figure 4.5(a), to a DP15 differential pressure transducer (a product of Validyne

Engineering VALIDYNE Engineering (2022)), as shown in Figure 4.4. The pressure differential

between ports 2 through 8 relative to 1 are recorded at steady state; the data is then used to

compute the pressure at all the ports by assigning an arbitrary pressure to port 8 - see Figure

4.6. The plot indicates that fully developed flow is achieved early and well ahead of port 1 . The

hot-film sensor can be mounted anywhere along the length of the channel in which the pressure drop

exhibits a linear behavior, because that is the region where the flow is fully developed. We chose

to mount the sensor, shown in Figure 4.5(c), between ports 4 and 5 on the polycarbonate sheet -

see Figure 4.1(a), Figure 4.5(a) and Figure 4.5(d). A sectional view of the assembled channel setup

through the sensor mount and sensor is shown in Figure 4.5(e).

62

4.3.2 Calibration procedure

The hot-film sensor 55R46, shown in Figure 4.2 and Figure 4.5(c), was calibrated using the

constant temperature anemometer MiniCTA 54T42, also a product of Dantec Dynamics Dantec

Dynamics (2023). The sensor is used with water as the working fluid and therefore the jumpers

inside the anemometer were adjusted according to the desired value of overheat and sensor resis-

tance. The sensor is flush-mounted on top of the polycarbonate sheet between pressure ports 4

and 5 , where the flow is fully developed - see our discussion in Section 4.3.1.

The calibration of the sensor was performed using 21 different flow rates through the channel;

a gate valve was installed between the prime mover and the inlet of the channel to control the flow

rate. For each flow rate, the pressure transducer and hot-film sensor measurements (voltages) were

recorded. The pressure measurement (voltage) provides the pressure drop between ports 4 and 5 ;

transducer calibration data is used to express it in Pa and the wall shear stress 𝜏w is then computed

using Eq.(4.2). The variation of the square of the hot-film sensor voltage (𝐸 2

a ) with the wall shear
stress 𝜏w is plotted in Figure 4.7 with the objective of computing the calibration coefficients in

Eq.(4.8). It can be seen that the wall shear stress increases monotonically with increase in the mass

flow rate. Also, the variation of 𝐸 2

a with 𝜏w depicts three distinct behaviors corresponding to three

distinct flow regimes: turbulent, critical and laminar. In the laminar and turbulent regimes, the data

points show a concave downward trend with increase in wall shear stress; the trend is reversed and

is concave upward in the critical regime.

The data points corresponding to the laminar and turbulent regimes are used to find a linear

a and 𝜏𝑛

fit between 𝐸 2

w, where the value of the exponent 𝑛 should lie in the range [0.1, 0.5]
Cameron Tropea (2007). A linear fit is found by choosing 𝑛 = 0.4 for the data in the laminar regime

and 𝑛 = 0.2 for the data in the turbulent regime. The plots are provided in Figure 4.8 and the

calibration coefficients in Eq.(4.8) are provided in Table 4.1; these coefficient were obtained with

an 𝑅2 value of 0.99 and a confidence level of 95%. The calibration coefficients in Table 4.1 will be

used to measure the wall shear stress generated by a Bernoulli pad in Section 4.5.

63

Figure 4.7 Calibration data showing the variation of 𝐸 2

a with 𝜏w

Table 4.1 Calibration coefficients and exponent for laminar and turbulent flow regimes

Flow regime
Laminar
Turbulent

𝐴
-0.216
-50.86

𝐵
0.798
32.70

𝑛
0.4
0.2

4.4 Bernoulli Pad Experimental Setup

An experimental setup was developed to measure the wall shear stress generated by a Bernoulli

pad. An important component of the setup is the Bernoulli pad assembly, comprised of the stem,

pad, and flow conditioning section - see Figure 4.9(a) and (b). The stem is a tube with an inside

diameter of 𝑑 = 25.4 mm (1.0 in), an outside diameter of 31.75 mm (1.25 in), and a length of 88.9

mm (3.5 in). The pad is a circular, flat plate with a diameter of 𝐷 = 203.2 mm (8.0 in) and a

thickness of 6.35 mm (0.25 in). A 25.4 mm (1.0 in) hole is located at the center of the pad along

with a 31.75 mm (1.25 in) counterbore to 75% depth in the pad. The stem interfaces with the

pad in the counterbore such that there are negligible steps for the fluid to encounter. Cast acrylic

was selected as the material for both the stem and the pad; its rigidity is adequate for the current

64

Figure 4.8 Linear calibration curves for laminar and turbulent regimes

purpose and its surface roughness is very low.

In particular, polished acrylic sheets typically

exhibit a surface roughness within 0.14 𝜇m < Ra < 0.23 𝜇m S et al. (2011), or equivalently,

0.62 𝜇m < Rz < 0.89 𝜇m Stähli (2013). Here, Ra is the arithmetic mean of the deviation from the

roughness centerline, and Rz is the mean roughness depth.

The flow conditioning section is used to obtain a low-turbulence-intensity and top-hat velocity

profile at the inlet of the stem. It is comprised of the following elements from upstream to down-

stream order: an inlet piece, a baffle that breaks down the incoming pipe flow, a honeycomb, three

wire meshes of decreasing hole size, and a flow contraction - see Figure 4.9(b). A cubic equation

was used to design the flow contraction profile Hussain and Ramjee (1976). The components

of the flow conditioning section were mated with the stem in a similar fashion as the stem-pad

interface to avoid disturbances in the flow. Except the inlet, baffle, honeycomb, and wire meshes,

the materials of the components of the flow conditioning section were chosen as cast acrylic and

polyurethane-coated 3D prints (fabricated with PLA - polylactic acid) for smooth surface finishes.

65

Figure 4.9 (a) Bernoulli pad assembly shown in its nominal configuration over the workpiece (b) exploded
view of Bernoulli pad assembly (c) sectional view of flush-mounted sensor in workpiece.

The primary components of the experimental setup are: workpiece, sensor mount (a detailed

CAD and a view below the workpiece in Bernoulli pad assembly is provided in Figure 4.10 and

Figure 4.11 respectively), linear stage, Bernoulli pad assembly mount and frame. An assembled

view of these components, which we will refer to as the Shear Test Station (STS), is shown in

Figure 4.12. A rectangular, flat plate of 7075-T6 aluminum was selected as the workpiece because

of its rigidity and electrical conductivity.

It has dimensions of 406.4 mm × 304.8 mm × 6.35

mm, and the surface of the workpiece was sanded and polished to reduce its surface roughness.

This process typically results in roughness levels in the range of 0.1 𝜇m < Ra < 0.4 𝜇m or

0.5 𝜇m < Rz < 1.4 𝜇m. Located at the center of the workpiece width and 254 mm down its

length is a counterbored hole that allows the hot-film sensor to mate with the workpiece such that

the sensing surface of the sensor is flush with the workpiece - see Figure 4.9(c). A grounding

66

Figure 4.10 A detailed cad of the sensor mount installed in the Bernoulli pad assembly

Figure 4.11 View of the sensor mount installed below the Bernoulli pad assembly

67

Figure 4.12 An assembled view of the Shear Test Station (STS), comprised of the Bernoulli pad assembly
mount, workpiece and linear stage.

lead necessary for the hot-film sensor was permanently fastened to one of the plate’s vertical faces.

Blind holes allow the sensor mount to be secured to the underside of the workpiece as shown in

Figure 4.13. The probe mount in Figure 4.13, contains steps in its internal geometry that constrain

the hot-film probe when it is fastened to the underside of the work surface.

A LTS300 linear translation stage manufactured by THORLABS (THORLABS (2023)) was

used to move the Bernoulli pad over the length of the workpiece such that the hot-film sensor can

measure the wall shear stress along the radial direction. The LTS300 has an absolute accuracy of

50 𝜇m and a repeatability of 2 𝜇m. The Bernoulli pad assembly mount, shown in Figure 4.12,

is a series of structural aluminum members that interface the Bernoulli pad assembly in Figure

4.9(a) to the linear stage. The joints of the Bernoulli pad assembly mount are comprised of slots

and threaded fasteners that allow for setting a uniform gap between the pad and the workpiece.

Structural aluminum extrusions form a frame that allows the workpiece to be positioned below

the Bernoulli pad assembly mount and ensure that the pad and the workpiece remain parallel with

68

Figure 4.13 The work surface, hot-film probe, and probe mount

travel of the linear stage. The fluid power assembly consists of a 250 W submersible pump, rubber

hoses, and a gate valve for flow rate control.

4.5 Wall Shear Stress Measurement

4.5.1 Procedure

The gap height between the pad and the workpiece was chosen to be 𝐻 = 1.3 mm for our

experiments. Shims were placed in the gap to ensure that the gap height was uniform; they were

removed after ensuring uniform desired gap height. The Shear Test Station (STS) was then placed

into a water tank with dimensions of 2120 mm × 749 mm × 762 mm - see Figure 4.14. The frame

suspended the workpiece in the center of the tank footprint and 356 mm from the bottom. The

dimensions of the tank and the location of the workpiece were chosen to ensure that the radial

outflow between the pad and workpiece was not affected by the water splashing from the tank

walls. Shims were placed between the frame and tank to level the STS before it was secured to

the tank with clamps. The fluid power assembly was connected, and the mass flow rate was set to

(cid:164)𝑚 = 0.046 kg/s. The mass flow rate was chosen such that the maximum wall shear stress was within

the range of the calibration setup. Sufficient water was added to the tank without submerging the

workpiece and the water was allowed to adjust to the ambient temperature (17.25◦ C) before data

collection began. The density and viscosity of water at this temperature are 𝜌 = 998.87 kg/m3 and

69

Figure 4.14 (a) A perspective view of the experimental setup (b) a side view showing the STS

𝜇 = 0.001073 Pa s.

The pump was powered on, and the system was allowed to reach steady-state. The linear stage

was used to move the Bernoulli pad from its nominal position, shown in Figure 4.9(a), in steps with

an average size of 5 mm, which is approximately twice the diameter of the hot-film sensor. The

step size was reduced to 1.0-2.0 mm close to the neck of the Bernoulli pad, where large variations

in the shear stress are expected. The total distance of travel was approximately ℓ = 83.0 mm; this

ensured that the sensor was never exposed to air. After each step, the sensor was powered on, and

the hot-film sensor output was obtained using an oscilloscope. After recording the data, the sensor

was powered off to avoid unnecessary heating. A view of the Bernoulli pad assembly during the

data collection is provided in the Figure 4.15.

The voltage outputs from the hot-film sensor were first corrected for temperature difference

between calibration and wall shear stress measurement conditions.

In particular, we used the

following relation for temperature correction:

𝐸corr =

(cid:18)𝑇f − 𝑇0
𝑇f − 𝑇a

(cid:19) 0.5

𝐸a

70

Figure 4.15 A side view showing the STS during the data collection process

where, 𝐸corr is the corrected voltage, 𝑇0 = 25.49◦C is the water temperature during calibration,

𝑇f = 40◦C is the film temperature, and 𝑇a = 17.25◦C is the water temperature during data acquisition.

The corrected voltage values were used to obtain the wall shear stress values using Eq.(4.8) with

the calibration coefficients provided in Table 4.1.

4.5.2 Comparison with numerical results

Numerical simulations were carried out using five different turbulence models, namely, Spalart-

Allmaras, 𝜅-𝜀, 𝜅-𝜔, Transition-SST, and RSM. These simulations were carried out for identical

flow domain and boundary conditions in the experiments - see Figure 4.16. The domain is axially

symmetric without variations in the azimuthal directions. Therefore, we used a two-dimensional

axisymmetric model to reduce computational time. Assuming incompressible flow, we imposed

a mass flow rate of (cid:164)𝑚 = 0.046 kg/s at the inlet, and exit pressure 𝑝 = 𝑝atm at the pad’s outlet.

Non-slip conditions were imposed on all solid walls.

Turbulence models implemented in commercial software require us to prescribe surface rough-

ness level in solid walls. The choice of surface roughness level could be very relevant in fully rough

flow, moderately important in transitionally rough flow, or completely unimportant in hydraulically

smooth flow. The maximum wall shear stress measured in the current experiments was 𝜏w, max =

71

Figure 4.16 A schematic of the Bernoulli pad used in simulations.

10.19 Pa. Hence, the maximum friction velocity was 𝑢𝜏max
thickness of the laminar sublayer in a turbulent boundary layer is 𝛿+

≜ √︁𝜏w,max/𝜌 = 0.1 m/s. Given that the
𝑙 ≜ 𝜌𝛿𝑙𝑢𝜏/𝜇 = 5, our flow would
𝑙 𝜇/(𝜌𝑢𝜏) = 49.5
𝜇m. As mentioned before, the roughness level in polished aluminum goes typically up to Rz = 1.4

be hydraulically smooth if its surface roughness’ representative height 𝜀w < 𝛿𝑙 = 𝛿+

𝜇m. Therefore, the flow through the pad’s gap is hydraulically smooth. For the simulations, any

choice of roughness height for which 𝜀+ < 5 will have no roughness effects. With this in mind,

the surface roughness of the pad and the workpiece were assumed to be 𝜀p = 𝜀w = 1 𝜇m to ensure

hydraulically smooth conditions as in the experiments.

In our experiments, the hot-film sensor voltage is obtained for the measurement region of length

ℓ = 83.0 mm shown in Figure 4.16. A comparison of these voltages with the voltage values in

Figure 4.7 indicate that the flow is laminar for 𝑟 ∈ [0, 12] ∪ [14, 83] and turbulent for 𝑟 ∈ [12, 14],

where the units are in mm. The change from laminar to turbulent flow close to the neck of the

pad is expected due to the constriction of the flow and the formation of a recirculation region due

to flow separation upon change in the direction of flow from axial to radial around a sharp corner.

To obtain the wall shear stress, we use the calibration coefficients in Table 4.1: the exponent and

coefficient for laminar flow are used for 𝑟 ∈ [0, 12] ∪ [14, 83] and those for turbulent flow are

used for 𝑟 ∈ [12, 14]. The wall shear stress and radial distance are non-dimensionalized using the

expressions in Anshul (Tomar et al. (2022)):

¯𝜏w =

(cid:19)

(cid:18) 𝜋𝜌𝑑3
4 (cid:164)𝑚𝜇

𝜏w,

¯𝑟 =

(cid:19)

(cid:18) 2
𝐷

𝑟

(4.9)

72

Figure 4.17 Variation of ¯𝜏w with ¯𝑟: a comparison of simulation and experimental results. The
experimental results, which were obtained at discrete values of ¯𝑟, are shown using ∗ marks. The top figure
shows the recirculation region and the flow around it Shi and Li (2016) - see figure inset. A magnified view
of the dotted portion of the top figure is shown in the bottom figure for comparison of the five different
turbulence models with experimental results.

The variation in the non-dimensional wall shear stress 𝜏w with non-dimensional radial distance

𝑟 is shown in Figure 4.17. No error bars are shown in these experimental results as the uncertainty

associated with the measurements was found to be negligible (the extent of the error bars is less

than the line thickness in the plot). For the purpose of model testing, the numerical simulation

results from the five different turbulence models are also presented in Figure 4.17. All models

coincide with each other for 0 ≤ 𝑟 ≲ 0.14 and 0.5 ≲ 𝑟 ≤ 1.0, but differ significantly from each

other in the reattachment region.

Overall, all models follow the qualitative behavior of the experimentally determined wall shear

stress after the reattachment region. Nevertheless, all five turbulence models over-predict the wall

73

Table 4.2 Dimensionless maximum wall shear stress values obtained from experiment, numerical
models, and power law Tomar et al. (2022)

Experiment
2656.1

Spalart-Almaras
2673.1

𝜅-𝜀
2741.4

𝜅-𝜔
2692.5

Transition-SST
2676.3

RSM Power Law
2696.9

2548.6

shear stress values in this region. This is consistent with the experimental results from Tomar

et al. (2022). All numerical simulations predict both the magnitude and location of the maximum

wall shear stress adequately, suggesting that all of them are able to capture the recirculation region

reasonably well. The maximum non-dimensional wall shear stress 𝜏w,max for all cases, including a

power law estimation from Tomar et al. (2022), is shown in Table 4.2. The simulation results and

the power-law prediction exhibit good agreement with the experimental results.

The experimental measurements show a secondary peak in wall shear stress right after the sharp

drop from the maximum wall shear stress. This secondary peak is followed by a gradual decay in

wall shear stress. This gradual decay in wall shear stress is easily explained by the 1/𝑟 proportional

decay in radial velocity as the cross-sectional area of the pad gap increases with radius Guo et al.

(2017). The secondary peak on the other hand is the result of reattachment phenomena. That

is, the flow initially “squeezes" underneath the separation bubble forming a vena contracta and

thinning the boundary layer. At the vena contracta, the flow is the fastest and the boundary layer

is the thinnest, inducing the highest level of wall shear stress. The flow expands relatively rapidly

toward the pad as it flows past the separation bubble. In this region, there is a localized change of

direction as shown in the inset of Figure 4.17. Here, part of the linear momentum transfers from the

radial direction to the wall-normal direction, with a resulting reduction in radial momentum. This

results in the sharp drop of wall shear stress right after the vena contracta. But shortly after that,

the momentum recovers its purely radial direction, and since the boundary layer is still very thin,

the shear stress experiences a small boost. This is likely what causes the secondary peak in wall

shear stress. Right after that, the boundary layer starts growing and the cross-section continues to

grow with radius, causing the subsequent gradual reduction on wall shear stress.

Most of the models fail to capture the wall shear behavior in the reattachment region. The

Spalart-Almaras model is the only scheme that exhibits the same qualitative behavior as the

74

experimentally measured wall shear stress, though the secondary peak is located slightly further

downstream than that in experiments. This result suggests that, in this model, the simulated

streamwise extent of the separation bubble downstream of the vena contracta is larger than that in

experiments. That the Spalart-Almaras performs better than the other models in the reattachment

region is likely because this model performs well in wall-bounded flows at moderate to low Reynolds

numbers under adverse pressure gradients.

All five turbulence models slightly over-predict the wall shear stress values for large values of

𝑟. This is consistent with the comparison between simulations and Particle Tracking Velocimetry

measurements in Tomar et al. (2022). A better turbulence scheme, such as LES or DNS, may

provide a better match for the entire domain of 𝑟; however, this will require significantly higher

computational effort and lies in the scope of future research. The most important aspect of a

Bernoulli pad for grooming applications is its maximum wall shear stress; the details of the

reattachment region are of relatively low importance. Henceforth, the results shown in Figure

4.17 suggest that our numerical solvers are appropriate to capture the most relevant features of

flow through a Bernoulli pad, with the Spalart-Almaras scheme exhibiting a slight advantage over

the others. It is then adequate to use this numerical model to optimize Bernoulli pad designs for

improving grooming efficacy in non-contact biofouling mitigation.

4.6 Conclusion

The experimental work presented here uses a constant temperature anemometer with a hot-film

sensor to quantify the wall shear stress generated by the action of Bernoulli pad over a proximally

located workpiece. An experimental setup, consisting of a rectangular channel, is designed to

calibrate the wall shear stress. The calibration of the sensor is carried out separately for laminar

and turbulent regimes. These calibration relations are subsequently used to measure the wall shear

stress generated by a Bernoulli pad. It should be mentioned that this experimental effort, which

quantifies the wall shear stress generated by a Bernoulli pad with water as the working fluid, is the

first of its kind.

Except for the reattachment region, numerical simulations exhibit good agreement with hot-

75

film wall shear measurements. The position of the maximum wall shear stress is found to be very

close to the neck of the Bernoulli pad, right below the belly of recirculation region. The predicted

position and magnitude of the maximum wall shear stress are very close to the measured values.

All tested computational models match each other beyond the reattachment region and capture the

general trend of the wall shear stress, albeit slightly overestimating its value. The Spalart-Almaras

model reasonably captures the wall shear-stress behavior in the reattachment region, while the rest

of the models fail to do so.

76

CHAPTER 5

SUMMARY AND FUTURE WORK

5.1 Summary

Numerical simulations have been carried out to develop a better understanding of the flow

physics of the Bernoulli gripper. The objective of the work is to employ the Bernoulli pad as a

shear-based grooming tool. Computational Fluid Dynamics, CFD simulations are carried out for a

range of parameters in the domain to achieve this. A power-law is obtained between the inlet fluid

power and maximum wall shear stress obtained at the workpiece. The power-law collapses the data

from all the parameter space on a single curve.

The data also indicate a power-law between the non-dimensional maximum shear stress and

inlet Reynolds number as: 𝜏∗

𝑤,𝑚𝑎𝑥 ∝ Re𝑛+1

𝑖𝑛 . The value of 𝑛 obtained from a direct fit into the

CFD data is broadly consistent with the value of 𝑛 found in the literature. Near the area where the

maximum shear stress occurs, the velocity and gradient are highest, and this spot is closer to the jet

exit 𝑑 rather than the outer diameter 𝐷 of the pad. According to the discussion, the power required

to counteract wall shear across the entire pad seems directly related to its maximum value.

The direct relationship between the inlet fluid power and maximum wall shear stress motivated

us to explore other ways of delivering the power to the system. For the first time in this work, the

Bernoulli pad is provided with rotation to increase the wall shear stress and the normal force. The

effect of rotation on the pad is first analyzed without inlet fluid power. The addition of inlet fluid

power, along with the rotation, results in stable equilibrium configurations. It was noticed that there

are two mass flow rate values for which the force generated on the workpiece is the same for a given

equilibrium height. Therefore, two equilibrium configurations exists, which in this work has been

referred as bi-stiffness equilibria. The stiffness values at these equilibrium points are different; the

stiffness is higher for the higher mass flow rates and vice versa.

The addition of rotation to the Bernoulli pad results in an increase in suction force. The work

analysis was carried out, and it was found that for the same inlet fluid power, the increase in

rotational speed from 0 rad/s to 60 rad/s results in a 14.93 % rise in the value of maximum wall

77

shear stress and a 13.85 % rise in the value of average wall shear stress. Similarly, for the same

inlet fluid power, the normal force increases from 1.4037 N (repulsive and stable) to -6.6827 N

(attractive and stable).

The numerical solver used in the present work is validated using the hot-film experiments. A

setup is designed to calibrate the hot-film sensor with a Constant Temperature Anemometer, CTA.

The calibration setup consists of a channel with flow conditioning parts installed at the inlet and

outlet. Flow conditioning parts are used to make sure that the flow is fully developed in the channel.

A fully developed flow in the channel ensures that a simplified analytical expression can be used to

calibrate the hot-film sensor. The calibration curve obtained using the channel depicts three distinct

regions: laminar, transition, and turbulent. Calibration coefficients are obtained for laminar and

turbulent regions.

The calibrated sensor is used to obtain wall shear stress in a Bernoulli pad assembly on a shear

stress station. The experimental setup, which consisted of a Bernoulli pad assembly and a shear

stress station, was designed to ensure that the flow at the inlet of the stem had a low-turbulence

intensity with a top-hat velocity profile. The wall shear stress obtained from the experiments is

compared with the results obtained through CFD simulations for different turbulence models. The

results from the experiments and numerical simulations show a very good match for the peak shear

stress value. This suggests that all the numerical models are able to capture the recirculation region

really well.

It was also observed that most of the models are unable to capture the wall shear

behavior in the reattachment region. Among all the turbulence models, the Spalart-Almaras model

performs better than the other turbulence models.

5.2 Future Scope

In this study, water was selected as the working fluid to simulate the flow field of the Bernoulli

pad. This choice aligns with the intended application, which involves utilizing the Bernoulli pad

to groom underwater surfaces without physical contact. By opting for water as the working fluid,

the study aimed to replicate real-world conditions accurately. However, it’s worth noting that the

decision to use water as the working fluid introduces potential complexities, particularly regarding

78

cavitation effects. Cavitation occurs when local pressure drops below the vapor pressure, leading

to the formation of vapor bubbles. These bubbles tend to form in regions where local velocity

is high, precisely where the highest wall-shear stress is observed. If not avoided, cavitation can

potentially harm the workpiece to be cleaned by resulting in erosion and pitting. The study did

not explore cavitation effects, but it acknowledges the importance of addressing them in future

research. Specifically, future work will focus on designing alternative flow geometries that can

generate sufficiently high levels of shear stress to groom surfaces effectively while also avoiding

cavitation. This aspect opens up avenues for further investigation and optimization of Bernoulli

pad applications in underwater surface grooming.

In the present work, we have designed an experimental setup to measure the wall shear stress

for a standard non-rotating Bernoulli pad and the numerical model is validated with five different

turbulence models. In the future, an experimental setup needs to be designed to measure the wall

shear stress for rotating Bernoulli pad. We expect this to be a challenging experiment as the gap

between the rotating pad and the stationary workpiece is very small and a small misalignment could

result in pad damaging against the workpiece. The CFD simulations in the present work for rotating

pads are limited to a maximum rotational speed, 𝜔 = 60 rad/s. The value of rotational speed was

chosen to avoid the convergence issues associated with the steep velocity and pressure gradients in

the domain. However, we acknowledge that in future to fully maximize the benefits of the rotating

pad, we need to simulate the rotating Bernoulli pad with higher rotational speeds.

79

BIBLIOGRAPHY

ANSYS (2022). Ansys fluent: Theory guide 2022 R2 (Section 24.4.3.3. Coupled Algorithm).

Bellhouse, B. J. and Schultz, D. L. (1966). Determination of mean and dynamic skin friction,
separation and transition in low-speed flow with a thin-film heated element. Journal of Fluid
Mechanics, 24(2):379–400.

Brun, X. and Melkote, S. (2009). Analysis of stresses and breakage of crystalline silicon wafers

during handling and transport. Solar Energy Materials and Solar Cells, 93:1238–1247.

Cameron Tropea, Alexander L. Yarin, J. F. F. (2007). Springer Handbook of Experimental Fluid

Mechanics. Springer Berlin, Heidelberg.

D’Agostino, R. B. and Stephens, M. A. (1986). Goodness-of-Fit Techniques. Marcel Dekker, Inc.,

USA.

Dahlgren, J., Foy, L., Hunsucker, K., Gardner, H., Swain, G., Stafslien, S. J., Vanderwal, L., Bahr,
J., and Webster, D. C. (2022). Grooming of fouling-release coatings to control marine fouling
and determining how grooming affects the surface. Biofouling, 38(4):384–400.

Dantec Dynamics (2023). Flush-mounted hot-film probe.

Dini, G., Fantoni, G., and Failli, F. (2009). Grasping leather plies by bernoulli grippers. CIRP

Annals, 58(1):21–24.

Donnelly, B., Sammut, K., and Tang, Y. (2022). Materials selection for antifouling systems in

marine structures. Molecules, 27(11):3408.

Erturk, S. and Erzincanli, F. (2020). Design and development of a non-contact robotic gripper
for tissue manipulation in minimally invasive surgery. Acta bio-medica : Atenei Parmensis,
91:e2020071.

Fantoni, G., Santochi, M., Dini, G., Tracht, K., Scholz-Reiter, B., Fleischer, J., Kristoffer Lien, T.,
Seliger, G., Reinhart, G., Franke, J., Nørgaard Hansen, H., and Verl, A. (2014). Grasping devices
and methods in automated production processes. CIRP Annals, 63(2):679–701.

Fernholz, H. H., Janke, G., Schober, M., Wagner, P. M., and Warnack, D. (1996). New developments
and applications of skin-friction measuring techniques. Measurement Science and Technology,
7(10):1396.

Frey, H. (U.S. Patent 5967578, 1999). Tool for the contact-free support of plate-like substrates.

Gibson, L., Galloway, L., Kim, S. i., and Spence, S. (2017). Assessment of turbulence model
predictions for a centrifugal compressor simulation. Journal of the Global Power and Propulsion

80

Society, 1:142–156.

Guo, J., Shan, H., Xie, Z., Li, C., Xu, H., and Zhang, J. (2017). Exact solution to Navier-Stokes
equation for developed radial flow between parallel disks. Journal of Engineering Mechanics,
143(6).

Houghtalen, R. J., Akan, A. O., and Hwang, N. H. (1996). Fundamentals of hydraulic engineering

systems. Prentice Hall.

Hu, P., Xie, Q., Ma, C., and Zhang, G. (2020). Silicone-based fouling-release coatings for marine

antifouling. Langmuir, 36(9):2170–2183.

Hussain, A. K. M. F. and Ramjee, V. (1976). Effects of the Axisymmetric Contraction Shape on

Incompressible Turbulent Flow. Journal of Fluids Engineering, 98(1):58–68.

Islam, S. M., Khan, M. T., and Ahmed, Z. U. (2020). Effect of design parameters on flow
characteristics of an aerodynamic swirl nozzle. Progress in Computational Fluid Dynamics, an
International Journal, 20(5):249–262.

Jawarneh, A. M. and Vatistas, G. H. (2006). Reynolds Stress Model in the Prediction of Confined

Turbulent Swirling Flows. Journal of Fluids Engineering, 128(6):1377–1382.

Jiang, F., Lee, G.-B., Tai, Y.-C., and Ho, C.-M. (2000). A flexible micromachine-based shear-stress
sensor array and its application to separation-point detection. Sensors and Actuators A: Physical,
79(3):194–203.

Kamensky, K., Hellum, A., and Mukherjee, R. (2019). Power scaling of radial outflow: Bernoulli

pads in equilibrium. ASME Journal of Fluids Engineering, 141.

Kamensky, K., Hellum, A., Mukherjee, R., Naik, A., and Moisander, P. (2020). Underwater shear-
based grooming of marine biofouling using a non-contact bernoulli pad device. Biofouling,
36:951–964.

Kamensky, K. M. (2020). A new paradigm for generating surface-normal forces for hull-cleaning
robots. PhD dissertation, Department of Mechanical Engineering, Michigan State University.

King, L. V. (1914). Xii. on the convection of heat from small cylinders in a stream of fluid:
Determination of the convection constants of small platinum wires with applications to hot-wire
anemometry. Philosophical transactions of the royal society of London. series A, containing
papers of a mathematical or physical character, 214(509-522):373–432.

Li, X. and Kagawa, T. (2014). Theoretical and experimental study of factors affecting the suction

force of a bernoulli gripper. Journal of Engineering Mechanics, 140:04014066.

Logue, J. (U.S. Patent No. 3523706, 1970). Apparatus for supporting articles without structural

81

contact and for positioning the supported articles.

McIlwraith, L. and Christie, A. (U.S. Patent No. 6601888, 2003). Contactless handling of objects.

Menesses, M., Belden, J., Dickenson, N., and Bird, J. (2017). Measuring a critical stress for
continuous prevention of marine biofouling accumulation with aeration. Biofouling, 33(9):703–
711.

Misimi, E., oye, E., Lillienskiold, A., Mathiassen, J., Berg, O. A., Gjerstad, T. B., Buljo, J.,
and Skotheim, O. (2016). GRIBBOT – Robotic 3D vision-guided harvesting of chicken fillets.
Computers and Electronics in Agriculture, 121:84–100.

Nakabayashi, K., Ichikawa, T., and Morinishi, Y. (2002). Size of annular separation bubble around
the inlet corner and viscous flow structure between two parallel disks. Experiments in Fluids,
32(4):425–433.

Neto, J. L. V., Martins, A. L., Neto, A. S., Ataíde, C. H., and Barrozo, M. A. S. (2011). CFD applied
to turbulent flows in concentric and eccentric annuli with inner shaft rotation. The Canadian
Journal of Chemical Engineering, 89(4):636–646.

Nikuradse, J. (1950). Laws of flow in rough pipes, ational advisory commission for aeronautics
(naca) technical memorandum 1292, (english translation of "stromungsgesetze in rauhen rohren."
vdi-forschungsheft 361, 1933.

Oliveira, D. R. and Granhag, L. (2020). Ship hull in-water cleaning and its effects on fouling-control

coatings. Biofouling, 36(3):332–350.

Olsson, R. G. and Williams, E. C. (U.S. Patent No. 3,438,668, 1969). Contactless lifter.

Owen, F. K. (1970). Transition experiments on a flat plate at subsonic and supersonic speeds.

Paivanas, J. A. and Hassan, J. K. (1981). Attraction force characteristics engendered by bounded,

radially diverging air flow. IBM Journal of Research and Development, 25(3):176–186.

Phares, D. J., Smedley, G. T., and Flagan, R. C. (2000). The wall shear stress produced by the
normal impingement of a jet on a flat surface. Journal of Fluid Mechanics, 418(1):351–375.

S, S., Nair, C. K., and Shetty, J. (2011). Effect of Different Polishing Agents on Surface Finish
and Hardness of Denture Base Acrylic Resins: A Comparative Study. International Journal of
Prosthodontics and Restorative Dentistry, 1(1):7–11.

Sheplak, M., Cattafesta, L., Nishida, T., and Mcginley, C. (2006). Mems shear stress sensors:

Promise and progress. IUTAM Symposium on Flow Control and MEMS, pages 67–76.

Shi, K. and Li, X. (2016). Optimization of outer diameter of Bernoulli gripper. Experimental

82

Thermal and Fluid Science, 77:284–294.

Song, C. and Cui, W. (2020). Review of underwater ship hull cleaning technologies. Journal of

marine science and application, 19(3):415–429.

Stähli, A. (2013). The technique of lapping. Pieterlen/Biel.

Stephens, M. A. (1974). EDF statistics for goodness of fit and some comparisons. Journal of the

American Statistical Association, 69(347):730–737.

Sun, B., Wang, P., Luo, J., Deng, J., Guo, S., and Ma, B. (2018). A flexible hot-film sensor array

for underwater shear stress and transition measurement. Sensors, 18(10).

THORLABS (2023). Lts300: 300 mm linear translation stage with integrated controller, stepper

motor.

Tomar, A. S., Kamensky, K. M., Mejia-Alvarez, R., Hellum, A. M., and Mukherjee, R. (2022). A
scaling relationship between power and shear for Bernoulli pads at equilibrium. Flow, 2:E29.

Tribou, M. and Swain, G. (2010). The use of proactive in-water grooming to improve the perfor-

mance of ship hull antifouling coatings. Biofouling, 26(1):47–56.

Trickey, J. S., Cárdenas-Hinojosa, G., Rojas-Bracho, L., Schorr, G. S., Rone, B. K., Hidalgo-Pla, E.,
Rice, A., and Baumann-Pickering, S. (2022). Ultrasonic antifouling devices negatively impact
Cuvier’s beaked whales near Guadalupe island, México. Communications Biology, 5(1):1005.

VALIDYNE Engineering (2022). Dp15 variable reluctance pressure sensor capable of range

changes.

Wagner, M., Chen, X., Nayyerloo, M., Wang, W., and Chase, J. G. (2008). A novel wall climbing
robot based on bernoulli effect. In IEEE/ASME International Conference on Mechtronic and
Embedded Systems and Applications, pages 210–215.

Wark, C. E. and Foss, J. F. (1984). Forces caused by the radial out-flow between parallel disks.

Journal of Fluids Engineering-transactions of The Asme, 106:292–297.

White, F. (1974). Viscous fluid flow. NY: McGraw-Hill, New York.

Winter, K. G. (1979). An outline of the techniques available for the measurement of skin friction

in turbulent boundary layers. Progress in Aerospace Sciences, 18:1–57.

Yeginbayeva, I. and Atlar, M. (2018). An experimental investigation into the surface and hydro-
dynamic characteristics of marine coatings with mimicked hull roughness ranges. Biofouling,
34(9):1001–1019.

83

APPENDIX

For this analysis, we will assume a system of cylindrical coordinates, (𝑟, 𝜃, 𝑧), with the 𝑧-axis

coinciding with the stem axis as shown in Figure 1.2, and the 𝑟 − 𝜃 plane coinciding with the

impingement wall. Further, we will assume that the origin of this coordinate system coincides with

the origin of a reference coordinate system, (𝑅, Θ, 𝑍), whose 𝑍-axis is aligned with and is positive

in the opposite direction of gravity. For simplicity, and without loss of generality, the system

(𝑟, 𝜃, 𝑧) is rotated with respect to the system (𝑅, Θ, 𝑍) around the reference axis (Θ = 0, 𝑍 = 0) by

an angle 𝛼 with respect to the 𝑍-axis. Finally, the reference axis (𝜃 = 0, 𝑧 = 0) makes an angle 𝜙

with the reference axis (Θ = 0, 𝑍 = 0).

In general, for any angle 𝛼 between the 𝑧- and 𝑍-axes, the gravitational field in the (𝑟, 𝜃, 𝑧)

coordinate system is given by:

g = −𝑔 sin 𝛼 sin(𝜃 + 𝜙) ˆr − 𝑔 sin 𝛼 cos(𝜃 + 𝜙) ˆθ − 𝑔 cos 𝛼ˆz

(A.1)

Note that in our simulations 𝛼 = 𝜋/2 and 𝜙 = −𝜋/2. Hence, the gravitational field in the

(𝑟, 𝜃, 𝑧) system is:

g = 𝑔 cos 𝜃 ˆr − 𝑔 sin 𝜃 ˆθ + 0ˆz

(A.2)

All the analysis will be carried out based on the more general form shown in equation (A.1),

and then the result will be simplified for the specific case shown in equation (A.2).

Specific Considerations

Under this condition, the outside hydrostatic pressure would change azimuthally at the pad’s

discharge due to changes in water depth as 𝜃 changes. This boundary condition would likely induce

a non-zero azimuthal change in the radial velocity component and possibly a non-zero azimuthal

velocity component. Nevertheless, our domain of interest is a disk sector defined by −4◦ ≤ 𝜃 ≤ 4◦,

for which the changes in depth are very small, making the azimuthal changes in pressure and

velocity (and hence the azimuthal component of velocity) negligible.

In this region, the radial

component of velocity dominates.

84

The flow needs to transition from purely axial to radial near the connection between the inlet

pipe and the pad. Near this region, the flow between the pad and wall still has an axial component,

but as the flow progresses through the gap, the velocity field gradually becomes purely radial:

u(𝑟, 𝑧) = 𝑢𝑟 (𝑟, 𝑧) ˆr + 0 ˆθ + 0ˆz. Here, we emphasize that the velocity field will asymptotically become

axisymmetric and a function of only the radial and axial coordinates. This is due to the fact that

the circumferential cross-section increases radially and that the two parallel plates impose non-slip

boundary conditions. This is the flow condition that we will use for our analysis.

Scaling

The following are appropriate dimensionless versions of each physical quantity:

𝑟 ∗ =

𝜕𝑟 ∗
𝜕𝑟

=

𝑟
𝑑

1
𝑑

𝑧∗ =

𝑧
ℎeq

𝜕𝑧∗
𝜕𝑧

=

1
ℎeq

𝑢∗ =

𝑢𝑟
𝑢in

𝑝∗ =

𝑝
𝜌𝑢2
in

Conservation Equations

Continuity

The general form of the continuity equation in cylindrical coordinates is:

85

𝜕 (𝑟𝑢𝑟 )
𝜕𝑟

1
𝑟

𝜕𝑢 𝜃
𝜕𝜃

1
𝑟

+

𝜕𝑢 𝑧
𝜕𝑧

+

= 0

For purely radial axisymmetric flow, we assume the conditions 𝑢 𝜃 = 𝑢 𝑧 = 0. The continuity

equation in cylindrical coordinates for incompressible-steady flow can then be simplified and

expanded as follows:

In dimensionless form:

𝑟

𝜕𝑢𝑟
𝜕𝑟

+ 𝑢𝑟 = 0

𝑟 ∗𝑢in𝑑

𝜕𝑢∗
𝜕𝑟 ∗ ·

𝜕𝑟 ∗
𝜕𝑟

+ 𝑢in𝑢∗ = 0

𝑟 ∗𝑢in𝑑

𝜕𝑢∗
𝜕𝑟 ∗ ·

1
𝑑

+ 𝑢in𝑢∗ = 0

𝑟 ∗ 𝜕𝑢∗

𝜕𝑟 ∗ + 𝑢∗ = 0

(A.3)

Then, partially integrating equation (A.3) gives the following result for the velocity field:

𝑢∗

∫

𝑢∗
𝑜 (𝜃,𝑧∗)

𝜕𝑢∗
𝑢∗ = −

𝑟 ∗

∫

𝑟 ∗
𝑜

𝜕𝑟 ∗
𝑟 ∗

(cid:19)

(cid:18)

𝑢∗
𝑜 (𝜃, 𝑧∗)
𝑢∗

ln

= − ln

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

𝑢∗(𝑟 ∗, 𝜃, 𝑧∗) =

𝑟 ∗
𝑜
𝑟 ∗

𝑜 (𝜃, 𝑧∗)
𝑢∗

(A.4)

Which shows the variation of velocity non-dimensionally.

86

Dropping the dependence on 𝜃 and writing the result in vector form, we obtain:

u∗(𝑟 ∗, 𝑧∗) =

𝑟 ∗
𝑜
𝑟 ∗

𝑜 (𝑧∗) ˆr
𝑢∗

Where 𝑢∗

𝑜 (𝑧∗) is the radial component of velocity at the reference radius, 𝑟 ∗

𝑜. The velocity field

𝑜 (𝑧∗) is a function of 𝑧∗ resulting from the non-slip boundary conditions and viscous effects.
𝑢∗

Momentum Equations

The next step in determining the asymptotic behavior for the pressure field is to use Navier-

Stokes equations for incompressible-steady flow considering the same simplifications imposed on

the continuity equation (𝑢𝜃 = 𝑢𝑧 = 0), plus one simplification related to the azimuthal derivative:

𝜕𝑢𝑟 /𝜕𝜃 = 0. This last simplification is imposed due to the narrow range of 𝜃 used in our simulations.

r-momentum equation:

After simplifications, the 𝑟-momentum equation becomes:

𝑢𝑟

𝜕𝑢𝑟
𝜕𝑟

= −

1
𝜌

𝜕 𝑝
𝜕𝑟

+ 𝑔𝑟 + 𝜈

(cid:20)

−

1
𝑟 2

(cid:18)

𝑟

𝜕𝑢𝑟
𝜕𝑟

+ 𝑢𝑟

(cid:19)

+

(cid:21)

𝜕2𝑢𝑟
𝜕𝑧2

(A.5)

Here, 𝑔𝑟 is the component of gravity in the radial direction as defined in equation (A.1).

Applying the continuity equation (A.3) in (A.5):

𝑢2
𝑟
𝑟

−

= −

1
𝜌

𝜕 𝑝
𝜕𝑟

+ 𝑔𝑟 + 𝜈

𝜕2𝑢𝑟
𝜕𝑧2

In dimensionless form:

(𝑢∗𝑢in)2
𝑟 ∗𝑑

−

(cid:1)

𝜕 (cid:0)𝑝∗𝜌𝑢2
in
𝜕𝑟 ∗

1
𝜌

𝜕𝑟 ∗
𝜕𝑟

·

= −

+ 𝑔𝑟 + 𝜈

𝜕2 (𝑢∗𝑢in)
𝜕𝑧∗𝜕𝑧∗

𝜕𝑧∗
𝜕𝑧

𝜕𝑧∗
𝜕𝑧

·

·

𝑢2
in
𝑑

−

𝑢∗2
𝑟 ∗ = −

𝑢2
in
𝑑

·

𝜕 𝑝∗
𝜕𝑟 ∗ + 𝑔𝑟 +

·

𝜈𝑢in
ℎ2
eq

𝜕2𝑢∗
𝜕𝑧∗2

·

87

𝜕 𝑝∗
𝜕𝑟 ∗ =

𝑔𝑟 𝑑
𝑢2
in

+

𝑢∗2
𝑟 ∗ +

𝜈𝑢in
ℎ2
eq

·

𝑑
𝑢2
in

𝑑
𝑑

·

𝜕2𝑢∗
𝜕𝑧∗2

·

𝜕 𝑝∗
𝜕𝑟 ∗ = − [sin 𝛼 sin(𝜃 + 𝜙)]

𝑔𝑑
𝑢2
in

+

𝑢∗2
𝑟 ∗ +

𝜈𝑢in
ℎ2
eq

·

𝑑
𝑢2
in

𝑑
𝑑

·

𝜕2𝑢∗
𝜕𝑧∗2

·

𝜕 𝑝∗
𝜕𝑟 ∗ = − [sin 𝛼 sin(𝜃 + 𝜙)] Riin +

𝑢∗2
𝑟 ∗ +

1
ℎ∗2Rein

𝜕2𝑢∗
𝜕𝑧∗2

·

(A.6)

Here, Riin = 𝑔𝑑/𝑢2

in is the Richardson number based on inlet conditions, Rein = 𝑢in𝑑/𝜈 is the
Reynolds number based on inlet conditions, and ℎ∗ = ℎeq/𝑑 is the dimensionless equilibrium gap.

Equation (A.6) suggests that viscous shear is scaled by the product ℎ∗2Rein.

We can now substitute 𝑢∗ from equation (A.4) into equation (A.6) (dropping 𝜃 dependence from

𝑢∗ to conform to our experimental conditions):

𝜕 𝑝∗
𝜕𝑟 ∗ = − [sin 𝛼 sin(𝜃 + 𝜙)] Riin +

𝑟 ∗2
𝑜
𝑟 ∗3

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2 +

1
ℎ∗2Rein

𝑟 ∗
𝑜
𝑟 ∗ ·

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

This equation can be integrated partially:

𝑝∗

∫

𝑜 (𝜃,𝑧∗)
𝑝∗

𝜕 𝑝∗ = − [sin 𝛼 sin(𝜃 + 𝜙)] Riin

𝑟 ∗

∫

𝑟 ∗
𝑜

𝜕𝑟 ∗ + 𝑟 ∗2
𝑜

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

𝑟∗

∫

𝑟 ∗
𝑜

𝜕𝑟 ∗
𝑟 ∗3

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

𝑝∗(𝑟 ∗, 𝜃, 𝑧∗) = − [sin 𝛼 sin(𝜃 + 𝜙)] Riin

(cid:0)𝑟 ∗ − 𝑟 ∗
𝑜

(cid:1) +

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

1
2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

+ 𝑝∗

𝑜 (𝜃, 𝑧∗)

We consider the case 𝛼 = 𝜋/2. Hence, equation (A.7) can be re-written as:

𝑝∗(𝑟 ∗, 𝜃, 𝑧∗) = − sin(𝜃 + 𝜙)Riin

(cid:0)𝑟 ∗ − 𝑟 ∗
𝑜

(cid:1) +

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

88

1
2
(cid:19)

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

+ 𝑝∗

𝑜 (𝜃, 𝑧∗)

𝑟∗

∫

𝑟 ∗
𝑜

𝜕𝑟 ∗
𝑟 ∗

(A.7)

(A.8)

𝜃-momentum equation:

After simplifications, the 𝜃-momentum equation becomes:

1
𝜌𝑟

𝜕 𝑝
𝜕𝜃

−

+ 𝑔𝜃 = 0

(A.9)

In dimensionless form:

𝜌𝑢2
in
𝜌𝑟 ∗𝑑

𝜕 𝑝∗
𝜕𝜃

= 𝑔𝜃

𝜕 𝑝∗
𝜕𝜃

=

𝑟 ∗

𝑔𝜃 𝑑
𝑢2
in

𝜕 𝑝∗
𝜕𝜃

𝜕 𝑝∗
𝜕𝜃

= − [sin 𝛼 cos(𝜃 + 𝜙)]

𝑟 ∗

𝑔𝑑
𝜌𝑢2
in

= − [sin 𝛼 cos(𝜃 + 𝜙)] Riin𝑟 ∗

(A.10)

Equation (A.10) can now be integrated partially:

𝑝∗

∫

𝑜 (𝑟 ∗,𝑧∗)
𝑝∗

𝜕 𝑝∗ = −Riin𝑟 ∗ sin 𝛼

𝜃

∫

𝜃𝑜

cos(𝜃 + 𝜙)𝜕𝜃

𝑝∗(𝑟 ∗, 𝜃, 𝑧∗) = −Riin𝑟 ∗ sin 𝛼 [sin(𝜃 + 𝜙) − sin(𝜃𝑜 + 𝜙)] + 𝑝∗

𝑜 (𝑟 ∗, 𝑧∗)

(A.11)

We consider the case for 𝛼 = 𝜋/2. Equation (A.11) can then be rewritten as:

𝑝∗(𝑟 ∗, 𝜃, 𝑧∗) = −Riin𝑟 ∗ [sin(𝜃 + 𝜙) − sin(𝜃𝑜 + 𝜙)] + 𝑝∗

𝑜 (𝑟 ∗, 𝑧∗)

(A.12)

In similar way we can obtain a relation with 𝑧-momentum equation. This will help in realizing the

flow field in all the three coordinates.

89

𝑧-momentum equation:

After simplifications, the 𝑧-momentum equation becomes:

In dimensionless form:

𝜕 𝑝
𝜕𝑧

= 𝜌𝑔𝑧

𝜌𝑢2
in

𝜕 𝑝∗
𝜕𝑧∗ ·

𝜕𝑧∗
𝜕𝑧

= 𝜌𝑔𝑧

𝜌𝑢2
in

𝜕 𝑝∗
𝜕𝑧∗ ·

1
ℎeq

= 𝜌𝑔𝑧

𝜕 𝑝∗
𝜕𝑧∗ =

𝜌 (−𝑔 cos 𝛼)
𝜌𝑢2
in

· ℎeq ·

𝑑
𝑑

𝜕 𝑝∗
𝜕𝑧∗ = −(cos 𝛼)Riinℎ∗

(A.13)

Equation (A.13) can now be integrated partially:

𝑝∗

∫

𝑜 (𝑟 ∗,𝜃)
𝑝∗

𝜕 𝑝∗ = −(cos 𝛼)Riinℎ∗

𝑧∗

∫

𝑧∗
𝑜

𝜕𝑧∗

𝑝∗(𝑟 ∗, 𝜃, 𝑧∗) = −(cos 𝛼)Riinℎ∗ (cid:0)𝑧∗ − 𝑧∗
𝑜

(cid:1) + 𝑝∗

𝑜 (𝑟 ∗, 𝜃)

(A.14)

Considering the case 𝛼 = 𝜋/2, we obtain:

𝑝∗(𝑟 ∗, 𝜃, 𝑧∗) = 𝑝∗

𝑜 (𝑟 ∗, 𝜃)

(A.15)

Equation (A.15) imposes the condition that 𝑝∗ must be independent of 𝑧∗ when 𝛼 = 𝜋/2. With

this in mind, equations (A.8), (A.12), and (A.15), can be re-written.

90

The equation obtained will remove the 𝑧 coordinate dependency.

𝑝∗(𝑟 ∗, 𝜃) = − sin(𝜃 + 𝜙)Riin

(cid:0)𝑟 ∗ − 𝑟 ∗
𝑜

(cid:1) +

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

1
2
(cid:19)

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

+ 𝑝∗

𝑜 (𝜃)

𝑝∗(𝑟 ∗, 𝜃) = −Riin𝑟 ∗ [sin(𝜃 + 𝜙) − sin(𝜃𝑜 + 𝜙)] + 𝑝∗

𝑜 (𝑟 ∗)

𝑝∗(𝑟 ∗, 𝜃) = 𝑝∗

𝑜 (𝑟 ∗, 𝜃)

(A.16)

(A.17)

(A.18)

Since all these equations must in principle describe the same pressure field, the results must be

compatible and satisfy all the differential equations for conservation of linear momentum. To find

the final form of the pressure field, let’s first subtract equation (A.17) from (A.16):

0 = − sin(𝜃 + 𝜙)Riin

(cid:0)𝑟 ∗ − 𝑟 ∗
𝑜

(cid:1) +

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

1
2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

+ 𝑝∗

𝑜 (𝜃) + Riin𝑟 ∗ [sin(𝜃 + 𝜙) − sin(𝜃𝑜 + 𝜙)] − 𝑝∗

𝑜 (𝑟 ∗)

0 = Riin𝑟 ∗

𝑜 sin(𝜃 + 𝜙) − Riin𝑟 ∗ sin(𝜃𝑜 + 𝜙) +

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

1
2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

+ 𝑝∗

𝑜 (𝜃) − 𝑝∗

𝑜 (𝑟 ∗)

(A.19)

For equation (A.19) to be satisfied, all the terms that are functions of 𝜃 must cancel each other,

and all the terms that are functions of 𝑟 ∗ must cancel each other as well:

0 = Riin𝑟 ∗

𝑜 sin(𝜃 + 𝜙) + 𝑝∗

𝑜 (𝜃)

⇒ 𝑝∗

𝑜 (𝜃) = −Riin𝑟 ∗

𝑜 sin(𝜃 + 𝜙)

(A.20)

91

0 = −Riin𝑟 ∗ sin(𝜃𝑜 + 𝜙) +

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

1
2

ln

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

− 𝑝∗

𝑜 (𝑟 ∗)

Further simplifying the above equation we will get,

⇒ 𝑝∗

𝑜 (𝑟 ∗) = −Riin𝑟 ∗ sin(𝜃𝑜 + 𝜙) +

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

1
2

ln

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

(A.21)

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

We now add equations (A.16) and (A.17):

2𝑝∗(𝑟 ∗, 𝜃) = − sin(𝜃 + 𝜙)Riin

(cid:0)𝑟 ∗ − 𝑟 ∗
𝑜

(cid:1) +

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

1
2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

− Riin𝑟 ∗ [sin(𝜃 + 𝜙) − sin(𝜃𝑜 + 𝜙)] + 𝑝∗

𝑜 (𝜃) + 𝑝∗

𝑜 (𝑟 ∗)

Substitute 𝑝∗

𝑜 (𝜃) and 𝑝∗

𝑜 (𝑟 ∗) from the results (A.20) and (A.21) respectively:

2𝑝∗(𝑟 ∗, 𝜃) = − sin(𝜃 + 𝜙)Riin

(cid:0)𝑟 ∗ − 𝑟 ∗
𝑜

(cid:1) +

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

1
2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

− Riin𝑟 ∗ [sin(𝜃 + 𝜙) − sin(𝜃𝑜 + 𝜙)] − Riin𝑟 ∗
𝑟 ∗2
𝑜
𝑟 ∗2

− Riin𝑟 ∗ sin(𝜃𝑜 + 𝜙) +

𝑜 (𝑧∗)(cid:3) 2

(cid:2)𝑢∗

1 −

(cid:18)

1
2

𝑜 sin(𝜃 + 𝜙)
(cid:19)

𝑟 ∗
𝑜
ℎ∗2Rein

+

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

⇒ 𝑝∗(𝑟 ∗, 𝜃) = −Riin𝑟 ∗ sin(𝜃 + 𝜙) +

(cid:2)𝑢∗

𝑜 (𝑧∗)(cid:3) 2

1
2

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

+

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

(A.22)

Finally, for condition (A.15) to be satisfied, equation (A.22) must not be a function of 𝑧∗. Which

implies:

𝜕 𝑝∗(𝑟 ∗, 𝜃)
𝜕𝑧∗

= 0

92

(A.23)

Figure A.1 Scatter plot of ˜(cid:164)𝑊 vs 𝐷∗ exhibiting no correlation.

Taking the partial derivative of equation (A.22) and making it equal to zero:

0 =

(cid:18)

1 −

(cid:19)

𝑟 ∗2
𝑜
𝑟 ∗2

𝑢∗
𝑜 (𝑧∗)

𝜕𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗

+

𝑟 ∗
𝑜
ℎ∗2Rein

· ln

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

(cid:19) 𝜕3𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗3

(A.24)

The solution to equation (A.24) guarantees that condition (A.8) is met, and gives the reference

velocity profile inside the Bernoulli pad’s gap.

Observations Related to Scaling Analysis

The dimensional analysis suggests a function of the form:

˜(cid:164)W = Φ (𝐷∗, ℎ∗, ˜𝜏, Rein)

(A.25)

However, as shown in Figure A.1, there is no correlation between ˜(cid:164)W and 𝐷∗, simplifying the

dimensionless function to:

˜(cid:164)W = Φ (ℎ∗, ˜𝜏, Rein)

(A.26)

On the other hand, the third term on the right-hand-side of equation (A.22) represents the effect

of viscous resistance on the pressure field. This term can be used to get an insight on how the

parameters from equation (A.26) relate to each other.

93

The equation now becomes:

𝑝∗(𝜃, 𝑧∗) = · · · + · · · +

𝑟 ∗
𝑜
ℎ∗2Rein

𝜕2𝑢∗

𝑜 (𝑧∗)
𝜕𝑧∗2

·

ln

(cid:19)

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

𝑝∗(𝜃, 𝑧∗) = · · · + · · · + 𝑟 ∗

𝑜 ln

𝑝∗(𝜃, 𝑧∗) = · · · + · · · + 𝑟 ∗

𝑜 ln

𝑝∗(𝜃, 𝑧∗) = · · · + · · · + 𝑟 ∗

𝑜 ln

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

(cid:19) 𝜕
𝜕𝑧∗

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

(cid:19) 𝜕
𝜕𝑧∗

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

(cid:19) 𝜕
𝜕𝑧∗

(cid:20)

(cid:20)

(cid:20)

1
ℎ∗2Rein

1
ℎ∗2Rein

1
ℎ∗2Rein

·

·

·

𝜇
𝜇

1
𝜇

·

·

ℎeq
𝑢in

ℎeq
𝑢in

·

·

(cid:21)

𝜕𝑢𝑜 (𝑧)
𝜕𝑧

(cid:21)

· 𝜏𝑜 (𝑧)

𝑑
𝑑

ℎeq
𝑑

·

(cid:21)

𝜏𝑜 (𝑧) 𝑑
𝜇𝑢in

𝑝∗(𝜃, 𝑧∗) = · · · + · · · + 𝑟 ∗

𝑜 ln

(cid:18) 𝑟 ∗
𝑟 ∗
𝑜

(cid:19) 𝜕
𝜕𝑧∗

(cid:21)

(cid:20) 𝜏∗
𝑜 (𝑧∗)
ℎ∗Rein

(A.27)

Equation (A.27) suggests that the local viscous shear, 𝜏∗

𝑜 (𝑧∗), scales with ℎ∗Rein. This relation

suggest that the maximum shear, 𝜏max, might also be scaled with ℎ∗Rein. With this new insight in

mind, the dimensionless function is even simpler:

˜(cid:164)W = Φ ( ˜𝜏, ℎ∗Rein)

(A.28)

94