# Sound generation mechanisms in a collapsible tube

Marco<EMAIL_ADDRESS>The following article has been submitted to
The Journal of Acoustical Society of America. After it is published, it will
be found at https://pubs.aip.org/asa/jasa. FLOW Research Center, Department of
Engineering Mechanics, KTH Royal Institute of Technology, Stockholm, SE-10044,
Sweden.  The Marcus Wallenberg Laboratory for Sound and Vibration Research,
Department of Engineering Mechanics, KTH Royal Institute of Technology,
Stockholm, SE-10044, Sweden  Elias Zea  The Marcus Wallenberg Laboratory for
Sound and Vibration Research, Department of Engineering Mechanics, KTH Royal
Institute of Technology, Stockholm, SE-10044, Sweden  Elias Sundström  FLOW
Research Center, Department of Engineering Mechanics, KTH Royal Institute of
Technology, Stockholm, SE-10044, Sweden.  Susann Boij  The Marcus Wallenberg
Laboratory for Sound and Vibration Research, Department of Engineering
Mechanics, KTH Royal Institute of Technology, Stockholm, SE-10044, Sweden
Mihai Mihaescu  FLOW Research Center, Department of Engineering Mechanics, KTH
Royal Institute of Technology, Stockholm, SE-10044, Sweden.

###### Abstract

Collapsible tubes can be employed to study the sound generation mechanism in
the human respiratory system. The goals of this work are (a) to determine the
airflow characteristics connected to three different collapse states of a
physiological tube and (b) to find a relation between the sound power radiated
by the tube and its collapse state. The methodology is based on the
implementation of computational fluid dynamics simulation on experimentally
validated geometries. The flow is characterized by a radical change of
behavior before and after the contact of the lumen. The maximum of the sound
power radiated corresponds to the post-buckling configuration. The idea of an
acoustic tube law is proposed. The presented results are relevant to the study
of self-excited oscillations and wheezing sounds in the lungs.

##  1 Introduction

In daily medical practice, auscultation is among the most cost-effective and
timely diagnostic tools for diseases of the respiratory and cardiovascular
system [29]. Different pathologies are characterized by specific sounds
(sometimes called adventitious sounds) that can guide the physicians toward
identifying and assessing the disease. From an acoustic perspective, the
problem of sound generation in the human body and its propagation represents a
challenging question that is far from being completely understood [37, 35].

From a biomechanical perspective, for respiratory and cardiovascular
pathologies, the generation of adventitious sounds is often correlated to the
interaction between the biological fluid (air or blood) and the human vessel
(airways, or veins and arteries) [30]. The interest in quantitatively modeling
the Fluid-Structure Interaction (FSI) phenomenology and the corresponding
sound generation mechanisms in the context of different diseases is justified
by the ambition to develop and improve sound-based diagnostic tools. In this
perspective, thanks to more accessible large-scale high performance
computational resources, Computational Fluid Dynamics (CFD) represents a
powerful tool able to furnish a 4D description of mass transportation in the
human body and the corresponding FSI dynamics [36].

The current computational power allows to study cardiovascular [46] and
respiratory [39, 54] diseases in both patient-specific [27, 53] and idealized
geometries [59, 57]. One remarkable example of an idealized geometry model of
a human vessel is a collapsible tube. Collapsible tubes can reproduce all the
relevant physical behavior of a human vessel and are widely employed in both
numerical models [51] and clinical practice [32]. The advantage of idealized
geometries for studying pathophysiological flows is not only found in the
lesser computational cost. Numerical models based on CFD simulation with the
ambition to provide a reliable quantitative analysis need to be experimentally
validated. The possibility to access the large amount of independent high-
quality experimental data [42, 14, 45] targeting collapsible tubes under
different conditions makes the validation step possible and reliable.
Moreover, by reducing the complexity of a biological system to its very
fundamental parts, it is possible to identify the main physical mechanisms
driving its dynamics. For these reasons, collapsible tubes play a major role
in analyzing many different pathologies of the cardiovascular and respiratory
systems.

One example of application of collapsible tubes in biomedical acoustics is the
study of respiratory wheezes. Wheezing is one of the most common adventitious
sounds generated in the lungs. It is connected to several different
pathologies of human lungs like asthma [21] and chronic obstructive pulmonary
disease [19]. Tonal features characterize it and can happen during both
inhalation and expiration. The sound generation mechanism of wheezing sounds
is associated with the onset of self-excited oscillations in the airways due
to the interaction between air and the conduct [22]. It is possible to study
such FSI-induced oscillations by implementing numerical [24] and experimental
[18] studies involving collapsible tubes. However, the connection between the
FSI and the corresponding generated acoustic signal is still to be completely
understood [60].

Another medical research field where collapsible tubes are widely employed is
the study of the collapse of the upper airways in pathologies like Obstructive
Sleep Apnoea (OSA). OSA patients suffer from recurrent pharynx collapses
during sleep, resulting in apnoea and consequent neurological disturbances
[52]. The diagnosis and assessment of the OSA severity is performed through an
expensive procedure called polysomnography. The main goal is to estimate the
level of collapsibility of the pharynx in terms of the so-called critical
pressure, i.e., the value of the intramural pressure corresponding to the
collapse of a compliant tube [47]. Collapsible tubes are also widely employed
in the study of the circulatory system [28]. One interesting research line is
connected to pulse wave imaging ultrasound in stenotic arteries [58], in which
the carotid is modeled using a collapsible tube undergoing a pulse
displacement. Other research directions include the analysis of stenotic
arteries [48], venous collapse prevention [49], the coupling between cardiac
cycle and coronary perfusion [50], and the study of choroidal blood flow [40].

Several studies in the acoustics literature have used collapsible tube models
for applications to voice generation. Collapsible tubes have been used to
model the vocal cords using an electrical analogy [5]. Since the tube behaves
as a flow-controlled nonlinear resistance for steady flows, similar to a
transistor, Conrad tested whether vocal cords have the same property. Titze
proposed a physical model for small-amplitude oscillations in the vocal folds
based on the body-cover hypothesis [7], which is supported by experimental
observations in excised larynges—a physical model of the vocal-fold mucosa,
and in human subjects [13]. In the context of laryngeal cancer, Tourinho et
al. have used a collapsible tube model for tracheoesophageal (TE) phonation
[44], a common speech rehabilitation method. The results showed that the tube
model reproduced the fundamental frequency and spectral characteristics of TE
speech. Recently, Li et al. trained a neural network to predict the glottal
flow waveform [43], using the collapsible tube model from Cancelli and Pedley
[6]. Li’s model was validated using experimental data from rabbits with
excised larynx.

In most of the studies in the literature, collapsible tubes have been employed
to study biological fluid flows, self-excited oscillations, collapse critical
pressure, or to produce acoustic reduced-order models. However, the nature of
the fluid-borne sound generation mechanisms in collapsible tubes relevant for
biomedical applications and their relation with fluid mechanics structures is
still an open question. This work addresses these problems by studying the
fluid behavior in three distinguished collapse phases: pre-buckling, post-
buckling, and post-contact. The methodology is based on CFD and Computational
Aeroacoustics (CAA) modeling. The two goals of this work are: 1) to determine
the main fluid characteristics related to sound generation under different
collapse conditions, and 2) to relate the overall sound power level generated
by the flow with the collapsible state of the tube. This last point is also
motivated by a similar experimental investigation on urethral obstruction
modeled as a collapsible tube [10].

The structure of the paper is as follows. Sec. 2 introduces the phenomenology
of a collapsible tube under physiological conditions. In Sec. 3, the details
of the numerical model implemented in this work are discussed. Sec. 4 concerns
the sound generation mechanisms due to the interaction of the fluid flow with
the geometries corresponding to three different stages of the collapse. In
Sec. 5, the analysis of the acoustic power radiated by the collapsible tube
under the three different collapse phases is discussed. The onset of acoustic
waves in the tube domain is investigated in Sec. 6. Finally, discussion and
conclusions can be found in Sec. 7.

##  2 Collapsible tube

The geometry of a collapsible tube is described by three non-dimensional
parameters (see Fig. 1). The length-to-diameter ratio $d=l_{0}/D$, the
thickness-to-diameter ratio $\gamma=h/D$, and the axial pre-stretch ratio
$l=L/l_{0}$.

Figure 1: (Color online) Top panel: radial cross-sectional view of the domain
pictured in its post-buckling configuration. The (half) diameter $D/2$ and the
thickness $h$ are sketched. Bottom panel: lateral view of the actual
configuration of the domain under the action of the boundary conditions
$p_{ext}$ (red dashed arrows) and $l$ (light blues arrow). The two black solid
lines indicate the clamped boundary. The blue box corresponds to the initial
configuration of the system.

The values of these parameters have been chosen to be relevant for biomedical
applications [1, 3] and are listed in Tab. 1.

In this study, one short side of the tube is clamped while the other is
stretched by an amount indicated by the pre-stretch parameter $l$ and then
clamped. Under these assumptions, the collapse is only driven by the so-called
intramural pressure i.e., the pressure difference between the interior and the
exterior of the tube:

$p_{intr}=p_{int}-p_{ext}\,.$ (1)

When this difference is negative, the tube starts to collapse. The collapse is
described by the so-called tube law [23], which relates the area of the
central cross-section of the tube and the intramural pressure (see Fig. 2).
The tube collapses for small negative values of the intramural pressure by
maintaining its axisymmetric shape (pre-buckling phase). When the intramural
pressure assumes larger negative values, the tube undergoes a buckling
phenomenon, resulting in a non-axisymmetric shape (post-buckling phase). Under
most physiological conditions, the tube’s cross-section is characterized by a
2-lobe shape [26]. Finally, for even larger negative values of the intramural
pressure, the internal walls of the tube touch each other (post-contact phase)
[2]. At first, the contact region is one point, and it becomes a line for more
negative intramural pressures (see Fig. 2).

Figure 2: (Color online) The tube law. The right axis indicates the intramural
pressure, while the top axis indicates the area of the central cross-section
normalised on the initial area. The green areas represents the regions of the
transition from the pre-buckling to the post-buckling phase and to the post-
contact configuration.

The transition between the three phases of the collapse (pre-buckling, post-
buckling, and post-contact) is marked by two particular values of the
intramural pressure called buckling critical pressures and contact critical
pressure. The value of such critical pressures depends on both geometric and
elastic parameters of the system [31, 45]. In [56], the authors have shown
that the buckling critical pressure can be estimated by treating the buckling
as a second-order phase transition (see also [41] for a rigorous proof for a
1D ring). In [55], a treatment for estimating the contact critical pressure is
derived.

The fluid behavior in a collapsible tube under different collapse scenarios
has been extensively studied from theoretical [4], numerical [12], and
experimental [42] points of view. However, understanding the relation between
the tube’s characteristic flow structures and the corresponding acoustic
features remains an open query. Some relevant studies have been conducted in
[33, 16], where the authors have studied the aeroacoustics of an orifice plate
in a duct using the two-port scattering method. This method, however, targets
the characterization of the scattering matrix of the system for an already
existing acoustic signal [11]. Moreover, the system’s geometry under analysis
does not reflect the 3D behavior of a human vessel.

