# Numerical investigation of fluid-structure interaction in a pilot-operated
microfluidic valve

Ahmed Aissa-Berraies  E. Harald van Brummelen  Ferdinando Auricchio

###### Abstract

The present paper is concerned with numerical investigation of the performance
of a pilot-operated control valve based on shape memory alloy actuation
control. The valve under investigation can be integrated into miniaturized
hydraulic systems and is developed to perform precise dispensing, mixing, or
dosing tasks while being able to withstand relatively high pressure
differences. The study evaluates the valve’s response under the current ON/OFF
and the desired proportional control regimes using numerical methods for
fluid-structure interaction. The computational model replicates the operation
of the valve, which requires an understanding of the complex interactions
between the fluid flow with the pressurized valve and the contact with the
valve seat during the opening and closing processes. In addition, the model
leverages advanced numerical techniques to overcome several complexities
arising mainly from the geometrical, material, and contact nonlinearities, and
to mitigate the shortcomings of the partitioned fluid-structure interaction
approach. Several simulations are conducted to examine the valve’s structural
and flow behavior under varying pressure conditions. Results indicate that the
valve is adequate for ON/OFF actuation control but is susceptible to flow-
induced vibrations during the proportional control regime that occurs due to
the sharp pressure drop in the valve-seat gap and the ensuing Venturi effect,
which counteract the opening of the main valve. The fluid-structure-
interaction simulations provide insight into the mechanism underlying the
flow-induced vibrations, which can serve to improve the design and enhance the
performance of the valve in microfluidic applications.

###### keywords:

Pilot-operated valve , Microfluidic system , Valve performance , Fluid-
structure interaction simulation , Contact , Viscous fluid , Hyperelastic
material , ANSYS Multiphysics

††journal: Journal of Fluids and Structures

[inst1]organization=Computational Mechanics & Advanced Materials Group,
Department of Civil Engineering and Architecture, Università degli Studi di
Pavia,addressline=Via Adolfo Ferrata, city=Pavia, postcode=27100,
country=Italy [inst2]organization=Faculty of Mechanical Engineering,
Multiscale Engineering Fluid Dynamics Section,addressline=P.O. Box 513,
city=Eindhoven, postcode=5600 MB, country=The Netherlands

Numerical modeling of a microfluidic pilot-operated valve using partitioned
fluid-structure interaction methods

Solution strategy for the multiphysics problem, taking into account high
added-mass effects and contact events

Valve performance evaluation under ON/OFF and proportional operation regimes

Assessment of flow-induced vibration during partial opening

##  1 Introduction

During the last four decades, the downscaling of piping systems has led to the
miniaturization of mechanical devices like valves and pumps. Extensive
research projects have been conducted to reduce the size of hydraulic devices
to the mesoscale and microscale levels. This trend has enabled the
commercialization of smaller products, such as coffee machines, dosing
systems, and dispensing machines, endowed with efficient and sophisticated
hydraulic circuits. Likewise, the leverage of microfluidic devices allowed for
advanced scientific development in medical research. As an example, it
supported the study of RNA sequencing [1] and intercellular interaction [2],
which has eventually accelerated the development process of various medicines
and drugs, such as the COVID-19 vaccines.

Valves are designed to be integrated into a circuit and control the fluid
flow, the pressure, or the flow direction. Several valve types are used in
various industrial applications since they can execute one or multiple
services, such as starting and stopping the fluid flow, flow throttling,
pressure regulation, relieving the circuit from overpressure, and redirecting
the flow. However, it may become challenging to accommodate multiple valves
with different features as required to deliver a specific flow service,
especially in miniaturized hydraulic systems. Such a limitation has encouraged
the realization of multi-tasking microfluidic valves that can accomplish two
or more functions at the same time. In this context, some devices can start,
stop, and change the flow direction, while others can control the flow and
flow direction and simultaneously regulate the fluid flow rate. Various
miniaturized valves have been successfully developed and commercialized [3].
In particular, the pilot-operated dual valve is one of the prominent types
nowadays. Such a valve can be perceived as modular, i.e., with two stages
mounted together to form a single unit. The pilot stage controls and regulates
the flow, while the main stage controls the direction of the flow. The device
is designed to improve the response time, reduce the power consumption
compared to single-stage valves, and ultimately relax the hydraulic circuit
density [4].

This paper presents a computational investigation of a smart pilot-operated
control valve (POCV), i.e., a two-way medium-separated valve comprising two
‘normally closed’ diaphragm-type stage valves [5]. The POCV is designed to
operate in two modes: The first mode is where the valve operates as an ON/OFF
control valve. That is, the actuator shifts into an open or closed position as
a function of a specific digital electronic input. The second mode regards the
POCV operating in a proportional control regime, where the actuator moves to
any intermediate position corresponding to a varying electronic (analog)
signal. The POCV is smart as it operates using a shape memory alloy-based
actuator controlled by an electronic system that can be connected to cloud and
mobile platforms for remote data exchange, management, and diagnostics. The
shape memory alloy (SMA) wire is a soft, lightweight, and highly flexible
material that can recover from different temporary shapes to its initial shape
under temperature control variations [6]. Ultimately, the microfluidic device
aims to provide a precise throughput and handle high-pressure differences
compared to other microfluidic valves of similar size.

Though the POCV has proven effective in controlling the flow for the ON/OFF
mode, limitations during the proportional operating mode were noticed.
Experiments during the proportional opening process revealed the emergence of
severe vibrations as soon as the actuation reaches a specific position. These
vibrations are associated with noise and turbulent outflow and are undesirable
for the performance of the valve. Several trial-and-error procedures were
conducted to mitigate the problem, but without success. The inability to
address these issues is mainly attributed to the miniature size limitations,
and the lack of advanced equipment to precisely locate these vibrations.

Recent advances in numerical methods have shown promising prospects for
industrial decision-making on complex production technologies supported by
simulation-based engineering [7]. Specifically, computational simulations have
been employed to investigate valve dynamics, providing a deeper understanding
of mechanical behavior and giving insights into performance within industrial
processes. For instance, a $3$D check valve model is developed and tested
using numerical methods for fluid-structure interaction (FSI) in [8] to
analyze the sealing characteristics of the diaphragm volumetric pumps. A
numerical tool using computational fluid dynamics (CFD) is developed in [9] to
compute the mass flow capacities and opening characteristics of spring-loaded
safety valves at operating pressures up to $3600$ bar. The study is later
extended in [10] to cover FSI numerical methods to analyze the opening
characteristic and eventual instability issues to be avoided. On the other
hand, several studies focus on analyzing the valve flow response as a function
of membrane deflection at mesoscale and microscale levels. An FSI simulation
is performed in [11] to evaluate the impact of the membrane material
properties on the flow behavior in a microvalve. A microfluidic device with a
single micro-channel is fabricated and tested in [12] using FSI computational
methods to investigate the fluid interaction with the deformable wall. The
work of [13] exemplifies the significance of the structural material
properties and the topology of the inlet region in the microvalve performance.
Furthermore, a few studies focus on the pilot-operated valve performance: A
dynamic model is developed in [14] to analyze the bifurcation within a
hydraulic system composed of a pilot-operated chamber valve and a piston. The
analysis shows the conditions leading to instability and provides a guideline
to initiate an optimal design. A simulation tool is developed by [15] and
consists of a $1$D fluid model coupled with a structural piston to predict the
instability of a pilot-operated pressure relief valve. The results were
validated after a confrontation with experimental data. Nevertheless, it is
interesting to mention that only scant attention has been given to pilot-
operated valves at the microfluidic scale. Although a pilot-operated miniature
digital hydraulic valve with a high flow capacity and a fast response was
introduced in [16], authors are unaware of additional research studies
covering computational simulations of microfluidic pilot-operated control
valves, let alone instability issues at smaller scales. This can probably be
attributed to the computational challenges of modeling such problems using
multiphysics simulation tools, the limited knowledge about these devices’
working operations, or exclusive innovation copyrights in designing and
developing these products to operate within microfluidic circuits.

The novelty in the present contribution pertains to leveraging advanced
simulation capabilities to overcome experimental limitations and understand
the performance and dynamics of the microfluidic pilot-operated valve under
investigation. In particular, numerical FSI methods are employed in this work
to reproduce the interaction of the fluid and the structural membranes during
the valve ON–OFF and proportional regimes. The latter includes the onset of
flow-induced vibrations.

Numerical methods for FSI have gained significant interest over the past two
decades as an analysis instrument in various engineering and scientific
fields, e.g. aerospace and civil engineering, transportation, medical
applications, and the beverage industry. Several applications are concerned
with multiphysics problems that can be solved using numerical FSI methods,
such as the motion of an aircraft’s wings due to the actions of the air at
high-speed [17], the fluid’s effect on valves during operational conditions
[18, 11, 19], airbag deployment for passenger safety in vehicles [20], and the
investigation of cardiovascular disorders [21, 22]. Although modeling the
transient behavior of a moving structure interacting with a fluid remains a
challenging task, considerable progress has been achieved following the
advancements in computational fluid dynamics (CFD), computational solid
mechanics (CSM), and high-performance computing (HPC).

Numerical methods for FSI can generally be categorized into monolithic and
partitioned approaches. Monolithic methods have been propounded as the more
robust [23]. However, these methods are in fact abstruse in view of the need
to develop dedicated codes. In addition, the matrix system of equations is
generally dense and ill-conditioned [24] and, hence, monolithic methods can be
memory expensive [25]. Partitioned approaches, on the other hand, provide more
versatility in compromising between robustness and efficiency [26]. Moreover,
and in practice most relevant, partitioned methods preserve the modularity of
the fluid and solid subsystems, thus enabling the reuse of advanced CFM and
CFD software, e.g. commercially available software such as ANSYS (ANSYS Inc.,
Canonsburg, PA, USA) [27] or open-source software such as OpenFOAM. The
advantages of modularity do not just pertain to the software, but also to the
computational models that have been developed, often over periods of many
years, for computational structural and flow analyses [25].

To assess the suitability of partitioned FSI approaches based on commercially
available software for complex FSI problems such as the considered valve
system, in the present contribution we regard the use of the ANSYS simulation
software. The valve system under consideration carries two essential
complications for partitioned methods. First, the fluid is incompressible and
has a relatively large density relative to the membrane material. As a result,
the fluid carries a relatively large added mass which, owing to the
incompressibility, cannot be controlled by means of the time-step size [28,
29]. This means that the FSI problem is intrinsically strongly coupled.
Accordingly, a strongly coupled partitioned approach must be applied, in which
the fluid and solid subsystems are solved several times within each time step.
This is to be contrasted to weakly coupled FSI problems, which can be solved
utilizing a so-called staggered method, in which the fluid and solid are
solved only once per time step [30]. Strongly coupled partitioned methods are
generally based on a simple fixed-point iteration approach referred to as
subiteration [26]. Subiteration has proven successful for, for instance,
butterfly-valve-actuation simulations [31]. Within each time step, the FSI
solution obtained by subiteration in principle coincides with that of a
monolithic approach. However, the robustness and efficiency of subiteration
deteriorate as the inertial effects of the fluid on the structure become more
dominant, i.e., as the FSI coupling strength as expressed by the fluid-solid
mass ratio increases. In cases where the fluid-solid mass ratio is too large,
subiteration is unstable or converges prohibitively slowly. In such cases,
auxiliary techniques are required to enhance the subiteration convergence
behavior, e.g. artificial compressibility (AC) [32], Aitken’s $\Delta^{2}$
method [33], or the interface quasi-Newton inverse least squares (IQN-ILS)
methods [34, 35]; see [26] for an overview.

The second fundamental complication of the considered valve FSI problem in the
context of partitioned solvers is that in the closed configuration of the
valve, the fluid volume upstream of the flexible membrane is fully closed. The
standard partitioned approach for FSI relies on a so-called Dirichlet–Neumann
partition of the FSI interface conditions, in which the fluid tractions are
imposed on the solid as a Neumann (natural) boundary condition, and the solid
displacement is imposed on the fluid as a Dirichlet (essential) boundary
condition. If the valve is closed, the upstream fluid volume is subjected to
Dirichlet boundary conditions on its entire perimeter, which, in combination
with the incompressibility of the fluid, leads to a compatibility condition on
the solid deformation. This problem is in the FSI literature referred to as
the incompressibility dilemma [36]. It is to be noted that in the valve
problem, also nearly closed configurations occur, which carry similar
problems. The compatibility condition underlying the incompressibility dilemma
can be avoided by reverting to a Robin–Neumann (or Robin–Robin) partitioning,
in which the fluid is subjected to a mixed-type boundary condition at the
fluid-solid interface [37, 38, 39]. However, this partitioning strategy is
typically unavailable in black-box FSI solution procedures. Another approach
to circumvent the compatibility condition is to introduce artificial
compressibility of the fluid in the iterative procedure [40], or simply assume
compressibility of the fluid in the FSI model [41].

