---
# System prepended metadata

title: 3.2 Disciplines in Hydrodynamics
tags: [Mathematical modeling and optimization]

---

---
tags: Mathematical modeling and optimization
title: 3.2 Disciplines in Hydrodynamics
---

## 3.2 Disciplines in Hydrodynamics

Hydrodynamic phenomena can all be traced back and described by the Navier-Stokes Equation. For incompressible ($\rho$ = const) flow under [Newtonian fluid](https://en.wikipedia.org/wiki/Newtonian_fluid) regime, the equation set yields:

$$
\begin{align}
\frac{\partial u_i}{\partial x_i} &= 0 \\
\\
\frac{\partial u_i}{\partial t} + u_j \frac{\partial u_i}{\partial x_j} &= \frac{1}{\rho}\frac{\partial p}{\partial x_i} + \nu\frac{\partial ^2 u_i}{\partial x_j} + g_i
\end{align}
$$

wehre $\nu$ stands for kinematic viscosity, which is a individual material characteristics. 

All the disciplines in the hydrodynamic framework can be derived from the equation above. In this section, we will focus on the certain engineering applications using the simplified behavior.  


### 3.2.1 Laminar pipe flow and Hagen-Poiseuille equation


Laminar flow is the classical flow type where both the inertial and viscous forces are evently matched. Due to its stability, analytical solutions can be achieved in some simplified circumstances.

For the pipe flow, the Navier-Stokes Equation shall be formulated the [cylinder coordinate system](https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Cylindrical_coordinates) to achieve plausible analysis. Based on this, a fully developed laminar pipe flows have some special feastures which simplify the governing equation. For the streamwise direction, the flow is assumed to be homogeneous, denoting $\displaystyle \frac{\partial u_x}{\partial x_x} = 0$. On the other hand, the axisymmetric feature of pipe flows restricts the streamwise velocity profile to be only dependent of radial direction $r$. In accompany with the steady state condition $\displaystyle \frac{\partial}{\partial t} = 0$, the momentum equation for streamwise direction ($x$) in the cylindrical coordinate system is simplified to
  
$$
\displaystyle \frac{1}{\mu}\frac{\partial p}{\partial x} =\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_x}{\partial r}\right)
$$

which has the analytical solution of velocity as

$$
\displaystyle u_x(r) = -\frac{1}{4 \mu}\frac{\partial p}{\partial x}\ \left(R^2 - r^2 \right)
$$

where $R$ is the radius of the pipe, $\displaystyle \frac{\partial p}{\partial x}$ is the pressure gradient in the streamwise direction, $\mu$ is the dynamic viscosity and $r$ is the dimension in the axial direction. 

Under above mentioned assumptions, this equation reveals the velocity at any point in the pipe with a given pressure drop. Upon this, volumetric flow rate is calculated by integrating the velocity profile over cross section. Also, thinking on the pressure gradient as the pressure drop ($\Delta p$) within the pipe length $L$. The volumetric flow rate yields to

$$
\displaystyle \dot{Q} \ = \ \frac{\pi}{8\mu} \frac{\partial p}{\partial x}R^4 \ =\  \frac{\pi}{8\mu} \frac{\Delta p}{L}R^4
$$