The present work lies in the larger project of establishing a relation between
the acoustics of a human vessel and its collapse level, with the potential
development of new diagnostic tools for respiratory and cardiovascular
diseases. The following section presents the numerical model implemented in
this study. The goal is to perform an in-depth analysis of the fluid behavior
under the three collapse phases (pre-buckling, post-buckling, and post-
contact) of a collapsible tube. The flow features responsible for the onset of
acoustic waves propagating in the tube are studied. Finally, the relation
between the acoustic power level generated in the system and the tube’s
collapse level is presented, supporting the existence of an acoustic tube law.

$d$ [-] |  $\gamma$ [-] |  $l$ [-]
---|---|---
$3$ $(2.5-6)$ |  $0.06$ $(0.03-0.1)$ |  $1.25$ $(1-1.6)$
Table 1: Values of the geometric parameters of the system. The values in
parenthesis correspond to the physiological range in the human airways.

##  3 Numerical model

It is possible to study the airflow under the three collapse states of a
physiological tube and the corresponding acoustic field by employing numerical
methods. The numerical model is implemented in the commercial software Siemens
Star-CCM+ (v. 2210). The simulation strategy consists of two steps:

  1. 1.

Solid-only simulation of a collapsible tube under the effect of an isotropic
external pressure. The goal is to compute the resulting deformed
configurations of the tube corresponding to the pre-buckling, post-buckling,
and post-contact phases of the collapse. These will be compared for validation
with experimental results and used as fluid domain geometries in the second
step.

  2. 2.

Simulations of the velocity and pressure field of air flowing in the three
fixed collapsed configurations obtained in the previous step. Since air is
treated as a compressible ideal gas, the resulting acoustic pressure
fluctuations are predicted at the same time.

###  3.1 Solid elastic model

The reference configuration geometry is a 3D cylindrical flexible tube with
finite thickness and it is implemented in the built-in CAD software in Star-
CCM+. The values of the geometric parameters $(d,\gamma,l)$ are listed in Tab.
1 and are typical for human airways. The system is subject to the following
boundary conditions (see Fig. 1). One short side of the tube is clamped
whereas the other one is first stretched by a quantity related to the
parameter $l$ and then clamped. The external walls of the domain are under the
effect of an isotropic inward pressure $p_{ext}$ which linearly increases in
time. Due to the absence of flow, the internal pressure vanishes, $p_{int}=0$,
and the intramural pressure depends only on the external pressure,
$p_{intr}=-p_{ext}(t)$. The following relation describes the time dependence
of the external pressure:

$p_{ext}(t)=\frac{P}{\tau}t\,,$ (2)

where $P$ is the maximum external pressure and $t\in[0,\tau]$ are the time
steps. By changing the value of $P$, it is possible to obtain, at time
$t=\tau$, the desired collapse state of the tube. If $P<P_{cr}^{b}$, i.e.,
smaller than the buckling critical pressure, the tube will be in its pre-
buckling phase. Analogously, if $P_{cr}^{b}<P<P_{cr}^{c}$ or $P_{cr}^{c}<P$,
the resulting collapse state will correspond to the post-buckling and post-
contact phase, respectively, where $P_{cr}^{c}$ is the contact critical
pressure. The values corresponding to the three different collapse states are
$1000$ Pa, $2000$ Pa, and $3000$ Pa for the pre-buckling, post-buckling, and
post-contact configuration, respectively (see also Fig. 3) and are compatible
with typical physiological conditions of the airways. Indeed, under normal
breathing conditions, the intramural pressure is approximately $2000$ Pa [9].
However, for asthma patients, for which the collapse of the airways is often
manifested, these values can reach up to $10000$ Pa [20], depending on the
size of the patient [17]. For values of $p_{ext}>P_{cr}^{c}$, the internal
walls of the tube (the lumen) touch each other. The contact is handled
employing a repulsive virtual plane [55], which prevents the intersection of
the solid domain. A sensitivity study of the ratio $P/\tau$ has been performed
in [56], and it is not reported here for brevity. For all the following
analyses, $\tau=2.5$ s is used.

Under such boundary conditions, the system undergoes large deformation that
cannot be adequately treated with linear elasticity [56]. Therefore, the tube
is modeled as a Neo-Hookean material whose strain energy potential is defined
as

$\mathfrak{U}=\Psi(I-3)+\frac{\Phi}{2}(\tilde{I}-1)^{2}\,,$ (3)

where $I$ and $\tilde{I}$ are the first and second invariants of the right
Cauchy-Green tensor, respectively. The coefficients $\Psi$ and $\Phi$ are
defined as

$\Psi=\frac{E}{4(1+\nu)},\qquad\Phi=\frac{E\nu}{(1-2\nu)(1+\nu)}\,,$ (4)

where $E$ is the Young’s modulus and $\nu$ is the Poisson ratio. The values of
these elastic parameters are the same as the samples used in [42], whose
experimental data will be used to validate the numerical results.

At each time step ($\Delta t=0.1$ s), the value of the intramural pressure and
the area of the central cross-section are registered. It is then possible to
determine the corresponding tube law. The validation of the numerical model
has been performed by comparison with a public experimental data set [34, 42]
obtained via 3D camera measurements of the area of the central cross-section
of the tube. A digital replica of the experimental specimen has been
implemented in Star-CCM+, and the value of the maximum pressure $P$ has been
set to match the experimental conditions. The comparison between the numerical
and experimental tube laws shows a fair match (see Fig. 3).

Figure 3: (Color online) Comparison between the tube law obtained via the
presented numerical model (red dashed line) versus the corresponding
experimental data from [42].

The output of these solid-only simulations consists of three experimentally
validated deformed configurations of a collapsible tube, corresponding to the
pre-buckling, post-buckling, and post-contact phase (see Fig. 4). Such
deformed configurations will be employed as numerical domains for the
subsequent fluid flow and acoustic field simulations.

Figure 4: From left to right, the cross-sections of the fluid domains obtained
from the solid model simulations for the pre-buckling, post-buckling, and
post-contact configurations.

###  3.2 Fluid model

The three geometrical computational domains depicted in Fig. 4 are used to
simulate the air velocity and pressure fields. The domain is extruded from the
collapsed tube in the upstream and downstream directions to prevent any
influence of the imposed boundary conditions on the region of interest. This
results in a division of the computational domain into three regions (see Fig.
5): inlet, sound generation and propagation region, and acoustic suppression
zone (ASZ) outlet. The sound generation and propagation region defines the
region of interest. The inlet and outlet are treated as free stream
boundaries, which allow to impose the flow velocity ($10$ m/s) and static
temperature ($308$ K). One advantage of free stream boundaries is that they
are perfect non-reflective boundaries for plane waves. As discussed in the
following, the geometry and boundary conditions of the problem allow only
plane waves to propagate in the acoustic frequency range. The flow pressure
fluctuations, characterized by a much smaller length scale, are dissipated by
the mesh stretching in the ASZ and, therefore, will not be reflected. All the
other walls are treated with no-slip boundary conditions.

Figure 5: (Color online) Lateral view of the computational domain of the fluid
numerical model. In green is highlighted the inlet region. The blue rectangle
shows the region of interest and represents the sound generation and
propagation region. Finally, the acoustic suppression zone is highlighted in
red.

The simulation solver and mesh grid are tuned to perform an unsteady
compressible flow Direct Noise Calculation (DNC), which allows the whole
aeroacoustics characterization of the system to be performed within one
simulation. Both the unsteady flow, which causes the noise generation, and the
resulting radiated sound waves are solved simultaneously. This approach
imposes precise requirements on the mesh grid, which is built in three steps:

  1. 1.

A coarse unstructured mesh of polyhedral elements is defined to discretize the
system. It allows to obtain a well-converged preliminary solution based on the
Reynolds Averaged Navier-Stokes (RANS) approach that will be used as the
initial condition for the following unsteady simulation.

  2. 2.

Since the time evolution of the fluid features is resolved via a segregated
Large Eddy Simulation (LES) model with a WALE sub-grid scale model, a first
mesh refinement is implemented in the region of interest to resolve the Taylor
micro-scale. Moreover, 10 prism layers are implemented on the walls to ensure
well-resolved boundary layers.

  3. 3.

To ensure that all the waves in the acoustic frequency range can propagate in
the domain, the mesh is further refined to ensure that at least $20$ grid
cells resolve the shortest acoustic wavelength.

Outside the region of interest, the mesh is coarser and is stretched in the
ASZ. In this way, the smaller-scale aerodynamic fluctuations are numerically
dissipated away from the region of interest and cannot reflect into the
domain.

The choice of the time step can also influence the numerical dissipation of
acoustic waves. An implicit second-order time marching scheme is implemented
in the model. A minimum of $15$ time steps per period of the highest frequency
sound wave are recommended for a DNC simulation. For this simulation, the time
step has been set to $\Delta t=3.1\times 10^{-6}$ s, corresponding to about
$16$ time steps for a $20000$ Hz acoustic wave.

The presented fluid model has been implemented in the commercial software
Star-CCM+ by Siemens. The fluid flow LES based solver has been previously
validated in similar geometries by comparing the pressure predictions to
corresponding experimental mid-line pressure data [38].

Velocity and pressure fields are extracted to calculate time-averaged and
fluctuating components on relevant cross-sections for the following post-
processing. The full-time history of pressure and velocity fields on the
domain surfaces are stored at each time step to compute the surface Fourier
transform. Several monitoring points for both velocity and pressure fields are
placed in the domain.

##  4 Sound generation mechanisms

The maximum velocity magnitude of the airflow under the three different
collapse states of the tube is of the order of $40$ m/s, corresponding to a
low Mach number regime ($\text{Ma}\simeq 0.12$). Under such conditions, sound
is generated mainly by the unsteady pressure loads caused by the flow vortical
structures washing the walls of the domain. Therefore, analyzing such
interaction in the three geometries under investigation is interesting. To
this end, the distribution of the mean axial component of the velocity is
analyzed on two orthogonal longitudinal cross sections (see Figs. 6-7),
corresponding to the lateral view (YZ-plane) and the top view (XZ-plane).

The Poiseuille-like flow in the pre-buckling configuration of the tube does
not show any interesting flow feature. It is irrelevant for sound generation
as it shows negligible time-dependent interactions with the walls.
Consequently, it will not be further discussed in this section. The flow in
the post-buckling configuration (Fig. 6) is characterized by a clear Coandă
effect. This effect describes the tendency of a jet flowing from an orifice to
attach to one side of the domain, developing a recirculation region on the
other side. As shown in the YZ-plane cross-section (lateral view, upper panel
in Fig. 6), the buckled configuration acts on the fluid like a sudden
contraction-sudden expansion zone. The resulting jet stays attached to the
lower part of the domain (due to the Coandă effect), producing a large
recirculation area in the YZ-plane. As discussed in Sec. 6, the Coandă effect
is responsible for the onset of acoustic waves in the domain. The cross-
section in the XZ-plane (top view, lower panel in Fig. 6) shows the symmetric
behavior of the airflow through the two lobes of the buckled configuration.

Figure 6: (Color online) Mean axial velocity in the two main longitudinal
cross-sections of the domain in post-buckling configuration. The flow is
characterized by the presence of Coandă effect. The arrows are obtained via
line integral convolution of the velocity field and represent the mean
velocity.