The numerical investigations in this paper focus on the performance of the
considered POCV valve for digital and analog signal inputs, i.e., in standard
ON/OFF operation and in a proportional-control setting. The computational
model takes into account the SMA-based actuator displacement acting on the
pilot valve and the main valve sealing properties over time. To account for
the sealing of the valve, in addition to the interaction of the fluid and
solid subsystems, the considered valve model also incorporates the contact
between the membrane and the valve seat and, hence, in fact corresponds to a
fluid-solid-contact interaction problem [42]. The fluid and structural models
(incl. contact) are simulated with the ANSYS Fluent and ANSYS Transient
Mechanical modules, respectively. The data mapping across each solver is
handled via ANSYS SystemCoupling. An advantage of the ANSYS software,
highlighting the benefits of reusing software modules that are enabled by the
partitioned FSI approach, is that it facilitates advanced geometrical modeling
of the valve system, provides functionality to handle contact, gives access to
nonlinear material models and means of calibrating such models, and provides
access to turbulence models. In addition, as the numerical model of the valve
system generally involves many degrees of freedom, one can exploit the various
parallel computing capabilities offered by ANSYS.

The remainder of the paper is organized as follows: Section 2 introduces the
microfluidic pilot-operated control valve and provides its characteristics and
specifications. The mathematical-physical model of the valve system and its
numerical approximation are elaborated in Section 3. This section also
considers the utilized coupling strategy. Section 4 presents the numerical
experiments and discusses the results. Finally, Section 5 presents concluding
remarks.

##  2 Pilot-operated control valve: Characteristics and operating conditions

The valve under investigation is a smart pilot-operated control valve (POCV)
developed by a high-tech start-up that designs and develops smart and
precision microfluidic systems. This smart valve is designed to provide
precise flow regulation and positioning accuracy during operation, with ON/OFF
and proportional control working regimes that can be effectively integrated
into hydraulic circuits for medical, food, and beverage applications. In
addition, the device is equipped with flow, pressure, or temperature sensors,
allowing closed-loop control of the customer’s desired function, such as
autonomous dosing, dispensing, profiling, or mixing, based on the signal from
the connected sensor. In the current state of the art, the control is
performed by the movement of a pilot membrane, activated by a shape memory
alloy-based actuator that changes the displacement of the moving part, which
enables an ON/OFF actuation with two positions (open or closed), or a
proportional actuation that can assume intermediate positions.

The valve is modular and consists of two medium-separated chambers, viz. the
main chamber and the pilot chamber. The pilot chamber part is connected to an
electronic circuit through a miniaturized shape memory alloy (SMA) material
wire. This smart material can change its shape under temperature variations
[4, 6]. The activated SMA material shape deforms the pilot membrane, resulting
in a pressure drop across the pilot stage which drives the opening of the main
valve. The pilot pressure variations and the spring stiffness contribute to
the actuation control in order to regulate the flow rate accurately and
deliver a predefined fluid throughput as required by the user.

Figure 1: Illustration of the actuation control based on shape memory alloy
material (Yellow). The left-hand side corresponds to a closed POCV; the pilot
membrane (Silver) is undeformed, and the upstream flow is blocked. The right-
hand side corresponds to an actuated valve; the pilot membrane deforms
according to the current SMA position, and the POCV switches to a partially or
fully open position.

To elucidate the SMA-based actuation control mechanism, Figure 1 shows the
positions of the pilot membrane with respect to the SMA actuator: If the
actuator is inactive, the POCV is in a normally closed mode (NC) [3].
Otherwise, the updated SMA shape drives the pilot membrane position to allow
fluid flow toward the outlet, which triggers the pressure change within the
valve. In addition, given the reversible SMA shape effect, the valve can be
opened and closed in a similar manner, and the two-way electro-thermal
mechanical actuation enables multitudinous cycles during its lifespan. It is
worth mentioning that the SMA-based actuation control yields remarkable
benefits for fluidic applications, e.g. noiseless actuation, compact design,
and lightweight products [6]. The schematics in Figure 2 depict the general
opening process of the valve in relation to the actuation control position;
Panel (b), in which the main valve opens partially, pertains to the
proportional actuation control for the POCV operation.

Figure 2: Schematic views of the pilot-operated valve under investigation for
different opening setpoints. The valve can be operated as ON/OFF valve ((a)
$\rightarrow$ (c)), and is designed to work as a proportional valve ((a)
$\rightarrow$ (b) $\rightarrow$ (c)). A: Inflow region; B: Outflow region; C:
Pilot chamber region; D: Pilot chamber outflow.

The POCV is designed to withstand high-pressure differences, ranging from
$0.2$ up to $5$ bar while ensuring efficient sealing properties within the
piping system. The technical information regarding dimensions and operating
parameters pertinent to the present investigation, in accordance with the
producer’s datasheet, are listed in Table 1.

Table 1: Technical information relevant to the pilot-operated valve under study. Geometrical data |
---|---
Inlet orifice diameter [mm] | $4$
Outlet orifice diameter [mm] | $8$
Main membrane diameter [mm] | $16.5$
Pilot membrane [mm] | $2.7$
Operating parameters |
Function | Normally closed
Activation | Pilot-operated
Control | ON/OFF
Pressure difference $\Delta P$ [bar] |  $0.2$ $\cdot\cdot$ $5$
Rise time [s] | $0.15$

##  3 Models and methods

In this section, we present the fluid-structure interaction (FSI) model of the
considered pilot-operated valve (POCV). In addition, we regard aspects of the
numerical treatment of the FSI model.

###  3.1 POCV model

The considered POCV configuration is illustrated in Figure 3. The POCV
comprises a fluid domain and a structural domain and is surrounded by rigid
exterior boundaries. The structural domain consists of two disjoint regions,
viz., the pilot membrane and the main membrane. In the open configuration of
the valve, the fluid domain is connected. In the closed configuration, the
fluid domain is disconnected, comprising two components: one connected to the
inlet of the valve, and one connected to the outlet. We denote by
$\Omega\subset\mathbb{R}^{3}$ the POCV domain, and by $\Omega_{\textsc{f}}(t)$
and $\Omega_{\textsc{s}}(t)$ the time-dependent domains occupied by the fluid
and solid, respectively. The solid domain comprises two disjoint parts,
$\Omega_{\textsc{s}}^{\textsc{m}}$ and $\Omega_{\textsc{s}}^{\textsc{p}}$,
corresponding to the main membrane and pilot membrane, respectively. The
fluid-solid interface $\Gamma(t)$ corresponds to the intersection between the
boundaries of the fluid and solid domains, according to
$\Gamma(t)=\partial\Omega_{\textsc{f}}(t)\cap\partial\Omega_{\textsc{s}}(t)$.
We denote by $\Gamma^{\textsc{m}}(t)$ (resp. $\Gamma^{\textsc{p}}(t)$) the
fluid-solid interface associated with the main (resp. pilot) membrane.

Figure 3: Illustration of the pilot-operated control valve.

###  3.2 Fluid subsystem

We consider a time interval $(0,T)$. Let
$\boldsymbol{u}:(0,T)\times\Omega_{\textsc{f}}(t)\to\mathbb{R}^{3}$ and
$p:(0,T)\times\Omega_{\textsc{f}}(t)\to\mathbb{R}$ denote the fluid velocity
and pressure, respectively. The fluid flow is described by the incompressible
Navier–Stokes equations:

$\displaystyle\rho_{\textsc{f}}\partial_{t}\boldsymbol{u}+\rho_{\textsc{f}}\boldsymbol{u}\cdot\nabla\boldsymbol{u}+\nabla
p-\nabla\cdot\boldsymbol{\tau}$ $\displaystyle=0$ $\displaystyle\text{in
}(0,T)\times\Omega_{\textsc{f}}(t)$ (1a)
$\displaystyle\nabla\cdot\boldsymbol{u}$ $\displaystyle=0$
$\displaystyle\text{in }(0,T)\times\Omega_{\textsc{f}}(t)$ (1b)

with $\rho_{\textsc{f}}$ as the fluid density. For closed or nearly closed
valve configurations, the fluid velocities are sufficiently small to consider
the flow as laminar. In this scenario, the stress tensor $\boldsymbol{\tau}$
coincides with the viscous stress and
$\boldsymbol{\tau}=\boldsymbol{\tau}_{\textsc{v}}:=\mu_{\textsc{d}}(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T})$
with $\mu_{\textsc{d}}$ the dynamic viscosity of the considered fluid. For
open configurations of the valve, the fluid tends to accelerate considerably
due to the high upstream pressure and the flow enters the turbulent regime,
viz., a regime in which the fluid exhibits flow separation and recirculation
regions. To account for turbulent effects in the open configuration, we invoke
the URANS model with Boussinesq closure. The stress tensor $\boldsymbol{\tau}$
in (1a) is accordingly decomposed into a viscous part and a turbulent part as
$\boldsymbol{\tau}=\boldsymbol{\tau}_{\textsc{v}}+\boldsymbol{\tau}_{\textsc{t}}$
where

$\boldsymbol{\tau}_{\textsc{t}}=\mu_{\textsc{t}}(\nabla\boldsymbol{u}+\big{(}\nabla\boldsymbol{u})^{T}\big{)}-\tfrac{2}{3}\rho_{\textsc{f}}k\boldsymbol{I}$
(2)

with $\mu_{\textsc{t}}$ the eddy viscosity, $k$ the turbulent kinetic energy
and $\boldsymbol{I}$ the identity tensor. The eddy viscosity is generally
closed in terms of the (averaged) fluid velocity by means of auxiliary
variables. In particular, we consider the $k$-$\omega$ shear stress transport
(SST) turbulence model, where $\omega$ represents the specific dissipation
rate [43]. In the $k$-$\omega$ SST model, the eddy viscosity is given by

$\mu_{\textsc{t}}=\frac{\rho_{\textsc{f}}ak}{{\rm max}(a\omega,\Omega F)}$ (3)

where $a$ and $\Omega$ are model constants that depend on the specific
implementation of the SST model, and $F:=F(k,\omega)$ is a specific function
to transition between the $k$-$\omega$ model and the $k$-$\varepsilon$ model
[44]. The SST turbulence model extends the fluid model by the aforementioned
two auxiliary variables $k$ and $\omega$, subject to:

$\displaystyle\rho_{\textsc{f}}\partial_{t}k+\rho_{\textsc{f}}\nabla\cdot{}k\boldsymbol{u}-\nabla\cdot(\Gamma_{k}\nabla{}k)$
$\displaystyle=G_{k}-Y_{k}$ $\displaystyle\text{in
}(0,T)\times\Omega_{\textsc{f}}(t)$ (4a)
$\displaystyle\rho_{\textsc{f}}\partial_{t}\omega+\rho_{\textsc{f}}\nabla\cdot{}\omega\boldsymbol{u}-\nabla\cdot(\Gamma_{\omega}\nabla\omega)$
$\displaystyle=G_{\omega}-Y_{\omega}$ $\displaystyle\text{in
}(0,T)\times\Omega_{\textsc{f}}(t)$ (4b)

where $\Gamma_{k}$ and $\Gamma_{\omega}$ represent the effective diffusivities
of $k$ and $\omega$, respectively. The nonnegative terms $G_{k}$ and
$G_{\omega}$ (resp. $Y_{k}$ and $Y_{\omega}$) represent certain production
(resp. dissipation) terms of turbulent kinetic energy and specific dissipation
energy, respectively. We refer to [44] for further details and interpretation
of the equations.

To accommodate the motion of the fluid domain in the fluid-structure
interaction problem, the equation of motion (1a) and the turbulent balance
laws (4) are reformulated in arbitrary Lagrangian-Eulerian (ALE) form; see,
for instance, [45, 46, 47, 48]. This part of the formulation is standard in
FSI and will not be further elaborated here.