This is called [Hagen-Poiseuille Equation](https://en.wikipedia.org/wiki/Hagen%E2%80%93Poiseuille_equation), stating the __flow rate__ yielded by the __defined pressure difference__ in the laminar and drag free condition.

Reformulating the equation, the pressure difference between the two ends ($P_s$ (source), $P_d$ (destination)) can be written in the form:
$$
\Delta p_{\text{steady}} = \frac{8 \mu \dot{Q} L}{\pi R^4}
$$


Extending the scope to the dynamic condition using Newton’s 2nd law to describe the [inertial force](https://en.wikipedia.org/wiki/Inertance), i.e. net force = mass × acceleration, 
$$
\begin{align}
\Delta p_{\text{Inertance}} \cdot A=& \ m\frac{du}{dt}= \rho A L \frac{1}{A}\frac{dQ}{dt}=\rho L\frac{dQ}{dt}\\
\\
\Delta p_{\text{Inertance}} =& \ \frac{\rho L}{A}\frac{d\dot{Q}}{dt}\\
\end{align}
$$
And the total pressure drop in a trasient system obtains both the steady and unsteady parts. Thus,

$$
P_s - P_d  = \Delta p_{\text{Inertance}} + \Delta p_{\text{Steady}} = \ \frac{\rho L}{A}\frac{d\dot{Q}}{dt} + \ \frac{8 \mu  L}{\pi R^4}\dot{Q}
$$

#### Hydraulic Actuator 
Hydraulic actuators are devices that convert the energy stored in pressurized fluid into controlled mechanical motion. In their simplest form, they consist of a piston or rotary vane enclosed in a cylinder, connected to a hydraulic supply and return line. When fluid is forced into one side of the actuator chamber, the resulting pressure produces a force on the piston proportional to the pressure and piston area, thereby generating linear displacement or torque.

<img src="https://mechanical-engineering.com/wp-content/uploads/2016/07/hydraulic-actuator-linear-actuator-420x228.jpg" width="46%">
<img src="https://hydrauliccylindermanufacturers.net/wp-content/uploads/2022/11/linear-actuators.jpg" width="46%">

Since the actuator is not fed by vicous fluid, the connecting hoses and valves can often be described by the Hagen–Poiseuille for flow within these narrow passages. This means the flow rate depends on pressure drop, fluid viscosity, and line dimensions. These restrictions create pressure losses and time delays that directly affect actuator response.

__<ins>Scenario: spring-return landing gear actuator:__

Single-acting spring return actuators are devices which are often implemeted in valve control systems, industrial machinery, and construction equipment, where reliable, repeatable motion is of central request. They convert hydraulic pressure into straight-line motion. When fluid pushes the piston out, the spring pulls it back once pressure is released.

<img src="https://www.researchgate.net/publication/335787001/figure/fig2/AS:802490052255745@1568339946849/Linear-hydraulic-actuator.png" width="65%">

<br/>
<br/>
    
    
Give parameters are:  
- Gear mass($M$) $200 \  kg$
- Stroke length $0.5\  m$
- Spring stiffness ($K$) $1.5×10^5\ N/m$
- Cylinder bore (d) 0.1 m, leads to A=$7.85×10^{−3} m^2$
- Hydraulic line: radius ($r$) 10 mm, Length ($L$) $0.3 m$
- Damping coefficient of the piston ($B$) $5000 N\cdot s/m$
- Input pressure $P_{in}(t)$ gradient to 15 MPa (typical hydraulic supply) within 0.01 second.
- Dynamic viscosity of the fluid ($\mu$) $1.0x10^{-3} \ Pa\cdot s$

_How is the piston displacement and the corresponding pressure inthe actuator?_
    
[**<ins>Sample solution**](https://cocalc.com/share/public_paths/64305df0acbc3f480732ad7d3c0c6c178c349e9f)
    
#### Frictional contribution

Viscous flow is always friction related. Friction effects leads to the loss of pressure and energy during the flow proceeding. To calculate the friction contribution on the laminar pipe flow, we can reformulate the flow rate $\dot{Q}$ to the form of bulk velocity $\bar{U}$. Thus,
$$
\bar{U} = \frac{\dot{Q}}{A} = \frac{R^2}{8\mu} \frac{\Delta p}{L}.
$$
Following the formulation of [pressure loss coefficient](https://de.wikipedia.org/wiki/Rohrreibungszahl) and setting $D = 2R$ as the diameter of the pipe, the drag coefficient shows:
$$
\begin{align}
C_d &= \frac{F_d}{\frac{1}{2}\rho \bar{U}^2 A}=  \frac{\Delta p }{\frac{1}{2}\rho \bar{U}^2} \\
\\
&= 64\frac{\mu}{\rho \ D\ \bar{U}}\frac{L}{D} = \frac{64}{\text{Re}}\frac{L}{D} 
\end{align}
$$

Following the concept of [Darcy-Weisbach Equation](https://en.wikipedia.org/wiki/Darcy%E2%80%93Weisbach_equation), the friction coefficient yields to $\displaystyle f_D = \frac{64}{\text{Re}}$, which will be the laminar part of the Moody's diagram later.

#### Applications in the medicine technology

The concept of Hagen-Poiseuille flow are widely applied especially in the microfluidic and medicine technology, where the Reyolds number stays low (laminar) and pressure drop as well as cross section area as the main driven effects for the flow patterns.
<br/>

<img style="float: right;" src="https://my.clevelandclinic.org/-/scassets/images/org/health/articles/21697-vasoconstriction-illustration" width=25%>

An example in the hemodynamics is the vein constriction relates to the high blood pressure. Based on the Hagen-Poiseuille equation, the vein contriction leads to descreasing of blood flow rate (in the power of 4). On the other hand, the heart will try to keep the same flow rate to assure the body working intact. So it pumps with higher pressure accordingly. This leads to vein or artery rapture, worst case as a stroke.   

<br/>
<br/>

<img style="float: right;"  src="https://upload.wikimedia.org/wikipedia/commons/thumb/f/fa/Mantoux_tuberculin_skin_test.jpg/220px-Mantoux_tuberculin_skin_test.jpg">

Injection of medical substance can be another application of Hagen-Poiseuille equation. The flow rate on the needle is determined by the the blood pressure in the vein, the dynamic visocity of medicine substance, the diamter of needle and pressure of substance reservoir, i.e. the pressure medical staff creates by pressing the syringe rod. Here, the concept of [pascal's law](https://en.wikipedia.org/wiki/Pascal%27s_law#:~:text=Pascal's%20principle%20is%20defined%20as,angles%20to%20the%20enclosing%20walls.) is applied again.
    
#### Temperature impact in fluid viscosity

Viscosity is a substitute for shear force of fluid. It can generally recognized as a material property, stating as the slope to velocity gradient with given shear $\tau$, such $\displaystyle \tau = \mu \frac{\partial u }{\partial y}$.
The viscosity doesn't have to be a constant value. In several material type, the viscosity can be dependent on pressure, shear, temperature, etc. The study of viscosity effect is in the discipline of [Rheology](https://en.wikipedia.org/wiki/Rheology).
Temperature is one of the impact for fluid viscosity. Distinguish by its states, different models are proposed the describe the behavior. For __gas phase__ medium, [Sutherland model](https://en.wikipedia.org/wiki/Temperature_dependence_of_viscosity#Sutherland_model) is widely used with reference viscosity and temperature ($\mu_{\text{ref}}$, $T_{\text{ref}}$) and a modeling parameter $S$ for the offset.
$$
\mu = \mu_{\text{ref}} \left( \frac{T}{T_{\text{ref}}}\right)^{3/2}\frac{T_{\text{ref}}+S}{T + S}
$$
For __liquid phase__, the description is more rudimentary. An [Arrehnius form](https://en.wikipedia.org/wiki/Arrhenius_equation) showing exponential behavior to temperature is used to applied for non-constent liquid viscoisty, 
$$
\mu = A\exp(\frac{B}{T})
$$
Some common materials with their modeling parameters can be found [here](https://en.wikipedia.org/wiki/Temperature_dependence_of_viscosity#Liquids).

__<ins>Scenario -- Syringe pump infusion through a microcatheter__

[Syring pump](https://en.wikipedia.org/wiki/Syringe_driver) is a medicine device specifically to deliver small amount of dose. Incorporated with microcatheter ($d \lt 0.5$ mm), drug delivery can be precisely conducted to body areas, wehere it is difficult to reach. 

<img src="https://www.terumo-europe.com/sites/productcatalog/Lists/ProductImages/terufusion-syringe-pump-smart-product-image.png" width="60%"><img src="https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcSPWf25slh1KQ8gan7YNXReV8OA7LmEgKqDrA&s" width="40%">
    
Construction and function of a syringe pump is straigtforward. A syringe is attached into a stepper motor, which presses the syringe piston to inject the fluid with defined displacement. The flow rate and dose can be precisely determined.

The flow type in the mircocatheter follows the Hagen-Poiseuille patterns, where the fluid viscosity determines the resistance. Besides, the liquid drug can be temperature depenedent. However, unexpected event might occurs such as tube kinking or blocking (which is called [__occlusion__](https://en.wikipedia.org/wiki/Vascular_occlusion) in the medicine vocabulary), such effects will cause sudden increase of flow resistance. i.e. danger in the system. 

 
    
_Say a syringe pump is delivering a liquid drug at a prescribed flow rate of 10 µL/min through a long, narrow microcatheter into a patient. Initially, the drug solution is at room temperature (20 °C) and relatively viscous, which increases the flow resistance described by the Hagen–Poiseuille law. As the infusion continues, the syringe and catheter warm toward body temperature (37 °C).
Midway through the infusion, at 150 seconds, the catheter becomes partially occluded, which increases the resistance by a factor of 10. This causes a rapid rise in syringe pressure and a collapse of flow._ 

_What is the viscosity and the pressure in the period?_    

[**<ins>Sample solution**](https://cocalc.com/share/public_paths/1594583eacd0133e50cf0bdad7791ba5c58acf0e)
    
---

### 3.2.2 Bernoulli equation and the Venturi flow types 

When the viscous effects are negligible for the focusing radar, the Navier-Stokes Equations are simplified to [Euler Equations](https://en.wikipedia.org/wiki/Euler_equations_(fluid_dynamics)). In a steady and irrotational state, the Euler equation simplifies to:
$$
u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{ρ}\frac{\partial p}{\partial x_i}+g_i
$$
Now, thinking on integrating the Euler equation along a piece of stream line $s$ from $A$ to $B$, taking $v$ as the velocity component along the stream line. Thus,
$$
\begin{align}
\int_A^B v\frac{\partial v}{\partial s} ds &= \int_A^B -\frac{1}{ρ}\frac{\partial p}{\partial s} ds + \int_A^B g_i\  ds \\
\\
\text{with a further step of}\qquad\qquad\qquad &&\\
\int_A^B d(\frac{1}{2} v^2) &= -\frac{1}{\rho} \int_A^B dp + g \int_A^B dh\\
\text{the integration results to}\qquad\qquad\qquad &\\
\\
\frac{1}{2} \rho v_A^2  + p_A +& \rho g h_A = \frac{1}{2}\rho v_B^2  + p_B + \rho g h_B
\\
\end{align}
$$
This equation is known as the [Bernoulli equation](https://en.wikipedia.org/wiki/Bernoulli%27s_principle). It describes the conservation of mechanical energy along a streamline in steady, irrotational flow of an inviscid fluid, neglecting frictional losses. The terms represent the kinetic energy per unit mass, pressure energy per unit mass, and potential energy per unit mass respectively.

In another words, for a given pressure difference on a streaming system, the velocity as well as flow rate can be evaluate using the concept of Bernoulli. Differ to Hagen-Poiseuille flow, the application of Bernoulli equation is not restricted in the laminar flow state. For all the Euler equation classified flow types (viscous effects negligible), the Bernoulli equation is applicable.

<img style="float: right;" src="https://dam.krohne.com/image/upload/t_ar43_pd_c/c_scale,w_540/q_auto/dpr_2.625/d_im-other:image-not-available.png/f_auto/v1/im-single-product/optibar/venturi-flowmeter-manifolds-pressure-transmitter-temperature-sensor?_a=ATAPpAA0" width=60%>

<br/>

[__Venturi fluid flow meter__](https://en.wikipedia.org/wiki/Venturi_effect) is a classical application of Bernoulli equation. The Venturi flow yield a speed up of flow velocity and decreasing of static pressure in the convergent section (section A), as well as slow down and pressure increment in the divengent section (section B). 

Measuring the pressure difference ($p_B - p_A$) from both side. The flow rate $Q$ yields:
$$
\displaystyle Q = A_1 \sqrt{\frac{2}{\rho}\frac{p_B - p_A}{(\frac{A_A}{A_B})^2 - 1}}
$$
where $A_A$, $A_B$ are the cross section area of the measuring position.

<br/>

<img style="float: right;" src="https://5.imimg.com/data5/SELLER/Default/2023/1/HR/YB/FP/120110231/2-pcs-nozzle-1-500x500.jpg" width=35%>

When __water flows through a nozzle__ device, such as a nozzle attached to a hose, its velocity can increase due to the constriction of the flow area. The Bernoulli equation can be applied to analyze the changes in pressure and velocity as the water flows through the nozzle. Setting point 1 in the reservoir and point 2 on the nozzle, the behavior from point 1 to point 2 reads
$$
\displaystyle gz_1+\frac{p_1}{\rho} =\frac{1}{2}v_2^2 +gz_2​+\frac{p_2}{\rho}
$$
velocity at nozzle yields
$$
\displaystyle  v_2 = \sqrt{2 \left(\frac{p_1 - p_2}{\rho} + g(z_1 -z_2)\right)}
$$

#### Head loss and frictional effects    
Although the Bernoulli equation is dervied from the momentum equation, the terms can rather be interpreted in the energy representation. In the hydraulic context, it can be rewritten in the __head (height)__ representation as
$$
H_{\text{total}} = \underbrace{\frac{v_A^2}{2g}}_{\text{Velocity Head}}  + \underbrace{\frac{p_A}{\rho g}}_{\text{Presure Head}} + \underbrace{h_A}_{\text{Elevation Head}} = \frac{v_B^2}{2g}  + \frac{p_B}{\rho g} + h_B 
$$

For systems with no energy loss, i.e. invicid flow, the conservation of total head is usually taken for analysis. The physical properties at two different points along the same streamline can transfer from pressure to velocity head or elevation head, or vice versa.
    
However, invicid flow doesn't exist in the reality. For every flow movement, viscous effect accompanies. Incorporated to the wall, vicous force becomes friction for the movement. On the Bernoulli equation basis, the energy loss due to frictional effects is declared as head loss. Thus, following the stream line ($s$) from point $A$ to point $B$, the Bernoulli equation is modified to
$$
    \frac{v_A^2}{2g}  + \frac{p_A}{\rho g}+ h_A = \frac{v_B^2}{2g}  + \frac{p_B}{\rho g} + h_B + h_{\text{loss}}
$$
with
$$
    h_{\text{loss}}  = \int_A^B \underbrace{\left(\iint_{\partial A}\mu \frac{\partial u}{\partial y} \ dx dz\right)}_{\text{frictional force}} \ ds   
$$

For pipe flow, [Darcy Weisbach Equation](https://en.wikipedia.org/wiki/Darcy%E2%80%93Weisbach_equation) provides empirical formulations incorporated with a frictional coefficient to adapt to different flow types. 

   
#### __<ins>Scenario: Pipe with diameter adaption__
    
<img src="https://live.staticflickr.com/65535/54846906237_b8b4b90fc9_b.jpg">

_A section in the pipe system transits from DN150 to DN100. In position $A$, the pressure measurement reads 1.7 bar(g) for transporting water with flow rate of 400 $m^3/h$. Neglecting the effects on the diamater adaptor, what is the expected pressure on position $B$?_ 

[__<ins>Sample solution__](https://cocalc.com/share/public_paths/72b442801bdd1456b05180cbeb640cad7d60fc37) 
<img src="https://live.staticflickr.com/65535/54846895482_2583ef640f_c.jpg">

A remarkable point in the calculation result appear when the pressure intersects the cavitation line. It denotes the gas vaperizes, which is absolutely critical in the water piping system. Adjust the design paramters will contribute to ensure the system safety. 

    
#### Discharge Coefficient
Following the scenario above, the flow rate out of the nozzle may not equal to the velocity by cross section. The geometrical burden can restrict the flow rate in certain degree. The quantification is the discharge coefficient, showing how well the flow can adapt to the geometrical change. Normally, this value is empirical defined, e.g.,
<img src="https://www.numeric-gmbh.ch/fig/discharge_coefficients.png">

The flow rate is then represented as
$$
\displaystyle Q = C_d \  A \ u
$$
where A is the cross section area, and $C_d$ is the discharge coefficient between $0$ (closed) and $1$ (cross section area totally utilized).

__<ins>Scenario: Surge tank fed by oscillation inlet__

Surge tanks are large tanks that serve as a reservoir for excess water near a turbine in a pipeline. Generally, pumps often generate pulsatile flows (especially gear pumps or piston pumps). The surge tank damps out the pulsations before fluid reaches sensitive valves or actuators, see [__video inside of a surge tank__](https://www.youtube.com/watch?v=RenWUexKaac).

   
Imagine a long pipeline from a reservoir to a turbine (or pump). Somewhere along the pipe, a vertical open tank (the surge tank) can be installed. i.e.
- __Normal steady operation__: flow bypasses the surge tank — the tank just has some water level equal to the hydraulic grade line.
- __Disturbance (valve closure, pump surge, load change)__:
    - If the pipe pressure rises suddenly, water flows into the surge tank, raising its level.
    - If the pipe pressure drops suddenly, water flows out of the surge tank back into the pipe.

<img src="https://www.shutterstock.com/image-illustration/water-supply-system-surge-tank-260nw-2434465969.jpg">
    
_Supposed the pipeline pressure head ($H_p$) oscillates around 10 m due to pump fluctuations:_
$$
H_p(t)=10 + 1.5\ \sin⁡(2\pi ft)
$$

_with the surge tank dimension:_
- _tank cross-sectional area: 1 m²_
- _connection area: 0.05 m²_
- _discharge coefficient: 0.6_

_What is the head and tank water level during the oscillation?_

[__<ins>Sample solution__](https://cocalc.com/share/public_paths/b64727c3169d8baac2ac25b1f7a11f4ee1128593)  

    
Implementing the surge tank functionality into a reservoir connected piping system (as the picture shows), a complex system is built. For one, the flow rate from the reservoir into the pipe ($Q_p$) will come both to the dwon stream and the surge tank, depending on the head height (a.k.a. total pressure) on each location.  
<img src="https://live.staticflickr.com/65535/54823644022_827dfbc38b_b.jpg">
Piece by piece analyzing the components of the system, we can start with the pipe section __from reservoir to surge tank junction__ (head in the junction position named as $H_j$). The force equilibrium is formulated as
    $$
    \begin{align}
    &m\frac{du}{dt} = F_\text{pressure} - F_\text{resistance} \\
    \\
    &m = \rho \ A_p L, \qquad u = \frac{Q_p}{A_p}, \qquad m\frac{du}{dt} = \rho L \frac{d Q_p}{dt}\\
    \\
    &F_\text{pressure} = A_p \ (p_R - p_j) = \rho g A_p (H_R -H_j) \\
    \\
    &F_\text{resistance} = R \ Q_p \\
    \\
    &\rho L \frac{d Q_p}{dt} = \rho g A_p (H_R -H_j) - R \ Q_p 
    \end{align}
    $$
and __surge tank__ as well we the flow rate discussed in the scenario.
    $$
    \begin{align}
    &\frac{d\ h_s}{dt}= \frac{Q_s}{A_s}\\
    \\
    &Q_s = C_d A_c \text{sgn}(H_p−h) \sqrt{2g∣H_p−h∣}
    \end{align}
    $$
To describe the downstream flow rate ($Q_d$), the head of the junction ($H_j$)is required, which
    $$
    Q_d = H_j /R_d
    $$
    where $R_d$ is the resistance of the downstream pipe, modeled as linear to the flow rate.
Finally, continuity equatuion has to be applied to keep the two characteristics consistent, which is
    $$
    Q_p = Q_s + Q_d
    $$

Here is a [<ins>__sample implemetation__](https://cocalc.com/share/public_paths/9124600648c819fc85d3c48c21e8f8bff70a8a1a) for the mechanism described above.

Incidents  such as rapidly closing up valve or shutting down downstream pump can be critical for large piping system. These indicate to a sudden discontinuity in the system, which create a local high head. Consequence of local high head is the reverse flow, which can create water hammer to breaks the pipe at it weak points. For such scenarios, surge tanks ease the pressure vialation by surving the addtional head level as buffer zone. A [<ins>__sample implemetation__]() describe this scenario and demonstrate the surge tank functionality by comparing to the no surge tank situation, showing in a sketch.
<img src="https://live.staticflickr.com/65535/54828337499_44789bf7c9_b.jpg" width="75%">

    
---    

### 3.2.3 Lubrication film and Reynolds equation

Flow patterns between two infinite parallel planes can be applied in diverse disciplines. In the mathematical point of view, the Navier-Stokes Equations simplified due to 2-dimensional, steady and incompressible flow patterns yield  
$$
\frac{\partial p}{\partial x}=\mu \frac{\partial^2u}{\partial y^2}
$$
where $p$ is the pressure, $\mu$ is the dynamic viscosity of the fluid, and $u$ is the velocity component in the x-direction.

The general solution of the equation above yield in a parabolic shape to
$$
u(y) = \frac{1}{2\mu}\frac{\partial p}{\partial x}y^2+\tilde{A}y + \tilde{B}\ .
$$
$\tilde{A}$ and $\tilde{B}$ vary from flow type to flow types, and are determined by the boundary conditions.
    
    


#### Couette and Poiseuille flow types

Couette and Poiseuille flows are the classical types. 

<img src="https://www.researchgate.net/publication/251857805/figure/fig11/AS:298025474904069@1448066219094/a-Plane-Couette-flow-upper-plate-moving-at-velocity-V-relative-to-the-lower-plate-b.png" width=65%>


In __Couette flow__, one plate is stationary and the other is in motion. The driving force for fluid motion is the relative velocity between the plates, denoting _no pressure gradient_ plays a role in the movement. As a result, the velocity profile across the gap between the plates is linear, with the velocity increasing linearly from the stationary plate to the moving plate. This linear velocity profile leads to a constant shear stress across the gap, which is proportional to the velocity gradient.  

Taking the $h$ as the height between two planes, the velocity profile yields
$$
u(y) = \frac{U}{h} y\quad,
$$
where $U$ is the velocity of moving plane.

Whereas, the planes are fixed for __Poiseuille flow__ and the pressure gradient is the main driving force for fluid movement. Thus
$$
u(y) = \frac{1}{2\mu}\frac{\partial p}{\partial x}(y^2-hy) \quad.
$$

The velocity profile in the Poiseuille flow across the channel's cross-section is parabolic, with the maximum velocity occurring at the center of the channel and decreasing to zero at the pipe walls. This parabolic velocity profile arises from the balance between the pressure gradient driving the flow and the viscous resistance of the fluid against the channel boundaries.

A general form combining the Couette flow and Poiseuille flow types yields to
$$
u(y) = \frac{U}{h} y + \frac{1}{2\mu}\frac{\partial p}{\partial x}(y^2-hy)
$$
And the shear stress ($\tau_{xz}$) on both walls has the representation
$$
\begin{align}
\displaystyle \tau_{xz}|_{y= 0} = \mu \frac{\partial u}{\partial y} |_{y= 0}&= -\frac{1}{2}\frac{\partial p}{\partial x}h +\mu\frac{U}{h}\qquad \qquad&&\text{fixed plane}\\
\\
\displaystyle -\tau_{xz}|_{y= h} = -\mu \frac{\partial u}{\partial y} |_{y= h}&= -\frac{1}{2}\frac{\partial p}{\partial x}h -\mu\frac{U}{h}\qquad \qquad&&\text{moving plane}\\
\end{align}
$$
As well as the flow rate over cross section $\dot{q}$ as
$$
\dot{q} = \int_0^h u(y) dy= -\frac{h^3}{12\mu}\frac{\partial p}{\partial x} + U\frac{h}{2}
$$
#### Reynolds Equation

Inserting the velocity representation in to continuity equation and integrating over the channel height $h$, it yields to
$$
\frac{\partial}{\partial x}\left(\frac{h^3}{6\mu}\frac{\partial p}{\partial x}\right) = U\frac{\partial h}{\partial x} 
$$
This equation shows the change of pressure gradient when the gap between two planes varies.


Extending to the span-wise direction, density variation and moving boundary for both sides, the [Reynolds Equation](https://en.wikipedia.org/wiki/Reynolds_equation) provides a comprehensive representation showing the pressure distribution within a thin fluid film. In corporated with the characteristics of [non-newtonian fluid](https://en.wikipedia.org/wiki/Non-Newtonian_fluid), the Rynolds equation serves as the fundamentals for the analysis in the lubrication applications.

#### Hydrodynamic Bearing lubrication

Journal bearings are components used to support rotating shafts within machinery. They operate by creating a thin film of lubricating fluid between the rotating shaft and the stationary bearing surface. Classical application are crankshafts of the internal combustion engine, where it transfer the piston generated energy to the kinematic of the vehicle.

<img src="https://www.ctconline.com/media/nz0nv52y/journal-bearing.png?width=634&height=208&mode=max" width=70%> <img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgibwdObG6maNLTpz4ZSt4Is5kl05IhiN7L9-5J1OKBrDX-8r0E5T_hQ8uHJGrgCn4mwo7R1Hn-wP-a0ujEEt_NFGzDtWKPN8kcW-PXmaRk5Zkr5qCB1OrSmPvTlX1zLRJJXLt2ovfBKGoq/s640/4.png" width="28%">

The Reynolds equation is employed to analyze the pressure distribution within this lubricating film and understand the bearing's performance characteristics. 

The clearance between liner and shaft presents in the formulation of angle $(\theta)$ as
$$
h(\theta) = c(1+ \epsilon \cos\theta)\qquad,
$$
where $c$ is the mean clearance ($R_l -R_s$), $\epsilon$ is the ratio of ecentricity $(e)$ to mean clearance $(c)$. 



#### Journal bearing with infinite length 
       
Imaginge an infinite long bearing, where the effects in $z$ direction is negligible. In the representation of $\theta$, the velocity $U$ becomes $\omega R_s$, the $dx$ becomes $R_s d\theta$.

The simplified Reynolds equation yields to 

$$\displaystyle
\frac{\partial}{\partial \ \theta}\left(\frac{h^3}{6 \mu}\ \frac{\partial p}{\partial \theta} \right) = R_s^2\  \omega\frac{\partial h}{\partial \theta}
$$
Integrating over $\theta$, the equation becomes
$$
\frac{\partial p}{\partial \theta} = 6 \mu \ \omega R_s^2\left(\frac{h(\theta)-h_m}{h^3(\theta)}\right)\quad,
$$
where $h_m$ is the clearance thickness by no pressure gradient situation, which can be solved by setting  $p(0) = p(2\pi)$. 

The pressure distribution over the angle yields
$$
p = \int 6 \mu \ \omega R_s^2\left(\frac{h(\theta)-h_m}{h(\theta)^3}\right) d\theta + \bar{C}\qquad
$$
and can be achieved by using [Sommerfeld's number](https://en.wikipedia.org/wiki/Sommerfeld_number).

#### Torque and Frictional coefficient

The frcition of a bearing operation comes from the shear stress. For an infinite long journal bearing, the shear stress and the resulting torque yield to
$$
\tau_{xz} = \mu \frac{\partial u}{\partial y}|_{y=h}\\
$$
and the torque ($T$) is the shear stress integrated over the entire shaft surface
$$
T = \tau_{xz} \ 2\pi R_s L\quad,
$$
where $L$ is the length of the bearing.
Assuming a friction coefficient $(C_f)$ defined as ratio of friction torque to axial load on bearing $(W)$, 
$$
\displaystyle C_f=\frac{T}{WR_s}\quad,
$$
The relation shows a linear behavior of rotational speed frictional coefficient.

__<ins>Scenario: infinite length journal bearing__
    
_A jounral bearing of the following dimensional parameters:_

- dynamic viscosity $\mu$ = 0.03      [$Pa\cdot s$]
- Shaft radius $Rs$: 0.03             [$m$]
- Roataional speed $\omega$: 1.0      [$rad/s$]
- Mean clearance $c_0$: 1    [$mm$]  (Rb -Rs)
- Length of the bearing $L$: 5 [cm]      
- eccentricity ratio $\epsilon$: 0.5  
    
__Numerical solution__

Discretizing the domain $\theta$, the derivative terms of $\displaystyle\frac{\partial p}{\partial \theta}$ and $\displaystyle\frac{\partial^2 p}{\partial \theta^2}$ can be linearized using the neighbouring points By utilizing finite difference methods, [a sample implementation](https://cocalc.com/share/public_paths/3612a76fa5c1c14ece4916d049757ecfe8cfbf5f) in python shows the pressure distribution around the shaft as

<img src="https://live.staticflickr.com/65535/54882018033_1e6300d132_b.jpg">

with

    Net lift force: 0.392 Nt in direction[-0.99996753 -0.00805841] 
    friction torque: -0.39222 Nt.mm

This shows the shaft from this position will be pushed to the left, and results to a new eccentricity.

    
    
#### Hydrostatic Bearing
    
Apart from hydrodynamic bearing, hydraulic bearing is a type of fluid bearing that relies on an external source of high-pressure liquid (typically oil, water, or air) to maintain a continuous, pressurized film of fluid between moving surfacesto prevent solid contact. It works even at no speed condition.

For constant film thickness and incompressible fluid in a hydrostatic bearing at steady state (no sliding-induced hydrodynamics), the Reynolds equation reduces to:
$$
\begin{align}
\quad \frac{\partial}{\partial x_i}(\frac{h^3}{6 \mu}\frac{\partial p}{\partial x_i})  =  0 \qquad\qquad 
⇒⇒ \quad \frac{\partial^2 p}{\partial x_i^2}&= 0\quad (\text{since}\  h=\text{const})
\end{align}
$$

with Dirichlet boundary conditions:

$p= P_\text{out}$ on the outer boundary,
$p=P_p$ (unknown) on the pocket boundary. 
    
In the numerics perspective, the classical Laplace or Poisson type equation can be handled both by FVM and FDM. However, the flux type methodology (FVM) has more flexibility in domain descretization and is rather recommanded. For the laplace equation of a two dimension case ($i: x, y$) for instance, 
<img src="https://live.staticflickr.com/65535/54990793926_2e53c06538_b.jpg" width="80%">
$$
\begin{align}
&\frac{\partial^2 p }{\partial x^2} + \frac{\partial^2 p }{\partial y^2} =      \int_{\partial S} \left(\frac{\partial p}{\partial x}dA_y + \frac{\partial p}{\partial y}dA_x  \right) = 0\\
\\
&= \frac{p_{i+1, j} - p_{i, j}}{x_{i+1, j} - x_{i, j}}\Delta y_E \ + \frac{p_{i-1, j} - p_{i, j}}{x_{i-1, j} - x_{i, j}}\Delta y_W \ +  \frac{p_{i, j+1} - p_{i, j}}{x_{i, j+1} - x_{i, j}}\Delta y_N \ + \frac{p_{i, j-1} - p_{i, j}}{x_{i, j-1} - x_{i, j}}\Delta y_S 
\end{align}
$$
In case of equidistant discretization i.e. $\Delta y_E = \Delta y_W = \Delta y_S = \Delta y_N = \Delta x_E = ...$,  the laplace equation yields:
$$
    p_{i,j} = \frac{1}{4} \left( p_{i+1, j} + p_{i-1, j} + p_{i, j+1}+ p_{i, j-1}\right)
$$

Incorporated with the boundary condition, this equation can be taken iteratively to approach the solution, e.g., by using [Jacobi method](https://en.wikipedia.org/wiki/Jacobi_method).

__Hydrostatic Thrust Pad Bearing__ 
Hydrostatic thrust pads are commonly applied in large turbines, compressors, and heavy-duty pumps to support significant axial loads with minimal friction. They are also used in precision machine tools and high-speed spindles, where high stiffness and vibration damping are critical for accuracy and reliability. A sketch is shown:
    
<img src="https://live.staticflickr.com/65535/54988970003_c8a3c15316_b.jpg">

__<ins> Scenario: Tiled hydrostatic thrust bearing__    
Say a thrust pad has the technical dimensionas followed:
- Pad dimensions: $L_x$=$L_y$=100 mm
- Pocket (recess): $a_x$=$a_y$=50 mm (centered)
- Supply pressure: $P_s$=3.0 MPa (gauge)
- Ambient pressure: $P_{\text out}$=0 (gauge)
- Lubricant viscosity $\mu\approx$ 0.040 Pa$\cdot$ps, $\rho\approx$ 860 kg/m3
- Restrictor: sharp-edged orifice, $d_r$ = 1.0 mm, discharge coefficient $C_d$ = 0.62

Due to the structural, material and manufactural constrain, thrust pad bearing may deform under specific circumstances. For such cases, pad is tiled to the loading surface, for which the performance can be very different. 
    
Under this aspect, how is the distribution if
a) __film thickness is uniform ($h$=300 μm) without pad movement (design state)__
    <img src="https://live.staticflickr.com/65535/54988800516_012f3eca09_b.jpg" width="50%">
([Sample implementation](https://cocalc.com/share/public_paths/a036b63fbf49809a954ab978af9c96b22989b008))
b) __non-uniform film thickness due to inclined of loading plate__ ($\alpha = 5e^{-3}$ m/m), 
<img src="https://live.staticflickr.com/65535/55033228607_5ab3ef86de_b.jpg" width="100%">
and
c) __non-uniform film thickness while the pad is moving__ ($U$ = -2 m/s)?
<img src="https://live.staticflickr.com/65535/55034308313_5d15f49a76_b.jpg" width="100%">
    


---
[__BACK TO CONTENT__](https://hackmd.io/@SamuelChang/H1LvI_eYn)