In the post-contact configuration of the tube, the lumen is partially closed.
Such a critical change in the geometry corresponds to a clear change in the
fluid behavior. Due to the contact of the internal walls of the tube, the
Coandă changes its direction (see the upper panel Fig. 7, lateral view). As
shown in the cross-section in the XZ-plane (lower panel in Fig. 7, top view),
the contact region acts like an adverse pressure gradient region, causing the
separation of the boundary layers from the wall. Consequently, a large
recirculation area is produced, constraining the flow towards the lateral
walls of the domain and, in turn, generating sound.

Figure 7: (Color online) Mean axial velocity in the two main longitudinal
cross-sections of the domain in post-contact configuration. The empty regions
correspond to the location of the contact of the lumen. The arrows are
obtained via line integral convolution of the velocity field and represent the
mean velocity.

Interestingly, the topology change from the post-buckling to the post-contact
configuration of the collapsible tube induces a rotation of $90^{\circ}$ of
the main recirculation area. The periodic flow separations in the post-
buckling and post-contact configurations cause an oscillating force on the
fluid in the perpendicular direction to the flow. Such unsteady force is the
primary sound generation mechanism in a collapsible tube under physiological
conditions. It is possible to estimate the unsteady pressure loads developed
on the walls in the three collapse configurations by analysing the pressure
root means square (RMS) on a probe line located on one boundary of the domain
(see the bottom panel in Fig. 8). Remarkably, the post buckling configuration
shows much higher values for the pressure fluctuations on the wall, suggesting
that this configuration generates higher sound power level with respect to the
post-contact configuration. To prove this point, an in-depth sound power
analysis is presented in the next section.

Figure 8: Top panel: Pressure RMS comparison between the three different
geometries. The post-buckling configuration shows the largest pressure
fluctuations on the wall. Bottom panel: position of the probe line where the
pressure RMS has been calculated.

##  5 Sound power analysis

The analysis of the velocity field discussed in the previous section shows
that the flow characteristics and the corresponding sound generation
mechanisms depend dramatically on the collapse state of the tube. The sound
power level can be computed by processing the results of the unsteady
compressible CFD simulations, analysis carried out for all three different
collapsed geometries under analysis. The analysis is based on the employment
of the source field Lighthill’s acoustic analogy. This entails that any
conversion of acoustic energy into flow kinetic energy is neglected, as well
as the sound field interference on the fluid flow. It is also assumed that the
effects of the fluid motion on the sound propagation in the source region are
not considered.

The present section contains an analysis of the sound power level generated by
the fluid flow in the near-fields under the different collapse phases. The
investigated domain can be interpreted as the acoustic near-field. Indeed, as
will be shown in the next section, the wavelength corresponding to the
acoustic pressure fluctuations is of the same order as the longitudinal
dimension of the domain. The following section will present the study of the
subsequent propagation of acoustic waves in the domain. It will exploit the
full power of the DNC simulations, which describe the coupling between the
sound and the flow fields.

Lighthill’s equation is a reformulation of the compressible Navier-Stokes
equations in the form of an inhomogeneous wave equation expressed in terms of
acoustic pressure fluctuation:

$\begin{split}&\frac{1}{c_{0}^{2}}\frac{\partial^{2}p^{\prime}}{\partial
t^{2}}-\frac{\partial^{2}p^{\prime}}{\partial x_{i}\partial x_{i}}=\\\
&\frac{\partial}{\partial t}\left(\frac{1}{c_{0}^{2}}\frac{\partial}{\partial
t}(p^{\prime}-c_{0}^{2}\rho^{\prime})\right)-\frac{\partial g_{i}}{\partial
x_{i}}+\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}(\rho
u_{i}u_{j}-\tau_{ij})\,.\end{split}$ (5)

Here, the prime variables $\rho^{\prime}$ and $p^{\prime}$ are acoustic
fluctuations of the density and pressure mean-field, respectively, $c_{0}$ is
the speed of sound, $g_{i}$ represents the forces acting on the fluid, and
$\tau_{ij}$ is the deviatoric stress tensor. The source term in Eq. 5 can be
decomposed into the following terms:

$\begin{split}&s_{m}=\frac{\partial}{\partial
t}\left(\frac{1}{c_{0}^{2}}\frac{\partial}{\partial
t}(p^{\prime}-c_{0}^{2}\rho^{\prime})\right)\\\ &s_{d}=-\frac{\partial
g_{i}}{\partial x_{i}}\\\ &s_{q}=\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}(\rho u_{i}u_{j}-\tau_{ij}).\end{split}$ (6)

These terms correspond to the $n=0,1,2$ terms of a multi-pole expansion of
order $2^{n}$, i.e., monopole, dipole, and quadrupole source terms. The
monopole term $s_{m}$ describes sound source mechanisms due to deviation from
adiabatic assumption, and it vanishes in the present analysis. The low Mach
number regime of the system implies that the contribution of the quadrupole
forces (which describe the sound produced by turbulence) is negligible
compared to the dipole term. Indeed, the ratio of sound power generated by a
dipole source to the one generated by a quadrupole source scales as Ma-2,
which implies that the contribution due to quadrupole sources to the total
sound power level is in this case of about $2\%$. Consequently, the main
contribution to the sound power generated by the flow in a collapsible tube
under physiological conditions comes from the dipole source $s_{d}$.

The dipole source term describes the sound produced by the fluctuating forces
due to flow separation from a rigid wall. As anticipated in the previous
section, the periodic character of the developed flow structures in the
generated shear-layer, as well as the flapping nature of the developed shear-
layer washing the wall, causes an oscillating force on the walls in the
perpendicular direction with respect to the flow. Another consequence of the
low Mach number regime is that the collapsible tube acts as a compact acoustic
source. Indeed, if $T$ is the typical temporal scale for sound production and
$U$ is the typical mean flow velocity, the typical acoustic frequency is
$f=1/T\sim U/D$, where $D$ is the diameter of the tube. The corresponding
wavelength is $\lambda=c_{0}/f\sim D/\text{Ma}$, which implies that the source
region is much smaller than the typical wavelength for low Mach numbers. Under
these assumptions, the sound power radiated by a dipole source can be written
as [15]

$W_{d}=\frac{\overline{\dot{F}^{2}}}{12\pi\rho_{0}c^{3}_{0}}$ (7)

where the dot represents the time derivative, and $F^{2}=|\textbf{F}|^{2}$,
i.e., the (squared) magnitude of the total force acting on the fluid,
$F_{i}=\int_{V}g_{i}dV$, where $V$ is the volume of the fluid domain. It is
possible to compute this term by extracting the total force acting on the
fluid from the CFD simulations for the pre-buckling, post-buckling, and post-
contact configurations. The resulting sound power levels are plotted in Fig.
9.

Figure 9: (Color online) Comparison of sound power level (SWL) defined as
$10\log(W_{d})$, where $W_{d}$ is defined in Eq. 7 in the three collapse
states under analysis. The value for the pre-buckling is smaller than the
reference and has been set to zero.

Interestingly, the maximum of the sound power level corresponds to the post-
buckling configuration. The explanation is related to the stronger interaction
between the flow and the walls in the post-buckling configuration. Indeed, as
discussed in the previous section, the numerical model predicts larger root
mean square pressure, with respect to the post-contact configuration, in the
region located downstream of the minimum cross-sectional area where the
developed shear-layers are interacting with the walls (see Fig. 8). Larger
pressure fluctuations result in higher values of $\overline{\dot{F}^{2}}$ in
Eq. 7 and, in turn, in larger sound power levels. Such observation implies an
optimal collapsed tube configuration, which yields maximum sound power between
the post-buckling and the post-contact phase. This suggests the existence of
an acoustic tube law, relating the intramural pressure (or the cross-section
area) to the sound power radiated by the collapsible tube. The derivation of
the acoustic tube law will be the object of future work, in which the system
will be studied in an FSI setting with a continuous collapse. It is important
to remark that the value computed using Eq. 7 represents the total acoustic
energy generated by the flow. However, it does not specify how such an energy
will be radiated. Partially, it will generate acoustic waves, partially it
will excite the domain walls, resulting in aeroelastic perturbations, and
partially will be dissipated. The far-field propagation of the sound, both
within the lumen and outside, will be investigated in future work.

##  6 Onset of acoustic waves

The sound generation mechanisms described in Sec. 4 induce pressure
fluctuations in the collapsible tube. The DNC setup implemented for the
numerical simulations allows to analyse at the same time the aerodynamic and
acoustic components of such pressure fluctuations. This section analyzes the
spectral properties of the pressure field in the tube under different collapse
states. Moreover, the onset of acoustic waves is observed.

Probing points for the pressure field are positioned in sensible positions in
the domain to capture the fluctuations associated with the vortical structures
developed in the shear-layer (see Figs. 10-12). From the time history of the
pressure signals, it is possible to obtain the corresponding power spectral
density (PSD) for the three different collapse states of the tube.

Figure 10: (Color online) Upper panel: Pressure PSD magnitude for the pre-
buckling configuration. Lower panel: instantaneous axial velocity and probe
positions in the domain. The probes’ colors correspond to the colors in the
PSD plots. Figure 11: (Color online) Upper panel: Pressure PSD magnitude for
the post-buckling configuration. Lower panel: instantaneous axial velocity and
probe positions in the domain. The probes’ colors correspond to the colors in
the PSD plots. Figure 12: (Color online) Upper panel: Pressure PSD magnitude
for the post-contact configuration. Lower panel: instantaneous axial velocity
and probe positions in the domain. The probes’ colors correspond to the colors
in the PSD plots.

A first remark is that the PSD for the pre-buckling configuration has
negligible amplitudes compared to the more collapsed states of the tube, as
expected from the negligible unsteadiness of the flow behavior in this
configuration. The second observation is that all the PSD plots show narrow-
band peaks. The corresponding frequencies are related to the small-scale
fluctuations in the shear layers associated with the separation mechanisms
described in Sec. 4. The shear layer fluctuations are visible in
correspondence of the probe points in the axial velocity snapshot plots in
Fig. 11 and 12.

The PSD for the post-buckling configuration shows clear asymmetric amplitudes
at the peak frequency $f=5630$ Hz for the two near-field probes, which is
explained by the asymmetric flow induced by the Coandă effect. Such asymmetry
is not detected in the PSD levels for the post-contact configuration due to
the disappearance of the Coandă effect, which induces a more symmetric flow
behavior. It is interesting to notice that although the power density at the
peak frequency in the post-contact configuration ($f=6250$ Hz) is higher than
the corresponding one in the post-buckling configuration, the analysis in Sec.
5 shows that the overall energy of the pressure fluctuation in the post-
buckling configuration is larger.

The value of the peak frequencies in the post-buckling and post-contact
configuration can be compared by means of the corresponding dimensionless
frequency (Strouhal number) defined as

$\tilde{f}=\frac{fD^{\prime}}{U}\,,$ (8)