The differential equations (1) and (4) must be furnished with suitable initial
and boundary conditions. The boundary of the fluid domain comprises the fluid-
solid interface at the main membrane, the wetted surface of the pilot
membrane, rigid walls, an inflow boundary and an outflow boundary; see Figure
3. The boundary conditions at the fluid-solid interface are considered
separately in Section 3.5. At rigid solid walls, we impose the usual no-slip
condition:

$\boldsymbol{u}=0\qquad\text{ at }(0,T)\times\Gamma_{\text{wall}}$ (5)

At the outlet, we impose a homogeneous traction condition according to:

$p\boldsymbol{n}-\boldsymbol{\tau{}n}=0\qquad\text{ at
}(0,T)\times\Gamma_{\text{out}}$ (6)

where $\boldsymbol{n}$ denotes the exterior unit normal vector on the boundary
of the fluid domain. It is to be noted that stability of the initial-boundary-
value problem requires that $\boldsymbol{u}\cdot\boldsymbol{n}\geq{}0$ during
the entire evolution of the flow field at parts of the boundary where (6) is
imposed. At the inlet, we impose pressure in accordance with the pressure jump
across the valve. Formally, this corresponds to a relation of the form (6)
with a non-zero right-hand side. However, as
$\boldsymbol{u}\cdot\boldsymbol{n}$ is generally negative at the inflow
boundary, the aforementioned stability condition associated with (6) is
violated at the inflow. Therefore, instead, at the inflow a scaled velocity
profile is imposed:

$\boldsymbol{u}(t,\boldsymbol{x})=U(t)\,\hat{\boldsymbol{u}}(\boldsymbol{x})\qquad\text{
at }(0,T)\times\Gamma_{\text{in}}$ (7)

where $\hat{\boldsymbol{u}}$ represents the prescribed inflow profile,
compatible with the no-slip condition at the intersections of
$\Gamma_{\text{in}}$ with the solid wall boundary $\Gamma_{\text{wall}}$. The
time-dependent scaling factor $U$ is such that the average normal component of
the traction at the inflow coincides with the prescribed inlet pressure
$p_{\text{in}}$:

$\frac{1}{\operatorname{mean}(\Gamma_{\text{in}})}\int_{\Gamma{\text{in}}}p-\boldsymbol{n}\cdot\boldsymbol{\tau{}n}=p_{\text{in}}$
(8)

At the section $\Gamma^{\textsc{p}}(t)$ of the boundary coinciding with the
wetted boundary of the pilot membrane (see Figure 3), the fluid domain is
deformed in accordance with the deformation of the pilot membrane, and the
fluid velocity coincides with the velocity of the pilot membrane:

$\boldsymbol{u}=\partial_{t}\boldsymbol{d}_{\textsc{s}}\qquad\text{at
}(0,T)\times\Gamma^{\textsc{p}}(t)$ (9)

where $\boldsymbol{d}_{\textsc{s}}$ stands for the deformation of the solid,
in particular, the pilot membrane. Boundary condition (9) can be conceived of
as so-called ‘one-way FSI’, in that the fluid domain deforms in accordance
with the deformation of the solid, and the fluid velocity coincides with the
solid velocity, but the solid does not respond to the traction exerted by the
fluid on the wetted boundary of the pilot membrane; see also Section 3.5.

A detailed treatment of boundary conditions for the SST turbulence model (4)
is beyond the scope of this work. At solid walls, including the pilot
membrane, homogeneous and non-homogeneous Dirichlet-type boundary conditions
are typically imposed on $k$ and $\omega$, respectively:

$\displaystyle k$ $\displaystyle=0\qquad$ $\displaystyle\text{at
}(0,T)\times(\Gamma_{\text{wall}}\cup\Gamma^{\textsc{p}})$ (10a)
$\displaystyle\omega$ $\displaystyle=\omega_{\text{wall}}\qquad$
$\displaystyle\text{at
}(0,T)\times(\Gamma_{\text{wall}}\cup\Gamma^{\textsc{p}})$ (10b)

where $\omega_{\text{wall}}$ represents an estimate of the dissipation rate at
the boundary. At the outflow boundary, homogeneous Neumann conditions are
imposed:

$\displaystyle\boldsymbol{n}\cdot\nabla{}k$ $\displaystyle=0\qquad$
$\displaystyle\text{at }(0,T)\times\Gamma_{\text{out}}$ (11a)
$\displaystyle\boldsymbol{n}\cdot\nabla{}\omega$ $\displaystyle=0\qquad$
$\displaystyle\text{at }(0,T)\times\Gamma_{\text{out}}$ (11b)

At the inflow part of the boundary, one generally imposes non-homogeneous
Dirichlet conditions on $k,\omega$, according to

$\displaystyle k$ $\displaystyle=k_{\text{in}}\qquad$ $\displaystyle\text{at
}(0,T)\times\Gamma_{\text{in}}$ (12a) $\displaystyle\omega$
$\displaystyle=\omega_{\text{in}}\qquad$ $\displaystyle\text{at
}(0,T)\times\Gamma_{\text{in}}$ (12b)

where $k_{\text{in}}$ and $\omega_{\text{in}}$ represent exogenous data,
corresponding to estimates of the turbulent kinetic energy and specific
dissipation rate at the inlet.

The initial conditions for (1) and (4) comprise a specification of
$\boldsymbol{u},k,\omega$ at $t=0$. In accordance with the fact that the valve
is initially in a closed configuration, we impose a homogeneous initial
condition on the velocity:

$\boldsymbol{u}(t=0)=0\qquad\text{in }\Omega_{\textsc{f}}(t=0)$ (13)

Because the fluid is stagnant in the initial closed configuration, the
turbulent kinetic energy and specific dissipation rate vanish initially:

$\displaystyle k(t=0)$ $\displaystyle=0\qquad$ $\displaystyle\text{in
}\Omega_{\textsc{f}}(t=0)$ (14a) $\displaystyle\omega(t=0)$
$\displaystyle=0\qquad$ $\displaystyle\text{in }\Omega_{\textsc{f}}(t=0)$
(14b)

###  3.3 Structure subsystem

The structure subsystem is illustrated in Figure 4. The structure comprises
three components, viz. the main membrane, the pilot membrane, and a quasi-
rigid contact dummy body. The dummy body represents the so-called valve seat,
which is the part of the casing of the POCV which comes into contact with the
main membrane. The dummy body is further elaborated in Section 3.6. The pilot
membrane is modelled as an independent actuator, in the sense that a
prescribed motion is imposed on part of its external (non-wetted) boundary,
representing the actuation, and the pilot membrane is ignorant to the presence
of the fluid at its wetted boundary. The main membrane, on the other hand,
does respond to the fluid-traction exerted on its wetted boundary. The main
membrane, therefore, represents the only deformable structure component in the
aggregated fluid-structure interaction problem. For the structural subsystem,
however, this essential difference between the main membrane and the pilot
membrane only manifest as a difference in the load. Hence, we present the
model for the structure subsystem, comprising both main and pilot membrane, in
a unified manner.

Figure 4: POCV structural system after spatial discretization. The main
membrane is constrained by a contact dummy body representing the valve seat
and a spring to guide the deformation. The pilot membrane is clamped and
controlled by a displacement boundary condition to replicate the SMA-based
actuator actions.

The membranes are composed of an incompressible rubber-like material. The
time-dependent deformation of the structural domain is rendered in an updated-
Lagrangian formulation in the numerical approximation. However, to facilitate
and condense the presentation, we present the model in the equivalent total-
Lagrangian formulation. The membranes exhibit dissipation under dynamic
deformation. The dissipation is modelled by Rayleigh damping. Because Rayleigh
damping is most conveniently expressed in a weak formulation, we regard the
structural model in weak form. Denoting by
$\hat{\Omega}_{\textsc{s}}=\Omega_{\textsc{s}}(t=0)$ the reference
configuration of the solid subsystem, the deformation
$\boldsymbol{d}_{\textsc{s}}:(0,T)\times\hat{\Omega}_{\textsc{s}}\to\Omega_{\textsc{s}}(t)$
satisfies:

$\displaystyle\int_{\hat{\Omega}_{\textsc{s}}}\rho_{\textsc{s}}\boldsymbol{w}\cdot\frac{\partial^{2}\boldsymbol{d}_{\textsc{s}}}{\partial
t^{2}}+\int_{\hat{\Omega}_{\textsc{s}}}\operatorname{Grad}\boldsymbol{w}:\boldsymbol{P}+\int_{\hat{\Omega}_{\textsc{s}}}\alpha\rho_{\textsc{s}}\boldsymbol{w}\cdot\frac{\partial\boldsymbol{d}_{\textsc{s}}}{\partial
t}\quad$
$\displaystyle+\int_{\hat{\Omega}_{\textsc{s}}}\beta\operatorname{Grad}\boldsymbol{w}:\boldsymbol{P}^{\prime}(\,\,\hat{\\!\\!\boldsymbol{d}}_{\textsc{s}},\partial_{t}\boldsymbol{d}_{\textsc{s}})$
$\displaystyle=l(\boldsymbol{w})\phantom{0}\quad\text{in}\
(0,T)\times\hat{\Omega}_{\textsc{s}}$ (15a)
$\displaystyle\int_{\hat{\Omega}_{\textsc{s}}}q(\det{}\boldsymbol{F}-1)$
$\displaystyle=0\phantom{l(\boldsymbol{w})}\quad\text{in}\
(0,T)\times\hat{\Omega}_{\textsc{s}}$ (15b)

for all admissible test functions $\boldsymbol{w},q$, where $\rho_{s}$ is the
solid mass density,
$\boldsymbol{F}=\operatorname{Grad}\boldsymbol{d}_{\textsc{s}}$ denotes the
deformation gradient, and $\boldsymbol{P}$ is the first Piola–Kirchhoff stress
tensor. The load functional in the right-hand side of (15a) vanishes for the
pilot membrane, and is further elaborated in Section 3.5 for the main
membrane. It is to be noted that the gradient operator in (15),
$\operatorname{Grad}$, acts in the reference configuration. The last two terms
in the left-hand side of (15a) represent Rayleigh damping, where
$\boldsymbol{P}^{\prime}(\,\,\hat{\\!\\!\boldsymbol{d}}_{\textsc{s}};\partial_{t}\boldsymbol{d})$
denotes the directional derivative of $\boldsymbol{P}(\boldsymbol{d})$ with
respect to $\boldsymbol{d}$, evaluated in the reference configuration
$\,\,\hat{\\!\\!\boldsymbol{d}}_{\textsc{s}}$ , in the direction
$\partial_{t}\boldsymbol{d}$. One may note that
$\,\,\hat{\\!\\!\boldsymbol{d}}_{\textsc{s}}(\boldsymbol{X})=\boldsymbol{X}$.
The coefficients $\alpha,\beta\geq{}0$ are referred to as the Rayleigh
coefficients. The damping term is such that in a finite-element approximation
of (15), it leads to a term $\mathbf{D}\,d_{t}\mathbf{d}$ with
$\mathbf{D}=\alpha\mathbf{M}+\beta\mathbf{T}_{0}$, where $\mathbf{M}$ and
$\mathbf{T}_{0}$ denote the mass matrix and tangent stiffness matrix at $t=0$,
respectively.

A suitable class of constitutive models to describe the elastic behavior of
the considered rubber-like material is provided by the hyperelastic material
models. Correspondingly, the first Piola–Kirchhoff stress tensor is the
derivative of an associated stored-energy density function
$W:=W(\boldsymbol{F})$ with respect to the deformation gradient. Denoting the
Green–Lagrange strain tensor by
$\boldsymbol{E}:=\boldsymbol{E}(\boldsymbol{F})=\frac{1}{2}(\boldsymbol{F}^{T}\boldsymbol{F}-\boldsymbol{I})$,
where $\boldsymbol{I}$ is the identity tensor, the material of the membranes
is characterized by a neo-Hookean constitutive relation, associated with the
stored-energy function:

$W(\boldsymbol{F})=C_{10}\operatorname{tr}(2\boldsymbol{E}_{\textsc{s}})-p_{\textsc{s}}\big{(}\sqrt{\det{}(2\boldsymbol{E}+\boldsymbol{I})}-1\big{)}$
(16)