where $D^{\prime}=\sqrt{2A_{min}/\pi}$ is the effective diameter, $A_{min}$ is
the area of the central cross-section of the tube, and $U$ is the bulk flow
velocity ($U=22$ m/s). The values of Strouhal number for the post-buckling and
post-contact cases are $\tilde{f}_{PB}=0.68$ and $\tilde{f}_{PC}=0.52$,
respectively. This difference is likely attributed to the different flow
mechanisms related to the post-buckling and post-contact collapse state.
Indeed, the flow behavior changes completely as the tube transitions from the
post-buckling to the post-contact configuration.

To study the onset of acoustic waves within the domain induced by the flow
features presented in Sec. 4, it is possible to compute the Fourier transform
of the pressure on the two main longitudinal cross-sections. The computation
is performed using a built-in implementation in Star-CCM+. The pressure data
is windowed with a Hann function, and the number of analysis blocks is $6$
with an overlap factor of $0.25$. Once the Fourier transform is obtained, it
is possible to compute the corresponding real and imaginary parts evaluated at
the two peak frequencies for the post-buckling and post-contact
configurations, $f_{PB}=5630$ Hz and $f_{PC}=6250$ Hz, respectively. These can
be used as the coefficients of the components of the Fourier series of the
pressure time evolution $p(t)$ at the frequency $f_{PB}$ and $f_{PC}$. In
formulae:

$p(t)=\text{Re}\left[\hat{p}(f_{i})\right]\cos\left(\frac{2\pi
t}{T}\right)+\text{Im}\left[\hat{p}(f_{i})\right]\sin\left(\frac{2\pi
t}{T}\right)\,,$ (9)

where $\hat{p}(f_{i})$ is the Fourier transform of the pressure signal
evaluated at the frequency $f_{i}$, where $i=\\{PB,PC\\}$. The corresponding
pressure fluctuations at specific peak frequencies on the two main
longitudinal cross-sections is displayed in Fig. 13 and in Fig. 14, for the
post-buckling ($f_{i}=f_{PB}$) and the post-contact ($f_{i}=f_{PC}$)
configurations, respectively.

Figure 13: (Color online) Snapshot of the pressure fluctuations in the post-
buckling configuration at $f=f_{PC}=5630$ Hz. The larger wavelength
fluctuations propagate at the speed of sound. Figure 14: (Color online)
Snapshot of the pressure fluctuations in the post-contact configuration at
$f=f_{PC}=6250$ Hz. The larger wavelength fluctuations propagate at the speed
of sound.

Both configurations show aerodynamic and acoustic pressure fluctuations,
indicating that acoustic waves are generated at the frequency corresponding to
the flow shear-layer fluctuations. The short wavelength pressure features
downstream of the constriction have a clear aerodynamic origin. Indeed, since
Eq. 9 shows the time evolution of one frequency component of the signal, it is
possible to estimate the phase velocity $c_{f}$, as $\lambda f=c_{f}$, where
$\lambda$ is the wavelength of the feature. At the same time, it is possible
to identify pressure fluctuations characterized by waves with larger
wavelengths moving at the speed of sound. These waves are most noticeable
upstream of the constriction, where the aerodynamic fluctuations are absent.
Hence, it is possible to conclude that the fluid dynamics features described
in Sec. 4, characterized by frequencies $f_{PB}$ and $f_{PC}$, induce the
onset of acoustic waves within the domain.

It is possible to further determine the nature of these acoustic waves by
analyzing the cut-off frequency $f_{n}^{c}$ of the system, i.e., the minimum
frequency that allows the propagation of higher-order modes of order $n$ in
the collapsible tube. The corresponding formula is [15]

$f_{n}^{c}=\frac{c_{0}k^{\perp}_{n}}{2\pi}\sqrt{1-\text{Ma}^{2}}\,,$ (10)

where $k^{\perp}_{n}$ represents the $n$-th eigenvalues of the Helmholtz
equation obtained by assuming a harmonic time-dependent solution in the
homogeneous wave equation. For a circular duct, the first eigenvalue is
$k_{1}^{\perp}=1.841/r$, where $r=0.003$ m is the radius of the non-deformed
part of the collapsible tube, where the wave propagation happens.
Consequently, the corresponding cut-off frequency is of the order of $30$ kHz,
which implies that only acoustic plane waves can propagate in the domain for
the frequency range considered in this study.

##  7 Discussion and conclusions

The analysis of the flow-induced sound generation mechanisms in a collapsible
tube have been performed using CFD and CAA methods. The different flow
features associated to three collapse states of the tube have been identified
and related to the generation of acoustic power. The analysis of the pressure
fluctuation on the walls of the tube and the corresponding sound power
analysis have shown that the maximum of the sound power is produced in the
post-buckling configuration. This suggest the existence of a non-trivial
acoustic tube law, with maximum sound power level located in the range of
pressures corresponding to the non-buckling configuration. Moreover, the
compressible LES model has shown the onset of acoustic plane waves at the
frequencies of to the shear-layer fluctuations. The results obtained in this
analysis open to new perspective research questions.

The relation between the airflow peak frequencies indicated in Figs. 10-12 and
the frequency of possible self-excited oscillations of the airways’ wall, in
turn related to wheezing [8], is not straightforward. Heil et al. [22] have
shown that a collapsible tube in the post-buckling configuration shows self-
excited oscillations for large enough Reynolds numbers. However, understanding
the underlying physical mechanisms for the onset of such oscillations is
currently an open question [25]. A possible research question is how the
energy content of the flow pressure fluctuations due to the unsteady flow
interacting with one side of the collapsible tube in the post-buckling
configuration is coupled with the solid structure. In other words, it is
relevant to determine the frequency response function of a collapsible tube,
and how is it related to the collapsible state of the system. The role of the
flow-induced acoustic waves on the onset of self-excited oscillations can be
also investigated. The analysis of the flow and acoustic features presented in
this work represents a first step in this direction, and it will be addressed
in a future investigation involving a two-way coupled FSI numerical model.

Another interesting perspective is related to the analysis of the non-
dimensional frequency associated to the peaks in the PSD plots 10-12 and the
physics of the fluid features in the post-buckling and post-contact
configuration of the tube. A more in-depth analysis can be designed to prove
this thesis. The strategy would be to simulate the fluid flow in several
instances of post-buckling and post-contact configurations and to compute the
corresponding Strouhal number. The thesis would be proven if the Strouhal
number is consistently equal to $\tilde{f}_{PB}$ for all the post-buckling
simulations and to $\tilde{f}_{PC}$ for all the post-contact simulations. This
analysis will be the object of a future investigation.

## Conflict of interest

The authors declare no conflict of interest.

## Data availability

The simulation data can be provided under reasonable request.

## Acknowledgments

The project is financed by KTH Engineering Mechanics in the thematic area of
Biomechanics, Health, and Biotechnology (BHB). Laudato is supported by the
Swedish Research Council, grant No. 2022-03032. Zea is partly supported by the
Swedish Research Council, grant No. 2020-04668. The authors acknowledge PRACE
for awarding access to the Fenix Infrastructure resources at CINECA, which are
partially funded by the European Union’s Horizon 2020 research and innovation
programme through the ICEI project under grant agreement No. 800858. The
simulations were partly run on the Swedish National Infrastructure for
Computing (SNIC) resources at the PDC Centre for High Performance Computing
(PDC-HPC).