where $C_{10}$ is a material constant and $p_{\textsc{s}}$ is pressure. It is
to be noted that $\operatorname{tr}(2\boldsymbol{E})=I_{1}-3$ with $I_{1}$ the
first strain invariant and
$\sqrt{\det(2\boldsymbol{E}+\boldsymbol{I})}=\det{}\boldsymbol{F}$. The
material parameter $C_{10}$, which characterizes the material response, is
calibrated on the basis of experimental data; see A.

In the POCV device, a spring is mounted to the main membrane to provide an
additional displacement-proportional load; see Figure 4. The spring is
modelled in the structure subsystem by a force
$-k(\boldsymbol{d}_{\textsc{s}}(t,\hat{\boldsymbol{X}}_{\textsc{spr}})-\hat{\boldsymbol{X}}_{\textsc{spr}})$,
where $k$ indicates the spring constant, and
$\hat{\boldsymbol{X}}_{\textsc{spr}}$ is the mount point of the spring in the
reference configuration. As elastic solids theoretically do not admit point
loads, the spring force is distributed on a region of the membrane surface,
$\Gamma_{\textsc{spr}}(t)\subset\partial\Omega_{\textsc{s}}(t)$. The spring is
hence represented by a traction:

$\boldsymbol{t}_{\textsc{spr}}(\boldsymbol{x})=\begin{cases}-k(\boldsymbol{d}_{\textsc{s}}(t,\hat{\boldsymbol{X}}_{\textsc{spr}})-\hat{\boldsymbol{X}}_{\textsc{spr}})/\operatorname{Area}(\Gamma_{\textsc{spr}})&\text{if
}\boldsymbol{x}\in\Gamma_{\textsc{spr}}\\\ 0&\text{otherwise}\end{cases}$ (17)

acting on the surface of the membrane.

The differential equations (15) must be provided with suitable initial and
boundary conditions. The boundary-conditions at the wetted boundary
$\Gamma^{\textsc{m}}(t)\subset\partial\Omega_{\textsc{s}}(t)$ of the main
membrane, corresponding to the fluid-solid interface, represented by the load
functional in the right member of (15a), are presented in Section 3.5. Let us
note that the aforementioned spring load is administered at a part of the
wetted boundary, i.e.,
$\Gamma_{\textsc{spr}}(t)\subseteq\Gamma^{\textsc{m}}(t)$. On parts of the
boundary where the membranes are fixed to the casing of the POCV,
$\Gamma_{\textsc{fix}}$, the structural configuration is identified with the
initial configuration, i.e.,

$\boldsymbol{d}_{\textsc{s}}(t,\boldsymbol{X})=\boldsymbol{X}\quad\text{for
all }(t,\boldsymbol{X})\in(0,T)\times\Gamma_{\textsc{fix}}$ (18)

For the main membrane, the fixed boundary and wetted boundary are
complementary. For the pilot membrane, the remaining part of the boundary,
including the wetted boundary, is subject to homogeneous traction (i.e.,
stress-free) boundary conditions, except for a small region at the center of
the external (non-wetted) boundary, where the deformation in accordance with
the actuation is imposed.

Suitable initial conditions for (15) are provided by a specification of the
initial deformation, $\boldsymbol{d}_{\textsc{s}}(t=0)$, and the initial
velocity, $\partial_{t}\boldsymbol{d}_{\textsc{s}}(t=0)$. Noting that the
configuration of the membrane at $t=0$ coincides with the initial
configuration, it holds that:

$\boldsymbol{d}_{\textsc{s}}(t=0,\boldsymbol{X})=\boldsymbol{X}\quad\text{for
all }\boldsymbol{X}\in\Omega_{\textsc{s}}(t=0)$ (19)

We assume that both membranes are initially at rest and, hence,

$\partial_{t}\boldsymbol{d}_{\textsc{s}}(t=0)=0\quad\text{in
}\Omega_{\textsc{s}}(t=0)$ (20)

###  3.4 Contact formulation

The deformation of the main membrane is confined by the rigid parts of the
POCV, notably the so-called valve seat. In the closed configuration of the
POCV, the main membrane is in contact with the valve seat in such a manner
that the fluid domain is separated into two disjoint parts. In the fully open
configuration, there is no contact between the membrane and the valve seat.
The deformation of the elastic main membrane is subject to a non-penetration
condition of the form

$\operatorname{dist}(\partial\Omega_{\textsc{s}}(t),\Gamma_{\textsc{vs}})\geq{}0,$
(21)

where $\Gamma_{\textsc{vs}}\subset\partial\Omega_{\textsc{f}}(t=0)$ denotes
the part of the rigid boundary of the POCV corresponding to the valve seat. In
the fully open configuration, the inequality in (21) is strict, i.e., the
distance between the membrane boundary and the valve seat is strictly
positive. Contact occurs if equality holds in (21). An extensive review of
contact mechanics is beyond the scope of the present work: the reader is
referred to [49] and [50] for a review of existing contact formulations.

A fundamental complication in numerical procedures for fluid-solid-contact
interaction problems [42], pertains to the fact that the topological change in
the fluid domain occurring at contact (from a simply connected set to
connected set with holes [partial contact] or a disconnected set [fully
closed]) manifests in the computational method as a collapse of the mesh in
the fluid domain. To avoid such a collapse of the mesh, it is generally
necessary to ensure a small-but-finite separation of the surfaces, i.e., that
the inequality in (21) holds strictly. To this purpose, we modify the contact
formulation in a way that it retains (aims to retain) a thin gap of thickness
$\epsilon>0$ between the valve seat and the main membrane; see Figures 4 and
5. This implies that the contact load is activated if the distance between the
main membrane and valve seat is less than or equal to $\epsilon$ thus ensuring
a small-but-finite separation between the main membrane and the actual valve
seat at contact. Simultaneously, the fluid domain in the contact region is
replaced by a Darcy-type flow resistance if contact is detected, to mimic the
sharp increase in the flow resistance that occurs as the separation between
the valve and the valve seat approaches zero.

Contact is modelled by means of an auxiliary contact traction on the elastic
structure. The valve seat is represented by a quasi-rigid dummy body; see
Figures 4 and 5. Noting that the fluid between the two contacting surfaces
effectively lubricates the contact, we ignore friction and, hence, the contact
force acts in the direction perpendicular to the contact surface. The
considered contact traction (force per unit area) is of the form:

$\boldsymbol{t}_{\textsc{c}}(\boldsymbol{x})=-\boldsymbol{n}_{\textsc{s}}\big{(}k_{\textsc{c}}\big{\lfloor}(\boldsymbol{x}-\mathcal{P}\boldsymbol{x})\cdot\boldsymbol{n}_{\textsc{s}}+\epsilon\big{\rfloor}+\lambda_{\textsc{c}}\big{)}$
(22)

where $\mathcal{P}\boldsymbol{x}$ denotes the projection of the point
$\boldsymbol{x}$ on the surface of the membrane
$\partial\Omega_{\textsc{s}}(t)$ (resp. the surface of the dummy body
$\Gamma_{\textsc{d}}(t)$) along the normal vector onto
$\Gamma_{\textsc{d}}(t)$ (resp. $\partial\Omega_{\textsc{s}}(t)$), the
operator $\lfloor\cdot\rfloor=\max(\cdot,0)$ is the positive part (or Macaulay
bracket) and $\lambda_{\textsc{c}}$ is the contact force, i.e., the Lagrange
multiplier associated with the non-penetration condition. The first term in
parentheses in (22) represents a penalty contribution, the penalty parameter
$k_{\textsc{c}}>0$ being analogous to an elastic-spring constant. By virtue of
the Macaulay bracket and the offset parameter $\epsilon$, the contact traction
is active only in regions where the signed distance between the membrane and
the dummy body is less than $\epsilon$.

Figure 5: Membrane to valve seat inter-body contact modelling. Blue: Contact
dummy body section. Grey: Main membrane section. The mechanical parts are
separated by a thin gap of thickness $\epsilon>0$ equal to the fluid gap width
for consistency.

###  3.5 Interface conditions of the FSI problem

The fluid and the structure are interconnected by kinematic and dynamic
interface conditions [45]. The kinematic conditions prescribe that the
boundaries of the fluid and solid domains coincide at the fluid-solid
interface, and that the fluid velocity coincides with the structure velocity
at the interface. The latter condition implies,

$\boldsymbol{u}=\partial_{t}\boldsymbol{d}_{\textsc{s}}\quad\text{in
}(0,T)\times\Gamma(t)$ (23)

Let us note that on the left-hand (resp. right-hand) side of (23), positions
are to be understood in the actual (resp. reference) configuration. The
kinematic condition (23) can be conceived of as an essential (Dirichlet-type)
boundary condition for the fluid subsystem.

The dynamic condition enforces the balance of tractions at the fluid-solid
interface. For the main membrane, the load functional in (15a) therefore
comprises the fluid traction, in addition to the traction exerted by the
spring and by contact, according to (17) and (22), respectively:

$l(\boldsymbol{w})=\int_{\Gamma(t)}\boldsymbol{w}\cdot(p\boldsymbol{n}_{\textsc{f}}-\boldsymbol{\tau}\boldsymbol{n}_{\textsc{f}})+\int_{\Gamma(t)}\boldsymbol{w}\cdot\boldsymbol{t}_{\textsc{c}}+\int_{\hat{\Gamma}}\boldsymbol{w}\cdot\boldsymbol{t}_{\textsc{spr}}$
(24)

It is to be noted that the first and second terms in (24), corresponding to
the traction exerted by the fluid on the solid and the contact traction, are
formulated in the actual configuration, while the third term, associated with
the distributed spring load, is expressed in the reference configuration. The
dynamic condition in (24) can be regarded as a natural (Neumann-type) boundary
condition for the structure subsystem.

###  3.6 Numerical approximation of the fluid and solid subsystems

The fluid calculations are executed using the commercial software package
ANSYS Fluent and are based on an approximation of the URANS transport
equations by means of a finite-volume method in an ALE formulation. A second-
order upwind scheme is used for the spatial discretization of the convective
terms and a least-squares cell-based formulation is used to approximate the
diffusive fluxes in the momentum equation.

Figure 6 illustrates the fluid-mesh partition from a plane-section-view
perspective. The mesh combines different element types. The bulk is composed
of octree hexagonal elements. Boundary-layer prisms are used in the vicinity
of rigid boundaries. The regions adjacent to the moving surfaces are mostly of
general polyhedral type. The region adjacent to the fluid-solid interface is
refined to capture the large gradients in this area, as shown in panels (b)
and (c) of Figure 6. This heterogeneous mesh layout serves to enhance the
accuracy while restricting the number of degrees of freedom in the discrete
approximation. The fluid mesh used in the numerical simulations comprises
$1.019.544$ cells.

Figure 6: Fluid mesh from a plane section view. The stationary (non-moving)
zones are refined using inflation layers to improve accuracy. The deforming
regions are refined using smaller cells to reduce the probability of mesh
folding due to excessive motion. (b) Detail 1 shows the fluid-structure
interfaces and the adjacent deforming zones. (c) Detail 2 depicts the pilot
membrane and adjacent moving domain. Full closure of the valve is performed
using a wall boundary condition (pilot region) and a gap model (main flow
section).

The pilot channel is closed using a wall boundary zone, while the main
membrane blocks the flow using a gap model. This gap model is active when the
distance between the FSI wet interface and the valve seat is below a
predefined contact distance threshold; see also section 3.4. A spring-based
smoothing method controls the deforming mesh dynamics. To improve the
robustness of the mesh-deformation algorithm, a strict tolerance is used on
the dynamic mesh convergence criterion, and under-relaxation of the mesh
motion parameter is applied.

The solution procedure for the discretized fluid equations is based on
concurrent pressure-velocity coupling, instead of segregated approaches such
as SIMPLE and PISO. The choice for a coupled solution procedure is motivated
by its improved robustness relative to segregated solution approaches, which
is particularly important in view of the moving mesh and the strongly varying
flow velocity (encoded by the local Reynolds number) within the flow field,
especially during the valve-closing event. In our experience, concurrent
pressure-velocity coupling yields a better compromise between computational
complexity, time-step-size restrictions, and robustness, compared to
segregated solvers. To retain accuracy and robustness throughout the evolution
of the valve configuration, it is important to ensure convergence of the CFD
solver to a strict tolerance. In particular, in addition to the ANSYS Fluent
default criterion that the scaled residual drops below $10^{-3}$, we require
that the relative pressure updates near the fluid-solid interface and near the
pilot membrane drop below $10^{-4}$.

The structure is discretized using tetrahedral Taylor-Hood elements with
quadratic displacements and linear pressures. The mesh partition of the main
membrane is displayed in Figure 4. The main-membrane domain is locally refined
near the valve-seat contact region and contains $24.551$ elements. The pilot-
membrane’s mesh is uniformly distributed and comprises $1.555$ elements.

Because of differences in their geometric representation, ANSYS Fluent does
not admit to directly monitor contact between an endogenous structural element
and an exogenous rigid surface, e.g. the main membrane and the casing of the
POCV. To facilitate the contact formulation in ANSYS Fluent, the valve seat,
which is the part of the casing of the POCV that comes into contact with the
main membrane, is therefore represented as a separate quasi-rigid dummy body
in the computational model.

The solution procedure for the structure subsystem, including contact,
consists of a standard Newton procedure, augmented with a time-step-reduction
procedure which acts as a fail-safe in case of divergence or excessively slow
convergence of the Newton process. If time-step reduction is invoked, the
structure subproblem is subcycled until it is synchronized with the fluid,
i.e., multiple smaller time steps are carried out. The linear tangent problems
in the Newton procedure are solved by means of a sparse direct solver.

###  3.7 FSI solution strategy

Approaches for solving fluid-structure interaction problems can generally be
classified as monolithic or partitioned [47, 45]. In the monolithic approach,
the aggregated system of equations is solved in a fully coupled manner at each
time step, i.e., the fluid and structure subsystems are solved simultaneously
by means of one solver. The main advantage of monolithic methods is their
robustness. However, a significant disadvantage is the inherent loss of
modularity, obstructing the reuse of established solution strategies and
software frameworks. In addition, monolithic methods generally carry a high
computational complexity [45]. Partitioned methods are iterative procedures,
in which the fluid and structural subsystems are solved separately, and which
(in principle) yield the solution of the coupled FSI problem upon convergence.
The standard iterative procedure is referred to as subiteration, and consists
in solving the fluid subsystem subject to the kinematic condition, and the
structure subsystem subject to the dynamic condition. By virtue of the
intrinsic partitioning, partitioned methods are inherently modular, and enable
reuse of existing methodologies that have been developed for the fluid and
solid subsystems. Consequently, partitioned methods are much more flexible and
versatile than monolithic methods. A further subcategorization of partitioned
methods can be considered, into loosely coupled and strongly coupled
techniques. In loosely coupled (or staggered) approaches [47, 51], the fluid
and solid subsystems are solved only once per time step. In strongly coupled
approaches [26], the iteration is continued until a prescribed convergence
criterion is satisfied.

The coupling strength in FSI problems is generally characterized by the so-
called added mass. This characterization is based on the notion that the force
exerted by the fluid on the solid is essentially proportional to the
acceleration of the solid, and the fluid therefore acts on the structure as an
added mass [28, 29, 52]. In partitioned methods, the fluid and the solid are
treated asynchronously and, hence, the added mass of the fluid is treated
explicitly. Consequently, strongly-coupled FSI problems, i.e., FSI problems
with large added-mass effects, form a vulnerability of partitioned methods,
leading to divergence or prohibitively slow convergence. Various techniques
have however been developed to stabilize and accelerate strongly-coupled
partitioned methods based on subiteration, e.g, underrelaxation, Aitken’s
method, Quasi-Newton methods [34, 35], or (interface) artificial
compressibility [32, 53] ; see [26] for an overview. It is to be noted that
artificial-compressibility methods pertain specifically to incompressible FSI,
based on the fact that the incompressibility of the fluid causes its added-
mass to be essentially independent of the time step [29].

An additional problem pertaining to enclosed incompressible FSI problems, is
the so-called incompressibility dilemma [36]. In enclosed incompressible FSI
problems, the structure can only assume deformations that are compatible with
the volume of the fluid. This compatibility condition [54] causes severe
instability in partitioned methods, if it is not properly handled. This
instability is essentially independent of the added mass of the fluid, and it
cannot be resolved by most of the stabilization techniques intended to address
added-mass effects, such as under-relaxation, quasi-Newton methods, and
Aitken’s method. Artificial-compressibility stabilization is effective, as it
replaces the incompressible FSI problems by a sequence of compressible ones.

The considered POCV problem represents a strongly coupled FSI problem which,
in addition, displays characteristics of the incompressibility dilemma. The
fluid carried by the POCV is typically aqueous, and is essentially
incompressible. The mass density of the material of the main membrane is
similar to that of the carried fluid. Hence, the POCV system displays a
significant added-mass effect. Moreover, in the closed configuration of the
POCV, the upstream fluid volume is essentially closed. In view of the
aforementioned properties of the POCV system, we opt for the so-called
solution stabilization technique provided by ANSYS Fluent, which is similar to
the artificial-compressibility method [12].

##  4 Numerical results

In this section, we analyze the POCV on the basis of numerical investigations.
Experimental tests have revealed that the POCV is dysfunctional in
proportional-operation mode, in that it exhibits severe vibrations as soon as
the flow rate reaches a specific setpoint, depending on the differential
pressure. Experimental investigation of the origins of these vibrations is
impracticable, due to the miniaturized valve size, and the perturbing effects
of measurement devices. Numerical simulations therefore provide the only
viable means of investigating the behavior of the POCV in proportional-
operation mode, to determine the source of the vibrations.

Before applying the numerical model to investigate the behavior of the POCV in
proportional-operation mode, we first regard it in the simpler settings of a
fully closed configuration, and in a binary (ON/OFF) scenario. The equilibrium
solution of the closed configuration also serves as initial condition for the
binary and proportional-control scenarios.

###  4.1 Normally closed valve

We consider the POCV in the closed configuration, for inlet pressures
$p_{\text{in}}\in\\{20,50,100,200,300\\}$ (kPa); see Equation (8) and Table 2.
The closed configuration at the various inlet-pressure levels is determined by
a sequence of events. The initial configuration of the POCV corresponds to an
undeformed main membrane, while the pressure difference across the channel is
set to zero. In addition, the dynamic condition (24) is disabled, so that the
main membrane only deforms due to the compressed spring and the contact
forces. While the main membrane deforms, the fluid domain deforms accordingly,
i.e. the main membrane and the fluid are one-way coupled. After approximately
$5\text{\,ms}$, the main membrane has closed and reached equilibrium. Next,
the inlet pressure is set to its prescribed value, and the dynamic condition
is enabled. A two-way FSI follows, which deforms the main membrane into its
equilibrium closed configuration corresponding to the imposed inlet-pressure
level. The end time is set to $T=100\,\text{ms}$, at which time the main
membrane has essentially reached equilibrium.

The time-step size that is used in the simulation depends on the dynamics of
the structure and the fluid which, in turn, depend on the imposed inlet
pressure. The time-step size is based on the Courant number of the fluid, the
mesh velocity, and the fluid and contact force magnitudes. Table 2 summarizes
the time-step size for each inlet pressure.

Table 2: Inlet pressures and corresponding time-step sizes for the FSI simulations. Pressure inlet (kPa) | $20$ | $50$ | $100$ | $200$ | $300$
---|---|---|---|---|---
Time-step size ($\mu$s) | $500$ | 100 | 100 | 50 | 10

Figure 7 presents the fluid pressure in the closed configuration of the POCV
for $p_{\text{in}}$ equal to $0,50,100$ and $200\,$kPa, in equilibrium. The
fluid and solid domains are displayed in the actual deformed configuration.
The plots show that the pressure is essentially piecewise constant in two
disjoint regions: upstream of the main membrane, the pressure is equal to the
inlet pressure, whereas downstream the pressure is equal to the outlet
pressure. This implies that the gap model adequately sustains the pressure
jump across the contact region. From Figure 7, one can moreover observe that
the deformation of the main membrane increases as the pressure difference
between the inlet and outlet increases.

Figure 7: Plane section view of the normally closed POCV at 0 bar (top left),
0.5 bar (top right), 1 bar (bottom left) and 2 bar (bottom right). Colors
indicate pressure distribution.

Figure 8 plots the deformation of the top surface of the main membrane, in a
cross-sectional view. One can observe that the largest displacements occur at
the center of the membrane and that the deformation increases with the inlet
pressure and, accordingly, the pressure difference across the main membrane.
Concomitant with these high-pressure differences are large contact forces to
prevent interpenetration of the membrane and valve seat or, more specifically,
the dummy body. Figure 8 conveys that at an inlet pressure of $300$ kPa the
implemented contact procedure is not able to completely prevent the
penetration of the membrane into the dummy body in the vicinity of the valve
seat, which is located at $\pm{}5\pm[0,0.5]$ mm. In this region, the maximum
vertical displacement at $300$ kPa is approximately $0.05$ mm. However, by
virtue of the dummy body, which has thickness $0.5$ mm, a thin but finite gap
remains between the main membrane and the valve seat, which is adequate to
avoid degeneration of the fluid mesh in the gap.

Figure 8: Membrane deflection of the normally closed valve for inlet pressure
values $p_{\text{in}}\in\\{0.5,1,2,3\\}$ bar.

###  4.2 ON/OFF operation

The POCV is operated in a binary ON/OFF mode by means of an SMA-based
actuator, which drives the pilot valve into either a closed or an open
position without intermediate setpoints. The POCV working principle is then
similar to that of a solenoid actuated control valve [55]. In this section, we
consider numerical simulations of the opening process of the POCV. The
simulation results are compared to the corresponding experimental results to
assess the accuracy of the FSI model. In practice, the opening process of the
POCV proceeds according to the following stages:

  1. 0.

The POCV is in normally closed mode;

  2. 1.

An input signal triggers the SMA actuator;

  3. 2.

The actuator shifts the pilot membrane into the open setpoint;

  4. 3.

The pressure in the pilot channel decreases due to the increasing flow speed;

  5. 4.

The main membrane deforms in the valve opening direction due to the changes in
the pilot pressure;

  6. 5.

The main membrane position and the pressure difference determine the flow rate
through the valve.

To parallel the aforementioned stages in the numerical simulation of the
valve-opening process, we proceed as follows: We start the simulation from a
normally closed valve state as described in Section 4.1. Next, the pilot
membrane and the adjacent fluid mesh are deformed using one-way FSI; see also
Section 3.3. The resulting geometry represents the open configuration of the
pilot membrane. However, at this stage, the pilot bypass is still closed by
means of an auxiliary wall boundary condition. In the next step, this
auxiliary wall boundary in the pilot bypass is suppressed. The fluid can now
flow through the pilot chamber. As a result of this flow, the pressure at the
upstream wetted boundary of the main valve decreases. By virtue of the dynamic
condition (24), the corresponding load on the main membrane deforms the
membrane in the valve-opening direction. The fluid domain deforms accordingly,
by virtue of the kinematic condition (23). The contact between the valve seat
and the main valve is released, and the flow can now pass directly from the
upstream region to the downstream region without traversing the pilot chamber.
The fluid flow and the main membrane interact in a dynamic fluid-structure-
interaction process until a steady state is reached.

To demonstrate the valve’s response to the opening event, Figure 9 displays
four sequential snapshots of the pressure distribution across the fluid-
structure interfaces and the pilot channel at various time instants. These
snapshots are observed for a pressure inlet of $p_{\text{in}}=3$ bar and
display, respectively, the pressure distribution in normally-closed (NC)
configuration; immediately after opening the pilot valve at $t=10\,m$s; at
$t=20\,m$s after opening the pilot valve; and after equilibration in the open
position, at $t=200\,m$s. Comparing the top-left panel and the top-right
panel, one can observe that as the actuator is switched to the open position
and the pilot membrane opens slightly, the pressure in the pilot tube that
connects the pilot chamber to the outlet (located on the right side of the
figures) immediately drops. Comparing the top-right panel to the bottom-left
panel, one can observe that due to the opening of the pilot membrane, the
pressure in the pilot chamber and, correspondingly, on the bottom side of the
main membrane, rapidly drops from $3$ bar to approximately $1.7$ bar, in less
than $20$ ms. The pressure drop induces a change in the load exerted on the
main membrane, resulting in the main membrane undergoing deformation and
transitioning into a fully open configuration (bottom-right panel). The open
configuration represents a new equilibrium state of the valve system which
depends on the applied inlet pressure.

((a)) Fully closed valve

((b)) Pressure at 0.10005s