## References

  * [1] Keith Horsfield and Gordon Cumming  “Morphology of the bronchial tree in man”  In _Journal of Applied Physiology_ 24.3, 1968, pp. 373–383
  * [2] Joseph E Flaherty, Joseph B Keller and SI Rubinow  “Post buckling behavior of elastic tubes and rings with opposite sides in contact”  In _SIAM Journal on Applied Mathematics_ 23.4 SIAM, 1972, pp. 446–455
  * [3] FREDERIC G Hoppin Jr, JM Hughes and JERE Mead  “Axial forces in the bronchial tree”  In _Journal of Applied Physiology_ 42.5, 1977, pp. 773–781
  * [4] Ascher H. Shapiro  “Steady Flow in Collapsible Tubes”  In _Journal of Biomechanical Engineering_ 99.3, 1977, pp. 126–147  DOI: 10.1115/1.3426281
  * [5] William A. Conrad  “A new model of the vocal cords based on a collapsible tube analogy”  In _Medical Research Engineering_ 13, 1980, pp. 7–10
  * [6] Claudio Cancelli and T.. Pedley  “A separated-flow model for collapsible-tube oscillations”  In _Journal of Fluid Mechanics_ 157 Cambridge University Press, 1985, pp. 375–404  DOI: 10.1017/S0022112085002427
  * [7] Ingo R. Titze  “The physics of small‐amplitude oscillation of the vocal folds”  In _The Journal of the Acoustical Society of America_ 83.4, 1988, pp. 1536–1552
  * [8] JAMES B Grotberg and NOAM Gavriely  “Flutter in collapsible tubes: a theoretical model of wheezes”  In _Journal of Applied Physiology_ 66.5, 1989, pp. 2262–2273
  * [9] STEPHEN J Lai-Fook and JOSEPH R Rodarte  “Pleural pressure distribution and its relationship to lung volume and interstitial pressure”  In _Journal of Applied Physiology_ 70.3, 1991, pp. 967–978
  * [10] H Teriö  “Acoustic method for assessment of urethral obstruction: A model study”  In _Medical and Biological Engineering and Computing_ 29 Springer, 1991, pp. 450–456
  * [11] Hans Boden and Mats Åbom  “Modelling of fluid machines as sources of sound in duct and pipe systems”  In _ACTA ACUSTICA-LES ULIS-_ 3 LES EDITIONS DE PHYSIQUE, 1995, pp. 549–560
  * [12] Matthias Heil  “Stokes flow in collapsible tubes: computation and experiment”  In _Journal of Fluid Mechanics_ 353 Cambridge University Press, 1997, pp. 285–312
  * [13] R.. Chan and I.. Titze  “Viscoelastic Shear Properties Of Human Vocal Fold Mucosa: Measurement Methodology and Empirical Results”  In _The Journal of the Acoustical Society of America_ 106, 1999, pp. 2008–2021  DOI: 10.1121/1.427947
  * [14] CD Bertram  “Experimental studies of collapsible tubes”  In _Flow Past Highly Compliant Boundaries and in Collapsible Tubes: Proceedings of the IUTAM Symposium held at the University of Warwick, United Kingdom, 26–30 March 2001_ , 2003, pp. 51–65  Springer
  * [15] Mats Åbom  “An introduction to flow acoustics”  Skolan för teknikvetenskap, Kungliga Tekniska högskolan, 2006
  * [16] Mats Åbom, Sabry Allam and Susann Boij  “Aero-acoustics of flow duct singularities at low mach numbers”  In _12th AIAA/CEAS Aeroacoustics Conference (27th AIAA Aeroacoustics Conference)_ , 2006, pp. 2687
  * [17] Christopher G Lausted et al.  “Maximum static inspiratory and expiratory pressures with different lung volumes”  In _Biomedical engineering online_ 5.1 BioMed Central, 2006, pp. 1–6
  * [18] V Oruç and MÖ Çarp  “An experimental investigation of collapsible tube flows at the onset of self-excited oscillations”  In _WIT Transactions on The Built Environment_ 92 WIT Press, 2007
  * [19] Peter J Barnes, Jeffrey M Drazen, Stephen I Rennard and Neil C Thomson  “Asthma and COPD: Basic Mechanisms and Clinical Management”  Elsevier, 2009  DOI: 10.1016/B978-0-12-374001-4.X0001-6
  * [20] John A Evans and William A Whitelaw  “The assessment of maximal respiratory mouth pressures in adults”  In _Respiratory care_ 54.10 Respiratory Care, 2009, pp. 1348–1359
  * [21] James E Gern  “The ABCs of rhinoviruses, wheezing, and asthma”  In _Journal of Virology_ 84.15 Am Soc Microbiol, 2010, pp. 7418–7426
  * [22] Matthias Heil and Jonathan Boyle  “Self-excited oscillations in three-dimensional collapsible tubes: simulating their onset and large-amplitude oscillations”  In _Journal of Fluid Mechanics_ 652 Cambridge University Press, 2010, pp. 405–426
  * [23] Robert J Whittaker, Matthias Heil, Oliver E Jensen and Sarah L Waters  “A rational derivation of a tube law from shell theory”  In _Quarterly Journal of Mechanics and Applied Mathematics_ 63.4 OUP, 2010, pp. 465–496
  * [24] Robert J Whittaker, Matthias Heil, Oliver E Jensen and Sarah L Waters  “Predicting the onset of high-frequency self-excited oscillations in elastic-walled tubes”  In _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 466.2124 The Royal Society Publishing, 2010, pp. 3635–3657
  * [25] Matthias Heil and Andrew L Hazel  “Fluid-structure interaction in internal physiological flows”  In _Annual review of fluid mechanics_ 43 Annual Reviews, 2011, pp. 141–162
  * [26] Hai-Chao Han et al.  “Artery buckling: new phenotypes, models, and applications”  In _Annals of Biomedical Engineering_ 41 Springer, 2013, pp. 1399–1410
  * [27] Goutham Mylavarapu et al.  “Planning human upper airway surgery using computational fluid dynamics”  In _Journal of Biomechanics_ 46.12 Elsevier, 2013, pp. 1979–1986
  * [28] Annunziato Siviglia and Marco Toffolon  “Steady analysis of transcritical flows in collapsible tubes with discontinuous mechanical properties: implications for arteries and veins”  In _Journal of fluid mechanics_ 736 Cambridge University Press, 2013, pp. 195–215
  * [29] Abraham Bohadana, Gabriel Izbicki and Steve S Kraman  “Fundamentals of lung auscultation”  In _New England Journal of Medicine_ 370.8 Mass Medical Soc, 2014, pp. 744–751
  * [30] Semra İçer and Şerife Gengeç  “Classification and analysis of non-stationary characteristics of crackle and rhonchus lung adventitious sounds”  In _Digital Signal Processing_ 28 Elsevier, 2014, pp. 18–27
  * [31] Pavel Kozlovsky, Uri Zaretsky, Ariel J Jaffa and David Elad  “General tube law for collapsible thin and thick-wall tubes”  In _Journal of Biomechanics_ 47.10 Elsevier, 2014, pp. 2378–2384
  * [32] Andrew Wellman et al.  “Test of the Starling resistor model in the human upper airway during sleep”  In _Journal of Applied Physiology_ 117.12 American Physiological Society Bethesda, MD, 2014, pp. 1478–1485
  * [33] Emma Alenius, Mats Åbom and Laszlo Fuchs  “Large eddy simulations of acoustic-flow interaction at an orifice plate”  In _Journal of Sound and Vibration_ 345 Elsevier, 2015, pp. 162–177
  * [34] AL Gregory, Anurag Agarwal and Joan Lasenby  “Collapse of Flexible Tubes-The Tube Law”, 2017
  * [35] T Hassan et al.  “A System for Measuring Sound Transmission Through Joints”  In _2019 IEEE Signal Processing in Medicine and Biology Symposium (SPMB)_ , 2019, pp. 1–4  IEEE
  * [36] Mohamed Yacin Sikkandar, Natteri M Sudharsan, S Sabarunisha Begum and EYK Ng  “Computational Fluid Dynamics: A Technique to Solve Complex Biomedical Engineering Problems-A Review”  In _WSEAS Transactions on Biology and Biomedicine_ 16, 2019, pp. 121–137
  * [37] Harish Palnitkar et al.  “Sound transmission in human thorax through airway insonification: an experimental and computational study with diagnostic applications”  In _Medical & Biological Engineering & Computing_ 58 Springer, 2020, pp. 2239–2258
  * [38] Lukas Schickhofer and Mihai Mihaescu  “Analysis of the aerodynamic sound of speech through static vocal tract models of various glottal shapes”  In _Journal of biomechanics_ 99 Elsevier, 2020, pp. 109484
  * [39] Stefan Schoder et al.  “Hybrid aeroacoustic approach for the efficient numerical simulation of human phonation”  In _The Journal of the Acoustical Society of America_ 147.2 AIP Publishing, 2020, pp. 1179–1194
  * [40] Richard F Spaide  “Choroidal blood flow: review and potential explanation for the choroidal venous anatomy including the vortex vein system”  In _Retina_ 40.10 LWW, 2020, pp. 1851–1864
  * [41] Stefano S Turzi  “Landau-like theory for buckling phenomena and its application to the elastica hypoarealis”  In _Nonlinearity_ 33.12 IOP Publishing, 2020, pp. 7114
  * [42] AL Gregory, A Agarwal and J Lasenby  “An experimental investigation to model wheezing in lungs”  In _Royal Society Open Science_ 8.2 The Royal Society, 2021, pp. 201951
  * [43] Zheng Li et al.  “A one-dimensional flow model enhanced by machine learning for simulation of vocal fold vibration”  In _The Journal of the Acoustical Society of America_ 149.3, 2021, pp. 1712–1723
  * [44] André Miazaki da Costa Tourinho et al.  “A study on tracheoesophageal phonation based on a collapsible channel model”  In _The Journal of the Acoustical Society of America_ 149.3, 2021, pp. 1979–1988
  * [45] M Amin F Zarandi et al.  “Effect of tube length on the buckling pressure of collapsible tubes”  In _Computers in Biology and Medicine_ 136 Elsevier, 2021, pp. 104693
  * [46] Keiichi Itatani et al.  “Hemodynamic parameters for cardiovascular system in 4D flow MRI: Mathematical definition and clinical applications”  In _Magnetic Resonance in Medical Sciences_ 21.2 Japanese Society for Magnetic Resonance in Medicine, 2022, pp. 380–399
  * [47] Elahe Kazemeini et al.  “Critical to know Pcrit: a review on pharyngeal critical closing pressure in obstructive sleep apnea”  In _Frontiers in Neurology_ 13 Frontiers Media SA, 2022, pp. 775709
  * [48] Kamlesh Kumar and Deepu Prabhakaran  “Dynamics of a collapsible tube with internal constriction”  In _Physics of Fluids_ 34.12 AIP Publishing, 2022
  * [49] Jan Liu, Maximilian Heumann, Kent W Stewart and Peter P Pott  “Development of a Venous Collapse Prevention Device for Blood Draw: Initial Proof of Concept”  In _Current Directions in Biomedical Engineering_ 8.2, 2022, pp. 213–216  De Gruyter
  * [50] Anneloes G Munneke, Joost Lumens, Theo Arts and Tammo Delhaas  “A closed-loop modeling framework for cardiac-to-coronary coupling”  In _Frontiers in Physiology_ 13 Frontiers, 2022, pp. 183
  * [51] Enrico Babilio, Ida Mascolo and Federico Guarracino  “From static buckling to nonlinear dynamics of circular rings”  In _Frontiers in Applied Mathematics and Statistics_ 9 Frontiers, 2023, pp. 1115227
  * [52] A-L Dubessy and I Arnulf  “Sleepiness in neurological disorders”  In _Revue Neurologique_ Elsevier, 2023
  * [53] Florian Kraxberger et al.  “On the Alignment of Acoustic and Coupled Mechanic-Acoustic Eigenmodes in Phonation by Supraglottal Duct Variations”  In _Bioengineering_ MDPI AG, 2023
  * [54] Martin Lasota et al.  “Anisotropic minimum dissipation subgrid-scale model in hybrid aeroacoustic simulations of human phonation”  In _The Journal of the Acoustical Society of America_ 153.2 AIP Publishing, 2023, pp. 1052–1063
  * [55] Marco Laudato and Mihai Mihaescu  “Analysis of the contact critical pressure of collapsible tubes for biomedical applications”  In _Continuum Mechanics and Thermodynamics_ Springer, 2023, pp. 1–12
  * [56] Marco Laudato, Roberto Mosca and Mihai Mihaescu  “Buckling critical pressures in collapsible tubes relevant for biomedical flows”  In _Scientific Reports_ 13.1 Nature Publishing Group UK London, 2023, pp. 9298
  * [57] Peng Li, Marco Laudato and Mihai Mihaescu  “Time-Dependent Fluid-Structure Interaction Simulations of a Simplified Human Soft Palate”  In _Bioengineering_ 10.11 MDPI, 2023, pp. 1313
  * [58] Nima Mobadersany et al.  “Pulse wave imaging of a stenotic artery model with plaque constituents of different stiffnesses: Experimental demonstration in phantoms and fluid-structure interaction simulation”  In _Journal of Biomechanics_ 149 Elsevier, 2023, pp. 111502
  * [59] Muhammad Shahzad Shabbir, Nasir Ali and Zaheer Abbas  “Impact of unsteadiness on the non-Newtonian flow of blood in a vascular tube with stenosis and aneurysm: analytical solution”  In _Waves in Random and Complex Media_ Taylor & Francis, 2023, pp. 1–15
  * [60] Sen Zhang and HaoFei Liu  “A fully-coupled three-dimensional fluid-structure interaction study on the externally-pressurized collapsible tube and the internal flow”  In _Journal of Theoretical and Applied Mechanics_ 61.2 Polskie Towarzystwo Mechaniki Teoretycznej i Stosowanej (PTMTS), 2023, pp. 395–406


# Soft finger dynamic stability and slip by Coulomb friction and bulk
stiffness

Hun Jang1∗, Valentyn Petrichenko1∗, Joonbum Bae2, Kevin Haninger1 ∗ denotes
equal contribution1Department of Automation, Fraunhofer IPK, Berlin, Germany
<EMAIL_ADDRESS>Department of Robotics and Mechatronics
Engineering, UNIST, Ulsan, Korea<EMAIL_ADDRESS>
###### Abstract

Soft robotic fingers can safely grasp fragile or non-uniform objects, but
their force capacity is limited, especially with less contact area: objects
which are smaller, not round, or where an enclosing grasp is not feasible. To
improve force capacity, this paper considers two types of grip failure, slip
and dynamic rotational stability. For slip, a Coulomb model for soft fingers
based on total normal and tangential force is validated, identifying the
effect of contact area, pressure, and grip position on effective Coulomb
coefficient, normal force and transverse stiffness. For rotational stability,
bulk stiffness of the fingers is used to develop conditions for dynamic
stability about the initial grasp, and a condition for when the rotation leads
to slip. Together, these models suggest contact area improves grip by
increasing transverse stiffness and normal force. The models are validated in
a range of grasp conditions, shown to predict the influence of object radius
and finger distance on grip stability limits.

##  1 Introduction

Soft fingers aim to improve generalization of grasping and manipulation,
especially with objects which are sensitive or have variation in surface
geometry. Soft fingers have been shown to grasp a large range of objects,
especially spherical or cylindrical objects when an enclosing grasp is used
[1, 2]. Such applications are well-suited to pick-and-place tasks in
unobstructed environments [3], where a large contact area between finger and
object can be employed.

However, soft fingers typically have reduced force transmission capacity
compared with rigid fingers [4]. The soft materials have limited ability to
build normal force at the contact surface [5], and compliance in the normal
direction makes force closure effects more complex [6]. This force
transmission limitation can be even more significant when contact area is
reduced; when smaller or flat objects are grasped, or when cluttered
environments limit the ability to realize large contact areas [3].

[width=0.3]svg-inkscape/intro_1_svg-tex.pdf_tex

(a) Instability with test objects

[width=0.55]svg-inkscape/intro_2_svg-tex.pdf_tex

(b) Instability with a mouse

Figure 1: Grasp instability in soft fingers, where (a) shows initial grasp,
yawing instability, and rolling instability from top to bottom. In (b), the
instability of a larger object is seen.

This can be a limitation in both the ability to manipulate heavy objects [7]
and in moving from pick-and-place to contact-rich manipulation, which often
requires more force [8]. Ongoing work has explored increasing force capacity
using stiffening elements [9], jamming techniques to increase stiffness [10],
incompressible actuating fluids [7], additional kinematic structure to
increase contact area [11], and structured compliance which is stiffer where
forces are needed [12]. The force capacity is also affected by grip parameters
such as degree of actuation, robot grip pose, and finger distance. However,
the optimal grip parameters are highly object specific. To determine good grip
parameters for a specific object, deep learning approaches are common for
grasp pose planning [13], but not yet applied to optimizing force capacity in
soft grasps.

Soft finger force capacity is limited by at least two phenomena, slip and
stability. Grip stability has a wide range of definitions [14, 15]. Planar
grip stability has been studied for force-closure [16], form-closure [17], and
stiff-but-underactuated fingers [18]. However, these analyses assume point
contact between object and gripper, as well as the gripper contact point being
perfectly stiff, allowing the expression of a static grip matrix $G$ which
forms the basis of most further analysis [14]. These analyses cannot be
directly applied to soft fingers, where contact is neither point-like and the
contact patches move in response to force. Some approaches do model finger
compliance [15], such as for tendon-driven fingers [19] or spherical
fingertips [6], considering Coulomb force conditions with linear stiffness
models. However, dynamic instability can be observed in soft fingers, where
rotation away from an initial object orientation occurs, as seen in Figure 1.
This effect is not initiated by slip, and requires modeling motion transverse
to the grip normal.

Slip is typically modelled with Coulomb friction in grasping [14, 20].
However, a Coulomb friction model is a point-contact model, raising questions
about its applicability when the contact area is larger. Finite-element based
contact models are being developed [21, 22], but not yet applied to grip slip
or stability. For spherical soft fingers, closed-form friction models which
include a torsional friction about the contact normal have been proposed [23,
24, 6], but not extended to general soft fingers. Deep-learning approaches for
slip detection or prediction [20] have been used for trajectory planning [25],
but not finger design or optimization of grip parameters.

To increase the force capacity of soft fingers, this paper proposes a model
for dynamic stability and validates the bulk Coulomb model. Compared with
friction modelling in soft fingers [20, 24], we examine the effect of contact
area, pressure, and offset, finding their impact on effective Coulomb
coefficient. This suggests increasing normal force is the primary way to
increase force capacity. However, we find that dynamic instability occurs as
grip force increases. This is then modelled, considering bulk finger
stiffness, developing dynamic stability conditions which depend on grip force,
transverse stiffness, and object radius. Compared with classical grip
stability [14, 16, 17, 18, 19], this is the dynamic stability of a compliant
system, analyzing the convergence to a certain pose. This paper first
validates the Coulomb model for pneumatic fingers, finding the influence of
contact area on finger slip force. Next, a model for dynamic rotational
stability is proposed, and conditions for this rotation to cause slip
developed. These models are validated experimentally on a range of object
sizes and grip conditions, validating the effect of object radius, contact
area and finger stiffness on maximum stable grip force.

##  2 Friction modelling

This section introduces the friction model, validating the bulk Coulomb model
for soft fingers. Coulomb friction models slip as a violation of the condition
$|f_{t}|<\mu f_{n}$, where $\mu\geq 0$ is the Coulomb friction coefficient,
$f_{n}$ normal force and $f_{t}$ tangential force. These models are popular,
even in soft robotics [26, 9, 27]. They can also be made stochastic, where
$\mu$ can be considered to vary with a distribution [20].

However, the Coulomb model is a point model, and it is unclear how it applies
to larger contact areas. One reason for this is that the pressure distribution
over the contact area not uniform, which can result in local slip, where local
relative motion occurs before bulk slip [23]. Computational methods to
estimate contact pressure have been proposed [22], but not yet extended with
frictional models. An additional effect on larger contact areas is torsional
friction, which can be modelled on spherical geometries [24], but not yet
extended to general contact. To investigate the accuracy of the bulk Coulomb
friction over soft grasp parameters (contact area, offset, pressure,
material), we develop a test environment as seen in Figure 2.

Figure 2: Experimental set up for the friction experiments.

###  2.1 Validation of bulk Coulomb

The experimental setup seen in Figure 2 is used with test objects 3D printed
from Tough PLA which are sections of a cylinder of radius $30$ mm. In this
way, the contact area is controlled by the test object, provided the finger is
in contact with the full arc of the object, which is verified visually per
experiment. The test objects are mounted onto a Force/Torque sensor (ME-
Meßsysteme MP11, load limit 500N, 20Nm, sampling rate 125Hz) to allow
measurement of the total normal and transverse force.

Pneumatic fingers are made with Wacker RT625 Elastosil, and the kinematic
structure is an adjusted form of the Pneunet fingers [28]. They are actuated
by controlled pressure valves (VEAB-L-26-D9-Q4-V1-1R1, Festo), where a
constant desired pressure is applied and real pressure measured. The fingers
are mounted on a Universal Robots UR10, which is used to initiate contact and
load the fingers. We then compare five different parameters: pressure, contact
area, relative position between finger and contact object (horizontal and
vertical offset) and contact surface material. The range of parameters tested
is shown in Table 1. The high friction condition applies a rubber to the
surface of the test object.

For each set of parameters, the following experiment process is followed: (i)
the finger is actuated in free space, (ii) the robot moves the finger into
contact, (iii) contact along the entire test object is visually verified, (iv)
linear motion in the positive y direction begins, with a velocity of $2$mm/s
and to a distance of $30$mm, (v) the robot moves back to the initial pose. As
the robot moves in the y direction, the total transverse force, which is
friction force, can be directly measured as $f_{t}=f_{y}$ from the
force/torque sensor. The total normal force is calculated as
$f_{n}=\sqrt{f_{x}^{2}+f_{z}^{2}}$, where $f_{z}$ is assumed to be from
asymmetrical contact pressure, not friction, as the contact has been
established with motion in x. The ratio of friction force to the total normal
force is plotted, which saturates at $\mu$ during slip. Experiments are
conducted for each parameter combination three times, and the graph displays
the mean value.

Parameter name | Value
---|---
Pressure [$bar$] | 0.4, 0.8, 1.2
Contact area [$cm^{2}$] | 2.1(S), 4.2(M), 8.4(L)
Horizontal offset [$mm$] | 0, -10, -20
Vertical offset [$mm$] | 40, 60, 80
Contact surface material friction | Ordinary, High
Table 1: Parameters of soft finger experiment for slip

[width=0.48]svg-inkscape/14Sept_pressure_G_svg-tex.pdf_tex

(a) Pressure

[width=0.48]svg-inkscape/14Sept_CA_G_svg-tex.pdf_tex

(b) Contact area

[width=0.48]svg-inkscape/14Sept_HOR_G_svg-tex.pdf_tex

(c) Horizontal offset

[width=0.48]svg-inkscape/14Sept_VER_G_svg-tex.pdf_tex

(d) Vertical offset

Figure 3: Experimental results for normal force

[width=0.48]svg-inkscape/14Sept_Fric_pressure_G_svg-tex.pdf_tex

(a) Pressure

[width=0.48]svg-inkscape/14Sept_FRIC_CA_G_svg-tex.pdf_tex

(b) Contact area

[width=0.48]svg-inkscape/14Sept_off_H_G_svg-tex.pdf_tex

(c) Horizontal offset

[width=0.48]svg-inkscape/14Sept_off_V_G_svg-tex.pdf_tex

(d) Vertical offset

Figure 4: Experimental results for friction force

[width=0.46]svg-inkscape/14Sept_Fc_1_G_svg-tex.pdf_tex

(a) Friction coefficient for parameters P: pressure, P(High Fc): high friction
surface, CA: contact area

[width=0.46]svg-inkscape/14Sept_Fc_2_G_svg-tex.pdf_tex

(b) Friction coefficient for green: horizontal offset, brown: vertical offset

Figure 5: Experimental results for friction coefficient

In all of the results for friction force in Figure 3 and ratio of friction to
normal force in Figure 5, three phases can be identified: the stick phase,
where the fingertip is not moving while the robot moves with constant
velocity, a transition region, and a constant friction phase as sliding
occurs. The slope in the stick phase represents the soft finger’s bulk
stiffness in y, $k_{y}$, as it stretches horizontally.

Increased pressure leads to higher normal forces, as seen in Figure 3(a), but
does not affect stiffness in the direction of motion. Similarly, contact area
changes the normal force significantly, but has negligible impact on $\mu$ in
Figure 5(a). However, contact area does increase this stiffness. Neither of
these parameters affects $\mu$ as much as material, shown in gray in 5(a). On
the other hand, variation in the robot’s horizontal and vertical position
cause larger changes in $\mu$, as seen in Figure 5(b), possibly due to the
different contact pressure distributions. Additionally, a higher vertical
offset increases the transverse stiffness as seen in the stick phase of Figure
4, as the object is closer to the finger basis.

The results are summarized in Table 2, showing the effects of a unit change of
each parameter on $f_{n}$ and $k_{y}$ and the range of $\mu$ realized. These
results suggests that normal force and material are more important than
contact area and pressure for the friction force limits. Variation in the grip
pose can change the normal pressure distribution, which can cause larger
variation in $\mu$. Stiffness in the direction of motion is mostly affected by
vertical offset and contact area.

Parameter |  $f_{n}$ [N/unit] |  $k_{y}$ [N/(mm unit)] |  $\mu$ (std)
---|---|---|---
Pressure | $3.69/bar$ | $0.046/bar$ | 0.58 ($0.025$)
Contact area | $0.47/cm^{2}$ | $0.024/cm^{2}$ | 0.63 ($0.017$)
Horizontal offset | $0.14/mm$ | $0.005/mm$ | 0.62 ($0.086$)
Vertical offset | $0.07/mm$ | $0.089/mm$ | 0.72 ($0.058$)
Table 2: Summary of the effect of a unit change in the parameter on normal
force and stiffness in direction of motion $k_{y}$, and the average $\mu$
values reached over the range of parameters tested.