((c)) Pressure at 0.12s

((d)) Pressure at 0.19995s

Figure 9:  Pressure contour snapshots of the POCV. The top right panel shows
the pressure drop at the vicinity of the pilot membrane shortly after opening
actuation, the bottom left and right panels display the pressure at different
times corresponding to a fully open valve. Note the valve’s fast time-to-open
response from $3$ bar to $1.7$ bar (around $20$ ms).

Figure 10 shows the streamlines of the flow in the POCV at different time
instances, for an inlet pressure of $3$ bar. At $t=75$ ms, the POCV is still
closed and the fluid velocity essentially vanishes throughout the entire POCV.
At $t=105$ ms, shortly after the opening of the pilot membrane, the fluid
exhibits significant flow velocities, primarily near the valve seat and along
the tube that connects the pilot valve to the main valve. The right-most panel
of Figure 10 corresponds to the steady state of the flow in a fully open
configuration. Noticing the similarity of the streamline patterns in the
center and right-most panels, one can infer that after opening, the flow very
quickly settles into a stable, steady regime. In the open configuration, most
streamlines connect the outlet to the inlet via the gap between the main
membrane and the valve seat, and only a few pass through the pilot stage. This
indicates that in the open configuration, the secondary flow via the pilot
stage is limited. Noting that in the open configuration, the velocity at the
inflow and outflow regions is essentially uniform, the simulation results
predict a smooth and stable flow at the outlet.

Figure 10: Streamlines during the opening of the microfluidic POCV at $t=75$
ms (left), $t=105$ ms (center) and $t=175$ ms (right) for an inlet pressure of
$3$ bar.

To validate the results of the numerical simulations in relation to the
response of the actual POCV device, Figure 11 presents a comparison between
the computed and experimentally obtained pressure differences across the
valve, for a differential pressure of $3$ bar. In view of the uncertainty of
the exact location of the pressure measurement in the experimental setup, e.g.
due to the size of the probe, in the simulation the pressure has been probed
at three representative distinct, nearby locations. These locations are
indicated in the inset in Figure 11. The computed pressure at the inlet indeed
coincides with the imposed inlet pressure of $3$ bar. In the closed
configuration, the pressure is not uniform in the inlet region, but instead,
it decreases as the gap between the main valve and the valve seat is
approached. This suggests that there is some leakage through the gap in the
numerical model. Bearing in mind the small dimensions near the gap, minute
flow velocities can lead to significant variations in the pressure. It is
notable that the experimentally obtained pressure is also non-uniform: in the
closed configuration, the pressure at the probe drops to $298$ kPa, i.e. $2$
kPa less than the imposed inlet pressure of $300$ kPa. This suggests that the
experimental setup also exhibits some minute leakage. The computed results
indicate that the steepness of the pressure drop during the opening of the
main valve and the post-opening pressure level depends sensitively on the
location. In addition, the results indicate that the pre-to-post opening
pressure drop is non-uniform with respect to the distance to the valve-seat
gap. The computed post-opening pressure level at the measurement point closest
to the gap agrees well with the experimental data. However, the pre-opening
pressure level at this location is significantly lower than in the
experiments. The strong sensitivity of the pressure readings to the location
and to the leakage through the valve-seat gap, evident from the simulation
results, in combination with the uncertainty in the location of the
measurement device, indicates that the pressure-drop comparison in Figure 11
can only serve a qualitative purpose.

Figure 11: Pressure evolution during opening of the POCV at differential
pressure $3$ bar: Computed pressure at 3 points near the valve-seat gap and at
the inlet, and experimental data.

To further assess the agreement between the numerical model and the actual
POCV device, Figure 12 presents the temporal variation of the flow rate
through the POCV in ON/OFF operation for differential pressures of $2$ bar and
$3$ bar. For the inlet pressure of $2$ bar (resp. $3$ bar), the sudden
increase in the flow rate from an initial value of $0$ l/h to approximately
$150$ l/h (resp. $200$ l/h), which occurs immediately after the control pilot
valve is opened, inducing the opening of the main valve, aligns well with the
experimental results. For inlet pressure $3$ bar, the computed result closely
follows the more gradual increase in the flow rate in the time interval
$100-120\,m$s, exhibited by the experimental data. For the experimental
results, the increase of the flow rate however extends to a slightly longer
time interval than for the computed results. Consequently, the experimentally
observed equilibrium flow rate exceeds the computed flow rate by approximately
10%. Conversely, for the $2$ bar case, the computed equilibrium flow rate
exceeds the experimental flow rate by approximately 10%. In view of the
sensitivity of the flow rate to the details of the valve deformation near the
valve seat, we assess the overall agreement between the computed and
experimental flow rate evolutions as good.

Figure 12: Computed flow rate and experimentally obtained flow rate for a
differential pressure of $2$ bar (left) and $3$ bar (right).

###  4.3 Proportional-control operation

The considered POCV in principle allows to be operated as a proportionally
controlled valve, providing the option to dynamically adjust the flow as
required by the user by means of an analog input signal instead of the binary
signal that is used to operate the valve in ON/OFF mode. Noting that the POCV
is controlled by a voltage applied to the SMA-based actuator, the analog
signal can assume intermediate values in the voltage range
$[0,V_{\mathrm{max}}]$, where $V_{\mathrm{max}}$ corresponds to the voltage at
which the control valve is fully open. Proportional operation greatly extends
the operational range of the microfluidic POCV to, for instance, dosing and
mixing tasks that require precise flow modulation. However, the POCV under
investigation is deficient in proportional-control operation, in that the
device exhibits severe vibrations once the differential pressure exceeds a
certain threshold. In this section, we apply the developed FSI model to
examine the vibrations in the POCV in proportional-control mode.

In Section 4.2 it was shown that after opening, the flow in the POCV very
quickly equilibrates. As the envisaged proportional-control operation is
relatively slow compared to this equilibration time, the proportional-control
operation can essentially be regarded as a sequence of quasi-stationary
operations. To facilitate the analysis, we therefore consider a step-wise
progressive opening (negative displacements) of the control valve, as depicted
in Figure 13. In view of the minimal gap that must be retained to accommodate
the fluid mesh, an initial gap of approximately $20\,\mu$m is applied,
corresponding to $5$% of the displacement of the pilot membrane in the fully
open position. However, the initial gap is closed by means of an auxiliary
wall boundary condition; see also Section 4.2. At $t=0.08$ s, this auxiliary
condition is suppressed, and the pilot valve essentially opens to $20\,\mu$m,
without displacement of the pilot valve. The opening is subsequently step-wise
increased in accordance with Figure 13, by imposing a displacement on the
pilot valve at the mounting area of the SMA-based actuator. Let us note that
the maximum negative displacement at $-0.4\,m$m corresponds to the opening of
the pilot valve considered in the ON/OFF scenario in Section 4.2.

Figure 13: Step-wise displacement of the pilot membrane in the proportional-
control scenario of the POCV. The opening process follows three steps: ($1$):
Establishing the closed state of the valve system until $t=0.08$ s (red
triangle); ($2$): Suppressing the wall boundary condition that models the
closed pilot spool and establish setpoint = $20\,\mu$m for $t\in(0.08,0.12)$
s; ($3$): opening of the control valve following the imposed step-wise
displacement of the SMA-actuator mount point for function $t\in(0.12,0.125)$
s. The valve displacement of $-0.4\,m$m corresponds to the fully open setpoint
in the ON/OFF scenario.

To establish that a minimal pressure differential is required for the POCV to
function, we first regard an inlet pressure of $0.5$ bar. Figure 14 displays
the evolution of the flow rate across the outlet and the pressure at a
specific point within the pilot chamber during the proportional opening
scenario in accordance with the step-wise displacement of the pilot membrane;
see Figure 13. One can observe that the opening of the pilot membrane at
$t=0.08$ s, via removal of the auxiliary wall boundary condition, induces a
pressure drop of approximately $5\,k$Pa. In addition, the flow rate
instantaneously increases to $1.5\,$l/h. The first step in the displacement of
the pilot membrane at $0.12$ s again yields a sharp decrease in the pressure
in the pilot chamber, to approximately $30\,k$Pa. Simultaneously, the flow
rate increases gradually to approximately $8\,$l/h. The subsequent step-wise
increments of the displacement of the pilot membrane, bear only a moderate
effect on the pressure in the pilot chamber, and a nonuniform effect on the
flow rate. At the maximum displacement of the pilot membrane at $t=0.25$ s,
the pressure in the pilot chamber (at the measurement point) is approximately
$27\,k$Pa, while the flow rate levels off at approximately $8.5$ l/h. It is
notable that this flow rate is disproportionally lower than the observed flow
rates for $2$ and $3$ bar in Section 4.2. To further elucidate the operation
of the POCV at $0.5\,$bar, we regard in Figure 15 snapshots of the velocity
streamlines before opening, at $t=0.075$ s, immediately after opening, at
$t=0.09$ s, and at full opening, at $t=0.25$ s. One can observe that after
opening the pilot valve, the streamlines pass through the pilot bypass only.
This indicates that the pressure differential of $0.5$ bar is insufficient to
unseal the main valve. This also explains why the flow rate at $0.5$ bar is
disproportionally less than the flow rate for $2\,$bar and $3\,$bar observed
in Section 4.2, and the unsteadiness of the flow in the intermediate and fully
open configuration.

Figure 14: Flow-rate and pressure evolution during the opening of the POCV in
accordance with Figure 13 at differential pressure $0.5$bar. Figure 15:
Streamlines during opening of the POCV according to Figure 13 at $t=75$ ms
(left), $t=90$ ms (center) and $t=250$ ms (right) at differential pressure
$0.5$ bar.

In the remainder of this section, we focus on the proportional operation of
the POCV for a pressure differential that suffices to unseal the main valve,
to examine the vibrations in the POCV in proportional-control mode. To
illustrate the valve response after the main valve opening, Figure 16 plots
the evolution of the difference in the pressure upstream and downstream of the
main valve, and the flow rate across the outlet, during the proportional valve
opening in accordance with Figure 13, for an inlet pressure of $1$ bar. Until
the opening of the pilot valve at $t=0.08$ s, the upstream pressure is
essentially uniform and equal to the specified inlet pressure. Accordingly,
the flow rate vanishes, except for some minute leakage through the valve-seat
gap. The pilot valve is opened to $5\%$ at $t=0.08$ s by suppressing the
auxiliary wall boundary condition near the pilot membrane. The pressure
difference across the main membrane caused by the $5\%$ opening of the pilot
membrane does not suffice to open the main valve. Hence, the upstream-
downstream pressure difference does not change, and the flow across the POCV
remains negligible. The pilot valve displacement at $t=0.12$ s to $25\%$ open
yields an abrupt reduction in the upstream-downstream pressure difference,
followed by a linear decay, and induces a flow of approximately $5$ l/h across
the POCV. At this setting of the pilot valve, the main valve is still closed,
and the flow passes through the pilot conduit; see also Figure 17 below.
Further displacement of the pilot valve at $t=0.14$ s to an approximately
$38\%$ open setpoint leads to a sharp increase in the flow rate from $5$ l/h
to $45$ l/h, and a concurrent sharp reduction of the upstream-downstream
pressure difference from $53\,k$Pa to $19\,k$Pa. These sharp changes in the
flow rate and pressure difference portend severe oscillations in the flow rate
and pressure difference. One can observe that the oscillations exhibit a
consistent frequency, with some modulations. The oscillations occur concurrent
with the opening of the main valve, which indicates that the considered POCV
cannot be used in proportional control mode. For completeness, we mention that
the FSI simulation terminated after $t=0.17752$ s due to failure to converge
in the subsequent time step.

Figure 16: Evolution of the flow rate and upstream-downstream pressure
difference during proportional operation of the POCV for an inlet pressure of
$1\,$bar.

To investigate the oscillatory behavior of the POCV in proportional-control
mode, we regard the flow in the POCV during proportional opening. Figure 17
displays the streamlines at different time instances. The top-left panel
displays the streamlines in the POCV at $t=128.8\,m$s, shortly after opening
the pilot valve to $25\%$. One can note that the streamlines between the
inflow and outflow pass through the pilot channel only, not through the valve-
seat gap, indicating that the main valve is still closed. The remaining panels
in Figure 17 display snapshots of the streamlines during the oscillation of
the POCV. The top-right panel displays the streamlines at $t=160.4\,m$s. This
moment corresponds to a trough in the differential pressure; see Figure 16.
The troughs in the differential pressure coincide with moments at which
contact between the main valve and the valve seat is established, i.e. when
the main valve closes. Simultaneously, a peak in the flow rate occurs. One can
observe that at $t=160.4\,m$s, no streamlines pass from the inlet to the
outlet via the valve-seat gap, in accordance with the fact that the main valve
is closed. The bottom-left panel displays the streamlines at $t=160.8\,m$s. At
this instant, the valve is still closed, the differential pressure is
increasing, and the flow rate is decreasing. The streamlines exhibit a similar
pattern as for $t=160.4\,m$s, in accordance with the fact that the main valve
is still closed. The bottom-right panel displays the streamlines at
$t=161.4\,m$s. In this plot, multiple streamlines pass from the inlet to the
outlet through the valve-seat gap, indicating that the main valve is open.
From the color of the streamlines, one can observe that the flow velocities in
the vicinity of the gap are very large. With reference to Figure 16, we note
that the flow rate exhibits a local minimum at $t=161.4\,m$s.

Figure 17: Streamlines during attempted proportional operation of the POCV at
a differential pressure of $1$ bar at $t=128.8\,m$s (top-left), $t=160.4\,m$s
(top-right), $t=160.8\,m$s (bottom-left) and $t=161.4\,m$s (bottom-right).

Closer inspection of the results in Figures 16 reveals that during the
oscillation, troughs in the differential pressure coincide with peaks in the
flow rate. These extrema correspond to moments at which the main valve closes.
Conversely, pressure peaks occur at moments at which the main valve opens
(results not displayed). It is further noteworthy that pressure peaks do not
coincide with local minima in the flow rate, which indicates that the
oscillation is not simply harmonic. Although it appears contradictory that a
significant outflow can be generated while the main valve is closed, one has
to bear in mind that the main membrane exhibits significant deformations
during the oscillation cycle, including the part of the cycle in which the
valve is closed. Hence, during this part of the cycle, the volume that exits
the POCV is balanced by the motion of the main membrane. Similarly, it is
noteworthy that the flow rate continues to decrease when the main valve opens.
This can be interpreted as that part of the volume that passes through the
gap, is stored in the volume which emerges from the motion of the main
membrane.

A main mechanism underlying the dynamic instability of the POCV in (attempted)
proportional operation, appears to be the significant pressure drop that
occurs in the thin valve-seat gap after opening. At the transition from the
thin valve-seat gap to the outlet channel, in the vicinity of the main
membrane, the flow sustains a strong divergence, which results in a
significant pressure drop, as well as strong vorticity; see Figure 17. The
strong pressure drop in the valve-seat gap counteracts the opening of the
membrane. Moreover, if the motion of the membrane is (locally) reversed due to
the negative pressure in the valve-seat gap, the Venturi effect [56, 20]
becomes operative, causing the valve-seat gap to collapse. The Venturi effect
pertains to a negative pressure singularity ($p\to-\infty$) that occurs in a
fluid between two approaching surfaces; see [20, Thm.29]. Hence, the collapse
of the valve-seat gap occurs forcefully, which is consistent with the strong
oscillations that have been observed in experimental investigations of the
POCV in proportional operation. The above elaboration of the mechanism that is
responsible for the oscillatory response of the POCV in proportional
operation, suggests that an important condition for successful, non-
oscillatory operation of the POCV is that during the opening of the main
valve, the thin-gap regime is transcended sufficiently fast. In the ON/OFF-
operation test cases considered in Section 4.2, the instantaneous complete
opening of the control valve yields a pressure drop across the main membrane
that apparently is sufficiently strong to fulfill this condition.

##  5 Conclusion

Valves play a central role in hydraulic and microfluidic devices. Extending
functionality of valves beyond their original operational range, can provide
an important contribution in miniaturization and cost reduction. From that
perspective, in the present paper we have numerically investigated a smart
pilot-operated control valve (POCV) that has originally been designed for
ON/OFF switching, with the aim of elucidating vibrations that have been
observed experimentally in proportional operation of the device. The numerical
investigations consider the fluid-structure-contact interaction that occurs in
the POCV during operation.

For the fluid-structure-interaction simulations, we have developed a numerical
model in the ANSYS framework, with the auxiliary aim of testing to what extend
such a generic simulation environment is suitable for simulating a complex
industrial FSI problem. On account of the incompressibility of the considered
fluid and the relatively high density of the fluid, the considered POCV
problem represents a strongly coupled FSI problem which, in addition, displays
characteristics of the incompressibility dilemma, because in the closed
configuration, the upstream volume is essentially closed. The so-called
solution-stabilization technique provided by ANSYS-Fluent, which is similar to
artificial compressiblity, allows to stabilize the iterations in the
partitioned solution strategy.

The numerical results convey that the valve exhibits reliable sealing
properties and allows fast switching from open to closed positions. The
simulations show that the considered POCV is reliable in ON/OFF control
operations. The corresponding numerical results display good agreement to
corresponding experimental data, in terms of flow rate and pressure
difference. Investigation of the proportional-operation regime, conveys that
the POCV is unsuitable for proportional control, and is susceptible to
vibrations in this regime. Self-excited vibrations occur as soon as the main
valve is unsealed, resulting in membrane oscillations and pressure pulsations.
The numerical results suggest that the vibrations are essentially caused by
the sharp pressure drop that occurs in the valve-seat gap, which counteracts
the opening of the main valve, and the ensuing Venturi effect, which causes a
forceful collapse of the valve-seat gap. An important condition for successful
operation of the POCV appears to be that during opening of the main valve, the
thin-gap regime is transcended sufficiently fast.

##  Appendix A Membrane material characterization

###  A.1 Constitutive hyperelastic material model

Experimental data are required for the appropriate calibration of a
hyperelastic constitutive model. In this context, three uni-axial tensile
tests were performed in the DICAr department at the University of Pavia [57].
The samples’ dimensions are shown in the table 3. The MTS Insight 10 kN
testing device was used to load the specimen with a $250$ N tensile load.

Table 3: The specimen dimension in terms of length, width, thickness, and Cross section [mm].  Length $l_{0}$ |  Width $W$ |  Thickness $t$ |  Cross section $A_{0}$
---|---|---|---
$45$ | $5$ | $0.45$ | 2.25

The obtained experimental plots shown in Figure 18 depict a quasi-linear
stress-strain mechanical response. Consequently, calibrating the uni-axial
tests using the incompressible neo-Hookean hyperelastic material model is
sufficient for data fitting, for which the constitutive model reads:

$\sigma_{h}=\mu_{s}(\lambda^{2}-\frac{1}{\lambda})$ (25)

with $\sigma_{h}$ is the neo-Hookean engineering stress, and $\mu_{s}=2C_{10}$
the initial shear modulus. For an in-depth overview of hyperelastic material
modeling, the interested reader is referred to [58]. The detailed steps
describing the material calibration using the neo-Hookean hyperelastic
material model are given in [59, 60].

###  A.2 Material calibration and curve fitting

The membrane’s material properties play a major role in the numerical analysis
of the POCV as the nonlinearity and the expected large deformations emanating
from the incoming pressure require a constitutive model to enable accurate
replication of the membrane dynamics. As previously mentioned in A.1, the test
results shown in Figure 18 show a quasi-linear trend implying that the neo-
Hookean material model would be a convenient choice to calibrate the rubber-
like material under study. ANSYS (ANSYS, Canonsburg, PA, USA) supports
material calibration using the neo-Hookean material model as well as various
other conventional hyperelastic models, such as the Mooney-Rivlin and Ogden
models, among others [61]. A curve-fitting approach is used for material
calibration, which incorporates a nonlinear regression analysis for the
experimental data as a function of the constitutive hyperelastic model and is
described as follows:

  1. 1.

Possession of experimental stress-strain data

  2. 2.

Prior knowledge about the structural deformation range expected

  3. 3.

Choosing the appropriate hyperelastic material model

  4. 4.

Performing a regression analysis using the least squares method

  5. 5.

Approximating the material constants associated with the hyperelastic material
model from the experimental stress-strain data and constitutive equations

  6. 6.

Comparing the experimental data against the obtained results. If the error is
large, a more accurate hyperelastic material model could be considered.

The resulting neo-Hookean stress-strain curve is plotted in Figure 18,
alongside the experimental data. Note that the curve fitting procedure is
executed taking into account only sample 1 specimen for simplicity purposes.
To further motivate the model selection, the stress versus strain graphs
corresponding to two commonly used hyperelastic material models are fitted:
The Mooney-Rivlin model with two ($2$) parameters and the Mooney-Rivlin model
with five ($5$) parameters. The stress-strain plot based on the neo-Hookean
model shows good alignment with the experimental data when the specimen’s
cross-head displacement corresponds to engineering strain values below
$100\%$. However, for larger tension strains, higher relative errors are
observed. This is expected, as the neo-Hookean constitutive model is better
suited for relatively low strain ranges [58]. Conversely, the results
associated with the Money-Rivlin models are more accurate for larger strain
ranges. Nevertheless, since the membrane deformations are expected to be below
$100\%$ strain ranges, the neo-Hookean model is chosen since it is the
simplest hyperelastic constitutive model available and is suitable for
capturing the material behavior with sufficient accuracy [62].

Figure 18: Curve-fitting results of the Hyperelastic material calibration.
Stress versus strain curve plots for experimental data versus different
hyperelastic material models. Red: Curve fitting using the neo-Hookean model.
Violet: Curve fitting using the Mooney-Rivlin model with $2$ parameters
material model. Cyan: Curve fitting using the Mooney-Rivlin model with $5$
parameters material model. Black: Experimental results.

###  A.3 Subsystems properties

The material calibration procedure gives the evaluated parameters in terms of
the initial shear modulus $\mu_{0}$ and the incompressibility parameter
$D_{0}$ with corresponding values equal to $0.946$ MPa, and zero
(incompressible material), respectively. The material density is measured
experimentally as $1140$ kg/m3. Because of material testing limitations, the
material-dependent Rayleigh parameters are borrowed from [63]. For
completeness, the material parameters used in the numerical investigations are
summarized in Table 4.

Table 4: Material properties of the fluid and solid used in the simulations.
The main and pilot membranes have identical material properties. The contact
dummy body material properties are consistent with the material of the POCV
casing. The spring is assumed mass-less. The Rayleigh parameters
$\alpha,\beta$ are taken from [63]. Fluid (Water)
---
Density |  $1140$ kg/m3
Viscosity |  $0.946\times 10^{6}$ Pa
Membranes (Rubber)
Density |  $1140$ Kg/m3
Initial shear stress |  $0.946\times 10^{6}$ Pa
Incompressibility parameter | $0$
$\alpha$ | $0.5$
$\beta$ | $0.005$
Dummy solid (Plastic PVC)
Density |  $1392$ kg/m3
Shear Modulus |  $1.02\times 10^{9}$ Pa
Poisson’s ratio | $0.4$
Young’s Modulus |  $2.86\times 10^{9}$ Pa
Bulk Modulus |  $4.78\times 10^{9}$ Pa
Spring
Longitudinal stiffness |  $0.1$ N/mm
Free length |  $14$ mm
Longitudinal damping |  $0$ N/mm
Spring length |  $4.868$ mm
Endpoints | Membrane-Fixed reference point

## Acknowledgments

The authors wish to express their gratitude to the European Union Horizon 2020
research and innovation program for funding Ahmed Aissa Berraies’ Ph.D.
fellowship, which is a part of the Marie Skłodowska-Curie ITN-EJD ProTechTion
actions (Grant Number 76463). The authors acknowledge also partial financial
support from the Ministry of Enterprise and Made in Italy and Lombardy Region
through the project “PRotesi innOvaTivE per applicazioni vaScolari ed
ortopedIChe e mediante Additive Manufacturing” - CUP : B19J22002460005, on the
finance Asse I, Azione 1.1.3 PON Businesses and Competitiveness 2014 - 2020.
The authors are grateful for the computing resources provided by the SurfSara
center. Special thanks are given to ANSYS Inc. for their generous contribution
in providing the necessary academic and HPC licenses.

## Data availability

The model and data generated and analyzed during the present investigation are
not publicly available due to commercial sensitivity. However, specific
information could be provided by the corresponding author upon a reasonable
request.

## Conflict of interest