##  3 Dynamic grip stability

If a bulk Coulomb model with a constant $\mu$ is accurate, increasing force
capacity should be as easy as increasing normal force $f_{n}$. However, in
some cases this leads to rotation of the object, in either yawing or rolling
as seen in 1(a), which can lead to slip as seen in Figure 1(b).

To understand why this happens, we define dynamic grip stability as the
convergence of the object rotation $\Theta$ to the initial rotation angle
$\Theta=0\degree$ over time. To evaluate the dynamic stability, the model in
Figure 6 is used. We model the object’s pose relative to the robot flange with
stiffnesses, where one side is fixed to the flange and the other to the
object. We assume that the stiffnesses are linear and equal between the
fingers, ignore any torsional stiffness at the contact area, and assume that
the direction of the normal and transverse force stays constant. We derive a
dynamics equation and analyze the stability about the equilibrium point
$\Theta=0\degree$ using the eigenvalues, where only a rotational motion of the
body with the torsion angle $\Theta$ is studied.

Let the bulk stiffness of the fingers be modeled with the corresponding
stiffness $k_{x}$ in the normal force direction and $k_{y}$ in transversal
direction, with a radius $r$ to the center of mass of the object. Grip force
of finger in x direction is modeled by the spring with a relative spring
displacement $\Delta x$ and preload in the y direction is assumed to be $0$.

Figure 6: View from above of the grasped object, fingers are modeled as
springs

###  3.1 Rotational dynamics

We model the pure rotation of the object seen in Figure 6. The dynamics can be
determined based on the energies of the system with Lagrangian mechanics

$\displaystyle\frac{d}{dt}\frac{\partial\mathcal{L}}{\dot{q}_{i}}-\frac{\partial\mathcal{L}}{q_{i}}=0,$
(1)

where $\mathcal{L}=T-V$. The total kinetic energy of the system $T$ in the
case of a pure rotation is given by

$\displaystyle T=\frac{1}{2}I\dot{\Theta}^{2},$ (2)

where $I$ is the rotational inertia and $\dot{\Theta}$ the rotation speed
about $z$. The potential energy of the system $V$ is derived from the relative
travel of the spring through the preload and through the distance covered by
the springs during object rotation. As the body rotates, the spring
displacement changes, giving the potential energy

$\displaystyle V=k_{x}(\Delta
x-r(1-\cos{\Theta}))^{2}+k_{y}(r\sin{\Theta})^{2},$ (3)

where $r$ is object radius, $k_{x}$ and $k_{y}$ the stiffness in respective
directions, and $\Delta x$ the preload displacement. With (2) and (3) inserted
into (1) the dynamic equation can be calculated to

$\begin{split}I\ddot{\Theta}=&2k_{x}r\sin{\Theta}\Delta
x-2k_{x}r^{2}(1-\cos{\Theta})\sin{\Theta}\\\
&-2k_{y}r^{2}\sin{\Theta}\cos{\Theta}.\end{split}$ (4)

###  3.2 Linearized dynamics and stability

With defining of the state vector $x=[\Theta,\dot{\Theta}]^{T}$, the problem
can be considered from the control systems point of view. The system is
nonlinear, but the linearized dynamics can be found and dynamic stability
about $x=[0,0]^{T}$ can be analyzed. When the system is dynamically stable
about a state, it converges to that state over time from starting conditions
near that state. When the system is dynamically unstable, it will diverge.

With the introduced state variables, the system matrix $A$ linearized about
the equilibrium point $\Theta=0$ can be written as

$\displaystyle A=\begin{bmatrix}0&1\\\ \frac{2r}{I}(k_{x}\Delta
x-k_{y}r)&0\end{bmatrix},$ (5)

where $\dot{x}\approx Ax$ about $\Theta=0\degree$. The $A$ matrix has two
eigenvalues which can be found as

$\begin{split}\lambda_{1/2}&=|\lambda I-A|=\pm\sqrt{\frac{2r}{I}(k_{x}\Delta
x-k_{y}r)}\\\ &=\begin{cases}\pm\sqrt{\frac{2r}{I}(k_{x}\Delta
x-k_{y}r)},\text{ for }k_{x}\Delta x>k_{y}r\\\ \pm
i\sqrt{\frac{2r}{I}(k_{y}r-k_{x}\Delta x)},\text{ for }k_{x}\Delta
x<k_{y}r\end{cases}.\end{split}$ (6)

If the term $k_{x}\Delta x=f_{p}$ exerted by the finger, is greater than the
stiffness in the transverse direction $k_{y}$ times the distance to the center
of masses $r$, then there will be a pole located in the right half-plane,
indicating instability. This gives a dynamic stability condition about the
equilibrium point $\Theta=0\degree$ of

$\displaystyle\underbrace{k_{x}\Delta x}_{f_{p}}<k_{y}r,$ (7)

and we denote the force when this condition is violated as $f_{p}^{i}=k_{y}r$.
As seen in Figure 1, instability can also happen about the x axis, in this
case the model results in the condition $f_{p}^{i}=k_{z}r$.

###  3.3 Rest position

The derived dynamic equation (4) can be used to predict the angle that occurs
after the force is applied by the fingers by setting $\ddot{\theta}$ to zero,
which has solutions for $\Theta_{r}=n\pi$ for $n\in\mathbb{Z}$, and additional
solutions at $\Theta_{r}\neq n\pi$ when

$\left|\frac{k_{x}\Delta x-k_{x}r}{k_{y}r-k_{x}r}\right|<1.$ (8)

The progression of the rest position angle $\Theta_{r}$ over the applied force
is shown for selected system parameters in the Figure 7. It can be seen that
the object rotates when a certain limit force is exceeded, where this limit
force follows from (7). It should be noted that if the condition (8) for the
analytical solution is not satisfied, the progression of the rest position
angle cannot be calculated and only the force $f^{i}_{p}$ can be determined.

The comparison of these results for different $k_{x}$ and $k_{y}$ values can
be seen in Figure 8. As indicated by (7), the initiation of rotation depends
on $k_{y}$, increasing as $k_{y}$ increases. However, the slope after rotation
is initiated decreases as $k_{x}$ increases.

[width=0.85]svg-inkscape/theta_svg-tex.pdf_tex

Figure 7: Rest position $\Theta_{r}$ over the changing of preload force
according to (4) and maximal friction angle $\Theta_{f}$ according to (12)
with $k_{x}=0.3N/mm$, $k_{y}=0.066N/mm$,$r=35mm$ and $\mu=0.6$.

[width=0.8]svg-inkscape/kxky_svg-tex.pdf_tex

Figure 8: Rest angle $\Theta_{f}$ progression as $k_{x}$ and $k_{y}$ vary.

###  3.4 Rotation leading to slip

As the object rotates, the normal and transverse forces change, possibly
leading to slip. The Figure 6 shows the normal force $f_{n}$ and static
friction force $f_{t}$, which from the nonlinear model can be defined as
follows:

$\displaystyle f_{n}$ $\displaystyle=$ $\displaystyle k_{x}(\Delta
x-r(1-\cos{\Theta}))$ (9) $\displaystyle=$ $\displaystyle
f_{p}-k_{x}r(1-\cos{\Theta})$ (10) $\displaystyle f_{t}$ $\displaystyle=$
$\displaystyle k_{y}r\sin{\Theta},$ (11)

where $f_{p}=k_{x}\Delta x$ is the preload force on the finger. These
equations are connected by the Coulomb friction condition $|f_{t}|\leq
f_{n}\mu$. With the assumption that the direction of the forces $f_{n}$ and
$f_{t}$ doesn’t change with rotation, the limiting case for the angle from
which the slip begins can be determined by solving $f_{t}=\mu f_{n}$, which
gives [29]

$\Theta_{f}=2(\arctan(\frac{a\pm\sqrt{a^{2}+b^{2}-c^{2}}}{b+c})+n\pi),n\in\mathbb{Z},$
(12)

where $a=k_{y}r,b=-\mu k_{x}r$ and $c=\mu f_{p}-\mu k_{x}r$. The angle
$\Theta_{f}$ indicates the rotation at which slip occurs. This can be
calculated for the same set of parameters and is also plotted in the Figure 7.
It follows that slip does not occur if the rest position of the system
$\Theta_{r}$ remains smaller than the critical angle for the friction
condition $\Theta_{f}$.

##  4 Experimental validation of grip stability model

To validate the model with stability conditions, silicon-based fingers are
used in a variety of conditions. The model in Section 3 makes two predictions:
when rotation away from initial orientation begins $f^{i}_{p}(k_{x},k_{y},r)$,
and when the rotation leads to slip at $f^{s}_{p}(k_{x},k_{y},r)$. For some
grip conditions, these two are close - when the slope of $\Theta_{r}$ over
$f_{p}$ is high. In other cases, $\Theta_{r}$ increases slowly.

###  4.1 Experimental set up

The first goal is to test the stability condition (7) about the equilibrium
point $\Theta=0\degree$ when grasping an object. The assumptions in the
derivation were that the distance between the fingers and the object center
$r$ does not change and the fingers have the same properties, i.e. the
stiffnesses $k_{x}$, $k_{y}$, $k_{z}$ are matched in both fingers.
Furthermore, it is assumed that the stiffness of the fingers in the transverse
direction $k_{y}$ does not change with changing in pressure, as indicated in
Section 2.

The model predicts that if the force applied by the fingers is less than the
term $k_{y}r$, then the object will remain in its initial rest position
$\Theta=0\degree$. If this limit force is exceeded, the rotation of the body
should happen, similar to the Figure 7 and a new rest position $\Theta_{r}$
should be assumed. As pressure increases, the slip force $f^{s}_{p,x}$ will be
exceeded which often results in smaller objects flying from the fingers.

The stiffness of the finger is needed to validate the stability model. In
order to calculate this, the test object is fixed on the force sensor and with
the measuring force and displacement of robot’s end-effector, the stiffness of
finger can be determined. First, the grip conditions from the stability
experiment are recreated - relative pose between object and gripper, an the
pressure in the finger. The end-effector is moved from $x_{0}$ to $x_{1}$, the
forces measured at each one as $f_{0}$ and $f_{1}$, and the stiffness $k$ can
be calculated with a simple relationship $k=\frac{f_{1}-f_{0}}{x_{1}-x_{0}}$.
To measure the stiffness in the y- and z- direction, a Cartesian motion of
$1-3$ mm is used to avoid slip and stay in the linear range.

###  4.2 Force instability experiments

In each set of experiments, the object is grasped from an experiment table,
with an initial grasp location flush with the bottom of the fingers. The
pressure is increased in steps of around $0.15$ bar, waiting for $10$ seconds
and stopping the experiment when instability and slip has occurred.

####  4.2.1 Validation of force stability limit

With the friction coefficient determined from the Section 2 and the calculated
stiffnesses from the description above, the determined and observed values can
be summarized in Table 3. The object and experimental setup can be seen in
Figure 1(a).