The authors declare that they have no known competing interests or personal
relationships that could have appeared to influence the work reported in this
paper.

## References

  * [1] M. Sarma, J. Lee, S. Ma, S. Li, C. Lu, A diffusion-based microfluidic device for single-cell rna-seq, Lab on a chip 19 1247.
  * [2] L. Pang, J. Ding, X.-X. Liu, Z. Kou, L. Guo, X. Xu, S.-K. Fan, Microfluidics-based single-cell research for intercellular interaction, Frontiers in Cell and Developmental Biology 9 (2021) 2182.
  * [3] K. W. Oh, C. H. Ahn, A review of microvalves, Journal of Micromechanics and Microengineering 16 (5 2006).
  * [4] A. Czechowicz, S. Langbein, Shape Memory Valves: Motivation, Risks, and Potentials, Springer International Publishing, Cham, 2015, pp. 101–112.
  * [5] G. P. Natarajan, S. J. Kim, C. W. Kim, Analysis of membrane behavior of a normally closed microvalve using a fluid-structure interaction model, Micromachines 2017, Vol. 8, Page 355 8 (2017) 355.
  * [6] G. Scalet, Two-way and multiple-way shape memory polymers for soft robotics: An overview.
  * [7] ProTechTion, Protechtion – msca-itn-ejd.
URL https://www.lacan.upc.edu/ProTechTion/

  * [8] M. B. Alberto, F. O. J. Manuel, M. F. Andrés, Numerical methodology for the cfd simulation of diaphragm volumetric pumps, International Journal of Mechanical Sciences 150 (2019) 322–336.
  * [9] A. Beune, Analysis of high-pressure safety valves, Phd thesis 1 (research tu/e / graduation tu/e), Mechanical Engineering (2009).
  * [10] A. Beune, J. G. Kuerten, M. P. van Heumen, Cfd analysis with fluid-structure interaction of opening high-pressure safety valves, Computers and Fluids 64 (2012) 108–116.
  * [11] Z. H. Lin, X. J. Li, Z. J. Jin, J. Y. Qian, Fluid-structure interaction analysis on membrane behavior of a microfluidic passive valve, Membranes 10 (2020) 1–14.
  * [12] S. K. Chimakurthi, S. Reuss, M. Tooley, S. Scampoli, Ansys workbench system coupling: a state-of-the-art computational framework for analyzing multiphysics problems, Engineering with Computers 34 (2018) 385–411.
  * [13] G. P. Natarajan, S.-J. Kim, C.-W. Kim, Analysis of membrane behavior of a normally closed microvalve using a fluid-structure interaction model, Micromachines 8 (12) (2017) 355.
  * [14] W. Wei, H. Jian, Q. Yan, X. Luo, X. Wu, Nonlinear modeling and stability analysis of a pilot-operated valve-control hydraulic system, Advances in Mechanical Engineering 10 (2018) 1–8.
  * [15] T. C. Allison, K. Brun, Testing and modeling of an acoustic instability in pilot-operated pressure relief valves (2016).
  * [16] T. Lantela, J. Kajaste, J. Kostamo, M. Pietola, Pilot operated miniature valve with fast response and high flow capacity, River Publishers 15 (2014) 11–18.
  * [17] M. Saeedi, R. Wuchner, K.-U. Bletzinger, Multi-fidelity fluid-structure interaction analysis of a membrane wing, International Journal of Aerospace and Mechanical Engineering 9 (1) (2015) 169 – 176.
  * [18] K. Ogawa, T. Kimura, Hydrodynamic characteristics of a butterfly valve — prediction of torque characteristics, ISA Transactions 34 (1995) 327–333.
  * [19] Y. Liu, L. Lu, K. Zhu, Numerical analysis of the diaphragm valve throttling characteristics, Processes 7 (2019).
  * [20] T. Opstal, van, Numerical methods for inflatables with multiscale geometries, Ph.D. thesis, Mechanical Engineering, proefschrift (Dec. 2013).
  * [21] L. Radtke, M. König, A. Düster, The influence of geometric imperfections in cardiovascular fsi simulations, Computers and Mathematics with Applications 74 (2017) 1675–1689.
  * [22] L. Cai, Y. Wang, H. Gao, X. Ma, G. Zhu, R. Zhang, X. Shen, X. Luo, Some effects of different constitutive laws on fsi simulation for the mitral valve, Scientific Reports 9 (1) (2019) 1–15.
  * [23] C. Michler, S. Hulshoff, E. van Brummelen, R. de Borst, A monolithic approach to fluid–structure interaction, Computers & Fluids 33 (2004) 839–848.
  * [24] T. Richter, A monolithic geometric multigrid solver for fluid-structure interactions in ALE formulation, Int. J. Num. Meth. Engng. 104 (2015) 372–390.
  * [25] S. K. Chimakurthi, S. Reuss, M. Tooley, S. Scampoli, Ansys workbench system coupling: a state-of-the-art computational framework for analyzing multiphysics problems, engineering with Computers 34 (2) (2018) 385–411.
  * [26] E. Van Brummelen, Partitioned iterative solution methods for fluid–structure interaction, International Journal for Numerical Methods in Fluids 65 (1-3) (2011) 3–27.
  * [27] Ansys, engineering simulation software.
URL www.ansys.com/

  * [28] P. Causin, J.-F. Gerbeau, F. Nobile, Added-mass effect in the design of partitioned algorithms for fluid–structure problems, Computer methods in applied mechanics and engineering 194 (42-44) (2005) 4506–4527.
  * [29] E. van Brummelen, Added mass effects of compressible and incompressible flows in fluid-structure interaction, Journal of Applied Mechanics: Transactions of the ASME 76 (2) (2009) 021206–1.
  * [30] C. Farhat, A. Rallu, K. Wang, T. Belytschko, Robust and provably second-order explicit–explicit and implicit–explicit staggered time-integrators for highly non-linear compressible fluid–structure interaction problems, International Journal for Numerical Methods in Engineering 84 (1) (2010) 73–107.
  * [31] M. Ezkurra, J. A. Esnaola, M. Martinez-Agirre, U. Etxeberria, U. Lertxundi, L. Colomo, M. Begiristain, I. Zurutuza, Analysis of one-way and two-way fsi approaches to characterise the flow regime and the mechanical behaviour during closing manoeuvring operation of a butterfly valve, International Journal of Mechanical and Materials Engineering 12 (2018) 409–415.
  * [32] J. Degroote, A. Swillens, P. Bruggeman, R. Haelterman, P. Segers, J. Vierendeels, Simulation of fluid–structure interaction with the interface artificial compressibility method, International Journal for Numerical Methods in Biomedical Engineering 26 (3-4) (2010) 276–289.
  * [33] D. P. Mok, Partitionierte lösungsansätze in der strukturdynamik und der fluid-struktur-interaktion.
  * [34] J. Degroote, K.-J. Bathe, J. Vierendeels, Performance of a new partitioned procedure versus a monolithic procedure in fluid–structure interaction, Computers & Structures 87 (11) (2009) 793–801, fifth MIT Conference on Computational Fluid and Solid Mechanics.
  * [35] R. Haelterman, A. Bogaers, K. Scheufele, B. Uekermann, M. Mehl, Improving the performance of the partitioned QN-ILS procedure for fluid–structure interaction problems: Filtering, Computers & Structures 171 (2016) 9–17.
  * [36] U. Küttler, C. Förster, W. A. Wall, A solution for the incompressibility dilemma in partitioned fluid–structure interaction with pure dirichlet fluid domains, Computational Mechanics 2006 38:4 38 (2006) 417–429.
  * [37] S. Badia, F. Nobile, C. Vergara, Fluid-structure partitioned procedures based on robin transmission conditions, Journal of Computational Physics 227 (2008) 7027–7051.
  * [38] S. Badia, F. Nobile, C. Vergara, Robin-robin preconditioned krylov methods for fluid-structure interaction problems, Computer Methods in Applied Mechanics and Engineering 198 (2009) 2768–2784.
  * [39] M. A. Fernández, J. Mullaert, M. Vidrascu, Explicit robin–neumann schemes for the coupling of incompressible fluids with thin-walled structures, Computer Methods in Applied Mechanics and Engineering 267 (0) (2013) 566–593.
  * [40] A. E. Bogaers, S. Kok, B. D. Reddy, T. Franz, An evaluation of quasi-newton methods for application to fsi problems involving free surface flow and solid body contact, Computers & Structures 173 (2016) 71–83.
  * [41] A. L. Spina, C. Förster, M. Kronbichler, W. A. Wall, On the role of (weak) compressibility for fluid-structure interaction solvers, International Journal for Numerical Methods in Fluids 92 (2020) 129–147.
  * [42] C. Ager, A. Seitz, W. Wall, A consistent and versatile computational approach for general fluid-structure-contact interaction problems, Int. J. Numer. Methods Eng. 122 (19) (2021) 5279–5312.
  * [43] F. R. Menter, F. R. Menter, Aiaa 93·2906 zonal two equation k·co turbulence models for aerodynamic flows. zonal two equation k-w turbulence models for aerodynamic flows (1993).
  * [44] F. R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA journal 32 (8) (1994) 1598–1605.
  * [45] E. H. van Brummelen, P. Geuzaine, Fundamentals of fluid-structure interaction, Encyclopedia of aerospace engineering (2010) 153–1.
  * [46] J. Donea, S. Giuliani, J. P. Halleux, An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions, Computer Methods in Applied Mechanics and Engineering 33 (1982) 689–723.
  * [47] C. Farhat, Cfd-based nonlinear computational aeroelasticity, Encyclopedia of computational mechanics (2004).
  * [48] Y. Bazilevs, V. M. Calo, T. J. R. Hughes, Y. Zhang, Isogeometric fluid-structure interaction: theory, algorithms, and computations, Comput. Mech. 43 (2008) 3–37.
  * [49] V. A. Yastrebov, Numerical methods in contact mechanics, John Wiley & Sons, 2013\.
  * [50] A. Popp, Mortar methods for computational contact mechanics and general interface problems, Ph.D. thesis, Technische Universität München (2012).
  * [51] S. Piperno, Simulation numérique de phénomènes d’interaction fluide-structure.
  * [52] C. Förster, W. Wall, E. Ramm, Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows, Comput. Methods Appl. Mech. Engrg. 196 (2007) 1278–1293.
  * [53] A. E. Bogaers, S. Kok, B. D. Reddy, T. Franz, Extending the robustness and efficiency of artificial compressibility for partitioned fluid-structure interactions, Computer Methods in Applied Mechanics and Engineering 283 (2015) 1278–1295.
  * [54] T. van Opstal, E. van Brummelen, G. van Zwieten, A finite-element/boundary-element method for three-dimensional, large-displacement fluid–structure-interaction, Computer Methods in Applied Mechanics and Engineering 284 (2015) 637–663.
  * [55] C. Warnett, Valve actuators (2015).
  * [56] T. van Opstal, E. van Brummelen, A finite-element/boundary-element method for large-displacement fluid-structure interaction with potential flow, Comput. Methods Appl. Mech. Engrg. 266 (2013) 57–69.
  * [57] Computational mechanics and advanced materials group, università degli studi di pavia, unipv, www.compmech.unipv.it.
  * [58] G. A. Holzapfel, J. Wiley, Nonlinear solid mechanics: A continuum approach for engineering science, Meccanica 2002 37:4 37 (2002) 489–490.
  * [59] I. Doghri, Mechanics of deformable solids: linear, nonlinear, analytical and computational aspects, Springer Science & Business Media, 2000.
  * [60] P. Martins, R. Natal Jorge, A. Ferreira, A comparative study of several material models for prediction of hyperelastic properties: Application to silicone-rubber and soft tissues, Strain 42 (3) (2006) 135–147.
  * [61] P. A. Martins, R. M. Jorge, A. J. Ferreira, A comparative study of several material models for prediction of hyperelastic properties: Application to silicone-rubber and soft tissues, Strain 42 (2006) 135–147.
  * [62] L. A. Mihai, A. Goriely, How to characterize a nonlinear elastic material? a review on nonlinear constitutive parameters in isotropic finite elasticity, Proceedings. Mathematical, Physical, and Engineering Sciences 473 (11 2017).
  * [63] R. K. Luo, Numerical prediction and experiment of impact on rubber hydro unit with multi-chambers using an effective dynamic co-simulation approach, Vehicle System Dynamics 58 (2020) 249–271.