Parameter | Value
---|---
Stiffness $k_{x}[N/mm]$ | 0.14
Stiffness $k_{y}[N/mm]$ | 0.11
Stiffness $k_{z}[N/mm]$ | 0.05
Object radius $r[mm]$ | 15
Expected instability preload force about z $f^{i}_{p}[N]$ | 1.65
Expected instability preload force about x $f^{i}_{p}[N]$ | 0.75
Observed slip force about x $[N]$ | 1.5
Table 3: Measured parameters for stability experiments, where $\mu=0.6$.

To avoid rotation, the finger normal force should not exceed $k_{y}r=1.65N$
for instability about z, and not greater than $k_{z}r=0.75N$ for instability
about x. In the experiments, the rotation about x was observed at force
$f^{s}_{p}=1.5N$. And this also should be expected from the model, since the
force limit about x is smaller than about z.

####  4.2.2 Finger distance and force limit

The effective finger stiffness is varied by varying the finger distance $d$
via an electrical parallel gripper, i.e. the offset between the finger and the
body at each measurement, which influences the stiffnesses and preload force.
A test object with a length of $30$ mm is grasped.

Table 4 shows a comparison of the derived model and the measured values of
slip force $f^{s}_{p,x}$ about x at the different distance $d$ between the
fingers. Furthermore, the expected instability force $f^{i}_{p,x}$ and
$f^{i}_{p,z}$ is written in the table at which the rotation should begin,
about the x and z axis, respectively.

$d$ [mm] |  $f^{s}_{p,x}$ [N] |  $f^{i}_{p,z}$ [N] |  $f^{i}_{p,x}$ [N] |  $k_{x},k_{y},k_{z}$ [N/mm]
---|---|---|---|---
50 | 1.5 | 1.65 | 0.75 | 0.14, 0.11, 0.05
40 | 2.15 | 1.2 | 1.05 | 0.11, 0.08, 0.07
30 | 2.9 | 1.2 | 6 | 0.1, 0.08, 0.4
20 | 3.38 | 1.2 | 9.3 | 0.1, 0.08, 0.62
Table 4: Validation of the model, where $f^{s}_{p,x}$ is the observed slip
force limit.

(a) Selection of possible contact objects

[width=0.8]svg-inkscape/val_2_svg-tex.pdf_tex

(b) Instability limit at different finger distances

Figure 9: Contact instability over variation in object properties, all
observed instabilities here were rolling

It can be seen that when the stiffness in z is increased, the force that leads
to dynamic instability and thus to slip of object increases, as predicted by
the model. While the trend is correct, at higher $k_{z}$ values, the observed
maximum force is lower than the model prediction. We attribute this to the
non-negligible $\Delta y$ in realistic gripping scenarios, which gives
additional transverse forces which can violate the slip condition earlier.

####  4.2.3 Changes in object properties

The grip stability limits are next assessed over a range of objects with
different lengths and contact areas. The objects and the resulting observed
slip force can be seen in Figure 9(a) and (b), respectively.

We can see that the object with a larger $r$ results in a higher $f^{s}_{p}$,
as predicted by the instability model. Additionally, the object with the
larger contact area allows a higher grasp force. The larger contact area was
shown to increase $k_{x}$ in Table 2, which would result in a higher
instability force $f_{p}^{i}$.

###  4.3 Rest angle validation

[width=0.75]svg-inkscape/theta1_svg-tex.pdf_tex

Figure 10: Experiment result for rest position with $k_{y}=0.05N/mm$, $r=25mm$

To validate the rest angle that the object takes as grip force is increased,
we grasp a test object as seen in Figure 2(a) of length $50$ mm. A long
indicating stick is glued to the bottom, and the object is lightly lifted off
the table and the finger pressure is increased. As the pressure increases, the
robot’s last axis is jogged back to align the indicator, providing a
measurement of the rotation angle. The results can be seen in Figure 10, where
the measured rest angle can be seen. While the initiation of rotation at the
instability force $f^{i}_{p}$ matches the model prediction well, the resulting
angles of the object are lower than predicted by the model by at least a
factor of $2$. We attribute this to the linearity assumptions of the finger
stiffness, where additional displacement of the fingers may change the
stiffness properties.

##  5 Conclusion

This paper investigated models to increase the friction force capacity of soft
fingers. It was shown that with a bulk Coulomb model, increasing contact area
mostly increases normal force and transverse stiffness, where coefficient
$\mu$ is mostly unchanged. The normal force was then shown to have a
rotational stability limit, where above a threshold grip force $f_{p}^{i}$,
rotation away from initial grasp point occurs. The instability model was
validated for a range of grip parameters and objects, showing increasing
object radius and transverse stiffness increase the stable normal force limit.

While the trends in object radius and bulk finger stiffness are validated, the
observed instability force is lower than anticipated at higher stiffness
values. Future work can apply more complex finger models which more accurately
describe the gripping force direction and nonlinear stiffnesses.

## References

  * [1] R. Ozawa and K. Tahara, “Grasp and dexterous manipulation of multi-fingered robotic hands: A review from a control view point,” _Advanced Robotics_ , vol. 31, no. 19-20, pp. 1030–1050, 2017.
  * [2] J. Hernandez, M. S. H. Sunny, J. Sanjuan, I. Rulik, M. I. I. Zarif, S. I. Ahamed, H. U. Ahmed, and M. H. Rahman, “Current Designs of Robotic Arm Grippers: A Comprehensive Systematic Review,” _Robotics_ , vol. 12, no. 1, p. 5, Jan. 2023.
  * [3] C. Piazza, G. Grioli, M. Catalano, and A. Bicchi, “A Century of Robotic Hands,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 2, no. 1, pp. 1–32, May 2019.
  * [4] M. Li, A. Pal, A. Aghakhani, A. Pena-Francesch, and M. Sitti, “Soft actuators for real-world applications,” _Nature Reviews Materials_ , vol. 7, no. 3, pp. 235–249, Nov. 2021.
  * [5] S. Liu, F. Wang, Z. Liu, W. Zhang, Y. Tian, and D. Zhang, “A Two-Finger Soft-Robotic Gripper With Enveloping and Pinching Grasping Modes,” _IEEE/ASME Transactions on Mechatronics_ , vol. 26, no. 1, pp. 146–155, Feb. 2021.
  * [6] H. Dong, C. Qiu, D. K. Prasad, Y. Pan, J. Dai, and I.-M. Chen, “Enabling grasp action: Generalized quality evaluation of grasp stability via contact stiffness from contact mechanics insight,” _Mechanism and Machine Theory_ , vol. 134, pp. 625–644, 2019.
  * [7] O. Azami, K. Ishibashi, M. Komagata, and K. Yamamoto, “Development of hydraulically-driven soft hand for handling heavy vegetables and its experimental evaluation,” in _2023 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2023, pp. 2577–2583.
  * [8] A. Billard and D. Kragic, “Trends and challenges in robot manipulation,” _Science_ , vol. 364, no. 6446, p. eaat8414, Jun. 2019.
  * [9] W. Park, S. Seo, J. Oh, and J. Bae, “A sensorized hybrid gripper to evaluate a grasping quality based on a largest minimum wrench,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 3243–3250, 2020.
  * [10] K. Singh, S. Gupta, A. Khosla, and H. Furukawa, “Transforming Soft Robotics: Laminar Jammers Unlocking Adaptive Stiffness Potential in Pneunet Actuators,” _ECS Journal of Solid State Science and Technology_ , vol. 12, no. 4, p. 047007, 2023.
  * [11] J. Zhou, S. Chen, and Z. Wang, “A Soft-Robotic Gripper With Enhanced Object Adaptation and Grasping Reliability,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 4, pp. 2287–2293, Oct. 2017.
  * [12] R. M. Hartisch and K. Haninger, “Compliant finray-effect gripper for high-speed robotic assembly of electrical components,” _arXiv preprint arXiv:2301.08431_ , 2023.
  * [13] R. Newbury, M. Gu, L. Chumbley, A. Mousavian, C. Eppner, J. Leitner, J. Bohg, A. Morales, T. Asfour, D. Kragic, D. Fox, and A. Cosgun, “Deep Learning Approaches to Grasp Synthesis: A Review,” _IEEE Transactions on Robotics_ , pp. 1–22, 2023.
  * [14] M. A. Roa and R. Suárez, “Grasp quality measures: Review and performance,” _Autonomous Robots_ , vol. 38, no. 1, pp. 65–88, Jan. 2015\.
  * [15] D. Prattichizzo, J. C. Trinkle, and B. Siciliano, “Grasping.” 2008.
  * [16] J. C. Trinkle, A. O. Farahat, and P. F. Stiller, “Second-order stability cells of a frictionless rigid body grasped by rigid fingers,” in _Proceedings of the 1994 IEEE International Conference on Robotics and Automation_. IEEE, 1994, pp. 2815–2821.
  * [17] A. Bicchi and V. Kumar, “Robotic grasping and contact: A review,” in _Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No.00CH37065)_ , vol. 1, Apr. 2000, pp. 348–353 vol.1.
  * [18] L. Birglen, T. Laliberté, and C. M. Gosselin, _Underactuated Robotic Hands_. Springer, 2007, vol. 40.
  * [19] M. Haas-Heger, G. Iyengar, and M. Ciocarlie, “Passive Reaction Analysis for Grasp Stability,” Jan. 2018.
  * [20] Z. Liu and R. D. Howe, “Beyond Coulomb: Stochastic Friction Models for Practical Grasping and Manipulation,” _IEEE Robotics and Automation Letters_ , 2023.
  * [21] V. Kurtz and H. Lin, “Contact-Implicit Trajectory Optimization with Hydroelastic Contact and iLQR,” Aug. 2022.
  * [22] J. Masterjohn, D. Guoy, J. Shepherd, and A. Castro, “Velocity Level Approximation of Pressure Field Contact Patches,” Jun. 2022.
  * [23] A. Fakhari, M. Keshmiri, and I. Kao, “Development of realistic pressure distribution and friction limit surface for soft-finger contact interface of robotic hands,” _Journal of Intelligent & Robotic Systems_, vol. 82, pp. 39–50, 2016.
  * [24] A. Fakhari, I. Kao, and M. Keshmiri, “Modeling and control of planar slippage in object manipulation using robotic soft fingers,” _ROBOMECH Journal_ , vol. 6, no. 1, pp. 1–11, 2019.
  * [25] W. Mandil and K. Nazari, “Action conditioned tactile prediction: A case study on slip prediction,” _arXiv preprint arXiv:2205.09430_ , 2022.
  * [26] N. Chavan-Dafle, R. Holladay, and A. Rodriguez, “In-Hand Manipulation via Motion Cones,” Feb. 2019.
  * [27] F. R. Hogan, J. Ballester, S. Dong, and A. Rodriguez, “Tactile Dexterity: Manipulation Primitives with Tactile Feedback,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2020, pp. 8863–8869.
  * [28] P. Polygerinos, Z. Wang, J. T. Overvelde, K. C. Galloway, R. J. Wood, K. Bertoldi, and C. J. Walsh, “Modeling of soft fiber-reinforced bending actuators,” _IEEE Transactions on Robotics_ , vol. 31, no. 3, pp. 778–789, 2015.
  * [29] K. Zelator, “The equation asinx+bcosx=c and a family of cyclic heron quadrilaterals,” p. 3, May 2008.