# Multiscale Modeling with Differential Equations

Clarissa Astuto, Giovanni Russo

###### Contents

  1. 1 Multiscale Modeling with Differential Equations
    1. 1.1 Introduction and Examples
    2. 1.2 Time multiscale in a damped oscillator
      1. 1.2.1 Dimensional analysis
      2. 1.2.2 Singular perturbation problems and differential-algebraic equations
    3. 1.3 Coupling conditions for sorption kinetics at the oscillating cell membrane
      1. 1.3.1 1D Multiscale model
      2. 1.3.2  Extension to more dimensions
      3. 1.3.3 Validation of the model and results
    4. 1.4 Coupling conditions for cardiovascular system
      1. 1.4.1 The 3D model
      2. 1.4.2 The 1D and 0D model
      3. 1.4.3 Coupling of 3D-1D, 3D-0D and 1D-0D models
    5. 1.5 Coupling conditions for shallow water equations at channel junctions
      1. 1.5.1 Mathematical model
      2. 1.5.2 Methodology
    6. 1.6 Appendix
      1. 1.6.1 Shallow water equations
      2. 1.6.2 Laplace-Beltrami operator

###### List of Figures

  1. 1.1 This figure shows that different sizes of venations are present in the same leaf.
  2. 1.2 Description of Saint-Venant problem.
  3. 1.3 Schematic setup of the real apparatus. The central sphere mimics the oscillating gas bubble. The detectors (in red) are located at distance P from the bottom of the vessel.
  4. 1.4 Representation of large arteries and thin capillaries in the circulatory system.
  5. 1.5 Shallow water equation is an approximation of NS equations, when the horizontal length scale is much greater than the vertical length scale, i.e. $\eta+h\ll$ x-direction of the flow.
  6. 1.6 Example of a channel network. Regions with two dimensional behaviour are zoomed in.
  7. 1.7 Schematic representation of system (1.38-1.39).
  8. 1.8 In these tree plots we see different oscillations for different elastic constants, $k$, namely $k=10^{2},10^{3},10^{3}$ (N/m). When $k\to\infty$, the oscillator behaves like a classical pendulum.
  9. 1.9 Scheme of surfactant disposition. They are composed by hydrophilic heads, that prefer to stay in the water, and and hydrophobic tails.
  10. 1.10 Representation of the external potentials $V(x)$ (in Eq. (1.49)), on the left, and $U(\xi)$ (in Eq. (1.50)), on the right, for $\varepsilon=0.05$. $U(\xi)$ is obtained from a change of variable of $V(x)$, with $\xi\in[0,L+1]$.
  11. 1.11 Scheme of the boundary condition in 2D where $\hat{n}$ is the outgoing normal vector to the domain $\Omega$, and $\varepsilon=\varepsilon(S)$ is the thickness parameter of the bubble surface $S$. The internal line near the bubble surface represents the boundary $\Gamma_{B}$.
  12. 1.12 Here there are two aspects of the same plot. On the left it is shown the effect of the attractive potential that becomes higher decreasing the value of $\varepsilon$. On the right there is the comparison between $c^{\varepsilon}$, and the dashed red line that represents the solution of the multiscale model, $c^{0}$.
  13. 1.13 The effect of the attractive potential in a 2D domain. Here the bubble has a reasonable circular shape. As before, in the first three panels (a,b,c) there are the solution $c^{\varepsilon}$ for three different values of $\varepsilon$. In panel (d) there is a 1D plot, that is a section of the previous ones for $y=0$ with an additional line. The black dashed line represents the solution of the multiscale model, $c^{0}$.
  14. 1.14 An example of a chopped domain. On the left, the purple object represents the time dependent fluid domain, while the red one represents the structure domain.
  15. 1.15 Classification of Newtonian and Non-Newtonian fluids
  16. 1.16 Systemic arterial tree consisting of larger arteries originates at the heart (A) and terminates at the points marked with (B) and structured trees representing smaller arteries originate at these terminals and provide main tree with outflow boundary conditions.
  17. 1.17 Scheme of domain splitting (the subscript ${\rm f}$ stands for fluid and ${\rm s}$ for solid): (1) $\Omega$ is the initial considered domain and $\Sigma$ the related boundary. (2) $\Omega_{\rm f}\cup\Omega_{\rm s}$ is a partition of $\Omega$ and $\Gamma$ is the common interface. (3) The 3D domain $\Omega_{\rm s}$ is replaced with a 1D domain $[0,L]$ in $z$-direction.
  18. 1.18 Coupling 2D and 1D domain.
  19. 1.19 (a): Plot of the functional defined in Eq. (1.129). (b): Plot of the functional defined in Eqs.(1.133-1.135).
  20. 1.20 (a): Scheme of the level-set function. (b): Scheme of a potential V(x) that represents a squared sink.

##  Chapter 1 Multiscale Modeling with Differential Equations

Prerequisites

  * •

Differential Equations (DE)

Learning Objectives

  * •

Learning about mathematical modelling

  * •

Describe mathematical modeling

  * •

Show how mathematical techniques can

  * •

Reduce the complexity of a system

  * •

Define multiple scales

  * •

Show how the boundary conditions are relevant in a PDE problem

  * •

Analyze and compare 3D, 2D, and 1D models

Learning Outcomes

Students will learn:

  * •

How to solve first order ODE, given the initial condition

  * •

How to use complex exponentials

  * •

About techniques based on asymptotic expansions

  * •

How to use references to build up new knowledge

  * •

How to practise new knowledge solving the proposed exercises and problems

### Contextual Introduction - Summary

Many physical systems are governed by ordinary or partial differential
equations (see, for example, Chapter ”Differential equations”, ”System of
Differential Equations”). Typically the solution of such systems are functions
of time or of a single space variable (in the case of ODE’s), or they depend
on multidimensional space coordinates or on space and time (in the case of
PDE’s).

In some cases, the solutions may depend on several time or space scales. An
example govered by ODE’s is the damped harmonic oscillator, in the two extreme
cases of very small or very large damping (see Section 1.3), the
cardiovascular system, where the thickness of the arteries and veins varies
from centimeters to microns, shallow water equations, which are valid when
water depth is small compared to typical wavelength of surface waves, and
sorption kinetics, in which the range of interaction of a surfactant with an
air bubble is much smaller than the size of the bubble itself. In all such
cases a detailed simulation of the models which resolves all space or time
scales is often inefficient or intractable, and usually even unnecessary to
provide a reasonable description of the behavior of the system.

In the Chapter ”Multiscale modeling with differential equations” we present
examples of systems described by ODE’s and PDE’s which are intrinsically
multiscale, and illustrate how suitable modeling provide an effective way to
capture the essential behavior of the solutions of such systems without
resolving the small scales.

###  1.1 Introduction and Examples

This chapter focuses on space and time multiscale challenges that physical
problems face in reality.

There are many cases in nature in which problems are characterized by the
appearance of structures, whose elements have sizes which span several orders
of magnitude.

One example that shows up in nature concerns the formation of venation
networks in leaves. The structure of the branches creates a very complicated
reticulum, with different ramifications and several different lengths. The
lengths of the leaf network can vary from $10^{-5}$ $m$, if we consider veins
diameters, to $10^{0}$ $m$, that is the length of some large leaf (see Fig.
1.1), and correspondingly the density of veins dramatically increases as the
diameter gets smaller and smaller, thus making it impractical to adopt the
same model to describe the flow through large and small veins. See for example
[2, 3, 4].

When several spatial scales appear in a physical system, for which it is not
possible or is impractical to resolve all the spatial scales, it is common to
refer to such cases as to space and/or time multiscale problems [1].

Another well-known example is given by pandemic behaviour. It is a phenomenon
that has, not only personal effects, but also, global consequences, such as
economic, social and political. To design a good model to describe a pandemic,
one has to take multiscale levels into consideration, from individuals to
worldwide population, from social distances to social policies, from hand
hygiene to the invention of a vaccine that can be distributed to several
countries.

Figure 1.1: This figure shows that different sizes of venations are present in
the same leaf.

In this chapter, the main focus is on four cases, in which different space and
time multiple scales are involved. The first application is an example of time
multiscale in a damped oscillator. The second topic is the capture rate of
substances freely diffusing around a trap; then we move on with the human
cardiovascular system, and, at the end, we describe how to reduce the
complexity of hydro flows in channels .

Figure 1.2: Description of Saint-Venant problem.

Solution of realistic models are often a challenge even for the most advanced
computer systems. Very often the complexity comes from the multiscale nature
of the processes and the variety of the parameters involved.

Several mathematical models, describing a large range of physical systems, are
described by a set of ordinary differential equations (ODEs) or partial
differential equations (PDEs), and then solved by suitable numerical methods.

There is a vast literature on mathematical models described by differential
equations, and an equally vast literature on the various techniques for their
numerical treatment [8, 26, 25].

Computational power of modern computers is not always sufficient to provide
accurate solutions, especially, in three space dimensions, because of the huge
number of unknowns that appears in the problem when multiple scales are
involved. Therefore, it is necessary to construct suitable multiscale
mathematical models of the physical system, that take advantage of the
specific nature of the problem, thus allowing a drastic reduction of the
number of degrees of freedom (DOFs), but still capturing the essential
features of the phenomena that one is interested in.

Finding the right boundary conditions can be very helpful in reducing the DOF
of a system, when the geometry allows it.

In 1855 Adhémar Barré de Saint-Venant, a French mathematician and engineer,
who gave fundamental contribution in the theory of linear elasticity,
introduced the following strategy, known as de Saint Venant’s principle to
reduce the complexity of a problem in linear elasticity [6] with free sides.
He ingeniously sidestepped the major difficulty of his problem by replacing it
with a simplified one: consider a cylinder composed of homogeneous isotropic
elastic material, with free sides, and in the absence of gravity. The cross
sections are arbitrary and, to satisfy static equilibrium, the needed loads
are applied, which cause the cylinder to deform, ending up with the problem of
de Saint Venant in which the goal is to calculate the stress and the
displacement throughout the cylinder. According to de Saint Venant principle,
the boundary conditions are relaxed, i.e., they are replaced by simpler
conditions at the end side of the body, which share with the original one the
same total force and torque. In this way, three out of the six components of
stress are set to zero, thus drastically simplifying the problem [7] (see Fig.
1.2). When considering torsion or flexion of the cylinder, the problem turns
into a boundary condition issue, and indicates how to choose the proper
conditions to impose (see, for instance, [7]).

This chapter will first show an example within the context of ODE’s (see
Section 1.2). We show a simple case in which multiple scales are involved in
time, and the parameters of the system vary in a very wide range. We describe
two different time scales, $t$, related to the time evolution, and $\tau$,
related to the undamped oscillator. By examining these two different time
scales, we can gain a deeper understanding of how the undamped oscillator
behaves and what its characteristics are. Moreover, comprehending how these
time scales interact with each other enables us to analyze and make
predictions about how the oscillator will behave under various conditions.
Some of the applications that could be studied are fascinating phenomena, such
as resonance, frequency locking, and other intriguing dynamical effects.

The second application will first focus on sorption kinetics application, in
which it is fundamental to couple systems with different dimensions, such as
modeling the capture rate of particles (3D) by moving traps (2D). This is a
common research topic among different fields of sciences, such as in
chemistry, when studying biomolecular reactions, or in biology, with complex
prey-predator ecological systems. The last example is given by living cells
and their moving membranes. The nature of this motion is associated to the
exchange of thermal and metabolic energy [13, 14, 16, 17] and it has been
studied how the oscillations of the cell membrane affect the capture rate of
substances freely diffusing around them. To obtain quantitative results, a
reproducible and tunable biomimetic experimental setup has been considered to
simulate the phenomenon. The system consists of an oscillating air bubble in
water and a diffusive flow of ions, see Fig. 1.3. When the negative ions meet
the surface of the bubble, they are partially adsorbed because of their polar
head, but at the same time, they also prefer to settle at the bubble surface
because of the hydrophobic repulsion between apolar tails and the water where
it is immersed (see Fig. 1.9). This sort of substance is called surfactant,
and in Section 1.9 we see it in details.

\begin{overpic}[abs,width=216.81pt,unit=1mm,scale={.25}]{Figures_new/setup_delta.pdf}
\put(21.0,14.0){{P}} \put(10.0,38.0){\footnotesize hydrophobic tail}
\put(16.0,36.0){\footnotesize hydrophilic} \put(18.0,33.0){\footnotesize head}
\end{overpic} Figure 1.3: Schematic setup of the real apparatus. The central
sphere mimics the oscillating gas bubble. The detectors (in red) are located
at distance P from the bottom of the vessel.

Modeling diffusion and trapping in presence of a bubble is very challenging
because of the multiple scales involved in space. The range of the attractive-
repulsive core of a trap is of the order of tens of nanometers ($10^{-9}m$), a
length that is prohibitively small for detailed simulations (see , e.g.,
[33]). while the size of the domain is a few millimeters ($10^{-3}m$), thus
almost 6 different orders of magnitude in space.

To overcome such difficulties, in [15] the authors propose a _multi-scale
model_ to effectively treat the problem describing the trapping by a suitable
boundary condition on the bubble surface. The potential is assumed to have a
range of small size $\varepsilon$. An asymptotic expansion in $\varepsilon$ is
considered, and the boundary conditions are obtained by retaining the lowest
order terms in the expansion. In the reduced multiscale model, the ions are
described by a volumetric concentration which satisfies a simple diffusion
equation in the domain of the liquid (3D), while the carriers that are trapped
on the surface satisfy an equation for the superficial concentration (2D),
that is coupled with the parabolic one. The two concentrations are related by
a mass balance equation which acts as a suitable boundary condition for the
concentration in the bulk. This model is described in detail in Section 1.9.

In Section 1.4 we refer to geometric multiscale approach for the numerical
simulation of blood flow problems, a very hot topic in applied mathematics
(see, e.g., [9, 11]). The geometric multiscale analysis refers to the idea of
coupling different dimensions in order to obtain a detailed study of the
phenomenon. It plays an important role when different scales or levels
interact in a system. With this approach, different features of the entire
cardiovascular system can described, focusing on the right scale that is more
suitable for the problem.

Figure 1.4: Representation of large arteries and thin capillaries in the
circulatory system.

Human circulatory system is well approximated by hyperbolic equations. Even if
blood can be considered an incompressible fluid, the elasticity of veins and
arteries introduces compressibility in the model, with finite speed of
propagaton, which is much smaller than the propagation of sound waves in
water. The first one who developed a mathematical model for liquid flow in an
elastic duct was Leonard Euler, as appears in a paper published after his
death [24]111It is conceivable that he wrote the paper even before his
celebrated papers on Euler equations of compressible gas dynamics, datin 1757
[23]. In this respect he may be considered the father of hemodynamics (as well
as of innumerable other mathematical discoveries!)

Let us restrict, for example, to the pulmonary blood circulatory system. The
vessels, have a wide range of different diameters. Pulmonary arterial and vein
trees are composed by many orders of branches, with diameters ranging from few
$\rm\mu$m for capillaries to about 2-3 cm for the main pulmonary artery [10].

This is what is called space multiscale challenges. Trying to solve so many
different orders of magnitude with the same space discretization, would be
very computationally expensive and, in certain cases, not effective, unless
all space scales are very well resolved. Mathematical modeling and
approximation are needed to investigate such cases.

Figure 1.5: Shallow water equation is an approximation of NS equations, when
the horizontal length scale is much greater than the vertical length scale,
i.e. $\eta+h\ll$ x-direction of the flow.

Going back to the idea of reducing the number of DOF, a common technique in
such cases consists in adopting different spatial models in different parts of
the system, and coupling them together. In particular, in regions where high
accuracy is not required, one can adopt a lower dimensional model.
Practically, the geometry of the problem reveals the region where the study is
focused on, describing it with a 3D model. The description of the background
is modeled with lower dimensional problems. As a result, the problem is solved
coupling 3D - 1D, 3D - 0D or 1D - 0D systems. A three-dimensional (3D) model
can be reduced to a one-dimensional (1D) model, or to a zero-dimensional (0D)
model, where only the dependence on time survives.

Matching domains with different dimensions can generate inconsistency at the
boundaries. For this reason, the boundary conditions play an important role in
the system. In particular, they play a key role in determining suitable
treatment of the junctions, i.e., those parts of the system which separate
domains of different dimensions. A similar problem appears in systems that
treat water channels (Section 1.3) and gas tubes.

In summary, the simulation of the entire cardiovascular system of human body
is a challenging multiscale problem. Adopting a multiscale approach, by
coupling domains of different size using different spatial dimensions, allows
numerical simulations of cardiovascular diseases, thus describing many
features of medical interest with progressively more accuracy and precise
quantification of the error.

The last example of space scales reduction is shown in Section 1.5. The
coupling of different dimensions is also widely used to describe flows in
artificial channels with applications to environmental problems [12, 19, 20],
where the basic model comes from the Shallow Water (SW) equations. These are
approximations of incompressible, irrotational Euler equation with free
boundary, valid when the depth of the water is small compared to the wave-
length. Under this assumption it can be shown that the pressure is basic the
hydrostatic pressure, i.e. in each point it is proportional to the height of
the column of water above that point. Section 1.6.1 in the Appendix is
dedicated to the extended derivation of the SW equations. In the 1D case such
equations are also called de Saint-Venant equations (yes, the same de Saint
Venant of linear elasticity!) [5].

Figure 1.6: Example of a channel network. Regions with two dimensional
behaviour are zoomed in.

In nature, many phenomena are described by SW equations, such as flows in open
channels and rivers. It has been proposed a 2D - 1D model for the junction of
shallow-water channel networks, with a localized 2D problem at the junctions
and a 1D model for the channels, because it is much more efficient from a
computational point of view. The difficulty is in finding the right expression
for the conditions at the junction between the channels; see Fig. 1.6.

###  1.2 Time multiscale in a damped oscillator

Multiscale problems can appear in time as well, when very different time
scales are involved in a phenomenon.

Let us consider again the example of the damped oscillator, seen in Chapter
Dimensional Analysis. The equations governing its motion are given by:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=v$ (1.1) $\displaystyle
m\frac{dv}{dt}$ $\displaystyle=-kx-\gamma v$ (1.2)

where $m,k$ and $\gamma$ are, respectively, the mass of the oscillator, the
elastic constant of the spring and the damping coefficient.

Eqs. (1.1-1.2) are a system of two differential equations of first order.
Differentiating the first equation (1.1) and making use of the second one
(1.2), the dynamics can be described by a differential equation of second
order:

$\frac{d^{2}x}{dt^{2}}=-\frac{\gamma}{m}\frac{dx}{dt}-\frac{k}{m}x.$ (1.3)

Eq. (1.3) is an example of ordinary differential equation of degree two, with
constant coefficients.

If $k=0$, there is no spring, and the equation for $v$ reduces to

$\frac{dv}{dt}=-\frac{\gamma}{m}v.$ (1.4)

The solution can be found by observing that the time derivative of the
function $v$ is proportional to itself. The functions that are proportional to
their own derivatives are the exponentials, therefore it is natural to look
for a solution of the form

$v=e^{\lambda t}.$ (1.5)

By inserting Eq. (1.5) into Eq. (1.4), one obtains

$\lambda e^{\lambda t}=-\frac{\gamma}{m}e^{\lambda t}$ (1.6)

from which we can deduce that

$\lambda=-\frac{\gamma}{m}$ (1.7)

thus, one solution of Eq. (1.4) is

$\displaystyle\bar{v}(t)=e^{-\frac{\gamma}{m}t}.$ (1.8)

Since the equation for $v$ is linear, if a given function $\bar{v}(t)$ is a
solution, then $C\bar{v}(t)$ is also a solution, with $C$ is an arbitrary
constant. If we impose that such a function satisfies the initial condition
$v(0)=v_{0}$, then this will fix the constant $C$, such that $C=v_{0}$, and
the solution will be

$v(t)=v_{0}e^{-\frac{\gamma}{m}t}.$ (1.9)

From Eq. (1.1) it is possible to compute the position of the mass by
integrating the velocity:

$\displaystyle x(t)$
$\displaystyle=x_{0}+\int_{0}^{t}v(t^{\prime})dt^{\prime}$
$\displaystyle=x_{0}+v_{0}\int_{0}^{t}e^{-\frac{\gamma}{m}t^{\prime}}dt^{\prime}$
$\displaystyle=x_{0}-v_{0}\frac{m}{\gamma}\left[e^{-\frac{\gamma}{m}t^{\prime}}\right]_{0}^{t}$
$\displaystyle=x_{0}-v_{0}\frac{m}{\gamma}\left(e^{-\frac{\gamma}{m}t}-1\right),$
(1.10)

where $x_{0}=x(0)$ and, if $t\to\infty$, then $x\to x_{\rm
f}=x_{0}+\frac{v_{0}\,m}{\gamma}$.

Notice that, dimensionally,
$\displaystyle\left[\frac{\gamma}{m}\right]=T^{-1}$, so
$\displaystyle\frac{m}{\gamma}$ is the decay time due to the damping, which we
denote by $t_{d}=\frac{m}{\gamma}$.

This result shows that, with no spring a mass would slow down to stop at final
position $x_{\rm f}$, starting at $x_{0}$ with initial velocity $v_{0}$.

Now let us consider the case in which $\gamma=0$. In such cases the equation
for the position $x(t)$ would be

$\frac{d^{2}x}{dt^{2}}=-\frac{k}{m}x.$ (1.11)

This equation tells us that the function $x(t)$ is such that it is
proportional to its second derivative, but with a different sign. What are the
functions whose second derivative are proportional to the opposite of
themselves? Those are the trigonometric functions! So we look for solutions of
the form

$x(t)=\sin(\omega t)\quad{\rm or}\quad x(t)=\cos(\omega t).$ (1.12)

In the first case we have

$\displaystyle x^{\prime}(t)$ $\displaystyle=\omega\cos(\omega t),$ (1.13)
$\displaystyle x^{\prime\prime}(t)$ $\displaystyle=-\omega^{2}\sin(\omega
t)=-\frac{k}{m}x(t)=-\frac{k}{m}\sin(\omega t)$ (1.14)

indicating that $\sin(\omega t)$ is a solution if $\omega^{2}={k}/{m}$.
Likewise, in the second case we have that $x(t)=\cos(\omega t)$ is a solution,
still providing $\omega^{2}={k}/{m}$.

Since Eq. (1.11) is linear, if $\sin(\omega t)$ and $\cos(\omega t)$ are
solutions, then any linear combination is a solution. The general solution is
therefore

$\displaystyle x(t)$ $\displaystyle=a\cos(\omega t)+b\sin(\omega t)$ (1.15)
$\displaystyle v(t)$ $\displaystyle=-a\,\omega\sin(\omega
t)+b\,\omega\cos(\omega t).$ (1.16)

The coefficients $a$ and $b$ can be determined by imposing the initial
conditions

$x(0)=x_{0},\quad v(0)=v_{0}.$ (1.17)

Observe that the expression of the two solution is not so different if we make
use of complex exponentials. Indeed, since we have Euler formula

$e^{i\omega z}=\cos(\omega t)+i\sin(\omega t)$ (1.18)

and its inverse:

$\cos(\omega t)=\frac{e^{i\omega t}+e^{-i\omega t}}{2},\quad\sin(\omega
t)=\frac{e^{i\omega t}-e^{-i\omega t}}{2i}$ (1.19)

where $i$ denotes the imaginary units, i.e., $i^{2}=-1$, then the solution of
the undamped harmonic oscillator can be written in terms of complex
exponentials:

$x(t)=C_{+}e^{i\omega t}+C_{-}e^{-i\omega
t},\qquad\omega=\sqrt{\frac{k}{m}}=\frac{1}{t_{0}}.$ (1.20)

The time constant $t_{0}=\omega^{-1}$ is a related to a full period $T$ of the
undamped oscillator by $T=2\pi t_{0}$, because after one period the argument
has changed by $\omega T=2\pi$.

Here we focus on the general case, i.e., $k\neq 0,\gamma\neq 0$. Eq. (1.3) for
the position has the following structure

$x^{\prime\prime}+\widetilde{\gamma}x^{\prime}+\omega^{2}x=0,\qquad\widetilde{\gamma}=\frac{\gamma}{m}=\frac{1}{t_{d}}.$
(1.21)

where for short we denote by $\leavevmode\nobreak\ {}^{\prime}$ the time
derivative. It is therefore natural to look for solutions of the form
$x(t)=e^{\lambda t}$, and inserting this expression in Eq. (1.21) we obtain
that $e^{\lambda t}$ is indeed a solution provided

$\lambda^{2}+\widetilde{\gamma}\lambda+\omega^{2}=0.$ (1.22)

This is an algebraic equation of degree 2. The discriminant is given by
$\Delta=\widetilde{\gamma}^{2}-4\omega^{2}$. If $\Delta>0$ there will be two
real negative roots, while if $\Delta<0$ there will be two complex conjugate
roots.

If $\Delta<0$, i.e. $\widetilde{\gamma}<2\omega$ or, equivalently, if
$t_{d}>2t_{0}$, the solution shows damped oscillations, if $\Delta>0$, i.e. if
$t_{d}<t_{0}$, the solution is damped with no oscillations, while if
$\Delta=0$ we have the so called critical damping222For a given fixed
oscillation frequency $\omega=\sqrt{k/m}$, the value of damping
$\widetilde{\gamma}^{*}$, that gives $\Delta=0$, is the one which has the
fastest damping..

####  1.2.1 Dimensional analysis

The Eqs. (1.1-1.2), governing the damped oscillator, apparently depend on
three parameters, $m,k$ and $\gamma$. However, we observe that after dividing
by the mass $m$, there are only two parameters the equation depends on. If we
use non-dimensional time to describe the evolution of the system, we can show
that there is only one parameter that determines the quantitative behaviour of
the solution.

Indeed, in an undamped oscillator, the quantity $\omega t$ can identify a non
dimensional time $\tau$. In term of this parameter the equations become:

$\displaystyle\omega\frac{dx}{d\tau}$ $\displaystyle=v$ (1.23)
$\displaystyle\omega\frac{dv}{d\tau}$
$\displaystyle=-\omega^{2}x-\frac{\gamma}{k}v.$ (1.24)

Differentiating (1.23), and replacing the expression of $dv/d\tau$ in the
second equation, after dividing by $\omega^{2}$ we obtain the following second
order differential equation for $x$:

$\frac{d^{2}x}{d\tau^{2}}+f\frac{dx}{d\tau}+x=0$ (1.25)

and the discriminant is $\Delta=f^{2}-4$, so the critical damping occurs if
$f=2$. If $f<2$, i.e., $t_{0}<2t_{d}$, there are damped oscillations, and if
$f>2$ there are no oscillations. If $f\to\infty$, the equation for $v$
suggests that $v\to 0$, and $x$ is constant, so the limit of infinite damping
is not very interesting.

If we want to observe a non trivial behaviour for $f\to\infty$, we have to act
differently.

Let us start again from the original system (1.1-1.2), and let us divide the
second equation by $k$, obtaining

$\frac{m}{k}\frac{dv}{dt}=-x-\frac{\gamma}{k}v.$ (1.26)

We can write this equation as

$t_{0}^{2}\frac{dv}{dt}=-x-t_{\gamma}v$ (1.27)

where $\displaystyle t_{0}=\sqrt{\frac{m}{k}}$ and $\displaystyle
t_{\gamma}=\frac{\gamma}{k}$. Notice that $t_{d}\,t_{\gamma}=t_{0}^{2}$.

We now define a new non-dimensional time $\tilde{\tau}$ as
$\displaystyle\tilde{\tau}=\frac{t}{t_{\gamma}}$, with this definition, Eq.
(1.23) becomes:

$\frac{t_{0}^{2}}{t_{\gamma}^{2}}\frac{d^{2}x}{d\tilde{\tau}^{2}}=-x-\frac{dx}{d\tilde{\tau}}.$
(1.28)

As $\displaystyle\frac{t_{0}}{t_{\gamma}}\to 0$ we obtain the equation

$\frac{dx}{d\tilde{\tau}}=-x$ (1.29)

whose solution is $x=x_{0}e^{-\tilde{\tau}}=x_{0}e^{-t/t_{\gamma}}$.

Note that, the limit is different form the previous one. Indeed, the original
second order system (1.1-1.2), which we can rewrite in the form

$\displaystyle\frac{dx}{dt}$ $\displaystyle=v$ (1.30) $\displaystyle
t_{0}^{2}\frac{dv}{dt}$ $\displaystyle=-x-t_{\gamma}v$ (1.31)

will relax to the system

$\displaystyle\frac{dx}{dt}$ $\displaystyle=v$ (1.32) $\displaystyle
x+t_{\gamma}v$ $\displaystyle=0,$ (1.33)

as $\tilde{\tau}_{0}\to 0$, which is a first order system of ODE’s.

The solution of such system requires fewer initial conditions: the initial
condition of the position, $x_{0}$, is sufficient, and the velocity is related
to the position by the algebraic equation $x+t_{\gamma}v=0$.

####  1.2.2 Singular perturbation problems and differential-algebraic
equations

This example is a particular case of a large class of systems in which there
is a small parameter that introduces a small scale. Such systems take the
form:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=f(t,x,y)$ (1.34)
$\displaystyle\varepsilon\frac{dy}{dt}$ $\displaystyle=g(t,x,y)$ (1.35)

where $x(t)\in\mathbb{R}^{n_{1}}$, $y(t)\in\mathbb{R}^{n_{2}}$.

As $\varepsilon\to 0$, the system degenerates into a system of $n_{1}$
differential equations, and $n_{2}$ algebraic equations of the form

$\displaystyle\frac{dx}{dt}$ $\displaystyle=f(t,x,y)$ (1.36) $\displaystyle 0$
$\displaystyle=g(t,x,y).$ (1.37)

This kind of systems are called differential-algebraic systems (DAE) (see,
e.g. [34]).

An example of such system is given by the pendulum obtained by a spring, with
elastic constant $k$ and rest length $L$, fixed on one end, and with a mass
$m$ on the other end (see Fig.1.7). The equations of motion of the coordinates
of the system are given by Newtons second law:

$\displaystyle m\ddot{x}_{1}$ $\displaystyle=-k(x_{1}-L\sin\theta)$ (1.38)
$\displaystyle m\ddot{x}_{2}$ $\displaystyle=-k(x_{2}-L\cos\theta)+mg$ (1.39)

Here

$x_{1}=L\sin\theta,\quad x_{2}=L\cos\theta$ (1.40)

from which we can deduce

$\sin\theta=\frac{x_{1}}{\sqrt{x_{1}^{2}+x_{2}^{2}}},\,\cos\theta=\frac{x_{2}}{\sqrt{x_{1}^{2}+x_{2}^{2}}},$

USing these relations, the equations of motion can be written as:

$\displaystyle m\ddot{x}_{1}$
$\displaystyle=-kx_{1}\left(1-\frac{L}{\sqrt{x_{1}^{2}+x_{2}^{2}}}\right)$
(1.41) $\displaystyle m\ddot{x}_{2}$
$\displaystyle=-kx_{2}\left(1-\frac{L}{\sqrt{x_{1}^{2}+x_{2}^{2}}}\right)+mg$
(1.42)

\begin{overpic}[abs,width=216.81pt,unit=1mm,scale={.25}]{Figures_new/axis_x1_x2.pdf}
\put(17.0,30.0){$\theta$} \put(31.0,20.0){$m$} \put(9.0,21.0){$x_{2}$}
\put(30.0,46.0){$x_{1}$} \end{overpic} Figure 1.7: Schematic representation of
system (1.38-1.39).

As $k\to\infty,$ then $\displaystyle
k\left(1-\frac{L}{\sqrt{x_{1}^{2}+x_{2}^{2}}}\right)$ becomes undetermined,
since $k\to\infty$ and
$\displaystyle\left(1-\frac{L}{\sqrt{x_{1}^{2}+x_{2}^{2}}}\right)\to 0$. Let
$\lambda$ denote such limit. Then the equations of motion become

$m\ddot{x}_{1}=-\lambda x_{1},\quad m\ddot{x}_{2}=-\lambda x_{2}+mg$ (1.43)

and $\lambda$ is determined by imposing that $x_{1}^{2}+x_{2}^{2}=L^{2}$. The
solution of this system provides not only the dynamics of the mass, but also
the tension $\lambda$ of the rigid spring. It is possible, with suitable
manipulations, the classical equations of motion of the pendulum

$L\frac{d^{2}\theta}{dt^{2}}+g\sin\theta=0$ (1.44)

starting from equations (1.43) and (1.40) (see Exercise 1).

Fig. 1.8 shows the trajectory of the mass of 1kg, oscillating in a pendulum
with rest length $L=$ 1m, starting from an angle $\theta_{0}=30^{\circ}$, with
an initial length $L$, for various values of the stiffness parameter. As the
stiffness $k$ of the spring increases, the trajectories resemble more and more
those of a classical pendulum.

Figure 1.8: In these tree plots we see different oscillations for different
elastic constants, $k$, namely $k=10^{2},10^{3},10^{3}$ (N/m). When
$k\to\infty$, the oscillator behaves like a classical pendulum.

###  1.3 Coupling conditions for sorption kinetics at the oscillating cell
membrane

Figure 1.9: Scheme of surfactant disposition. They are composed by hydrophilic
heads, that prefer to stay in the water, and and hydrophobic tails.

The term surfactant comes from surface active agent. They are composed by
anions and cations, that have two different configurations. The negative ions
are amphiphilic molecules, i.e., they both have hydrophilic and hydrophobic
parts, and, for this reason, they are absorbed at the air-water interface. The
hydrophobic part prefers to stay in the air, while the hydrophilic part stays
in water. This will cause a decrease in the surface (or interfacial) tensions.

In this section we see how to treat different orders of magnitude in space, in
the same mathematical model. How we were saying in the Introduction, an
expansion in $\varepsilon$ (the range of the attractive-repulsive potential)
will be considered, in order to obtain a new boundary condition for the
reduced multiscale model.

Introducing the local concentration of ions $c=c(\vec{x},t)$, its time
evolution in a fluid is governed by following the conservation law

$\displaystyle\displaystyle\frac{\partial c\left(\vec{x},t\right)}{\partial
t}$ $\displaystyle=-\nabla\cdot J(\vec{x},t)\qquad{\rm in}\,\Omega$ (1.45)
$\displaystyle\displaystyle J$ $\displaystyle=\ -D\left(\nabla c+\
\frac{1}{k_{B}T}c\nabla V\right)\qquad{\rm in}\,\Omega$ (1.46)

where $k_{B}$ is the Boltzmann’s constant, $T$ is the absolute temperature
(assumed to be constant) and $V(x)$ is a suitable potential function that
simulates the attractive-repulsive behavior of the bubble with the particles.

A part from the attraction, the potential models also the impermeability of
the bubble, thus, at very short distances, it is designed to repulse the
particles from the surface of the trap, and, for this reason, its derivative
becomes negative in the proximity of the bubble radius.

The multiscale model has been derived and tested in various dimensions $d$.

####  1.3.1 1D Multiscale model

Figure 1.10: Representation of the external potentials $V(x)$ (in Eq. (1.49)),
on the left, and $U(\xi)$ (in Eq. (1.50)), on the right, for
$\varepsilon=0.05$. $U(\xi)$ is obtained from a change of variable of $V(x)$,
with $\xi\in[0,L+1]$.

This subsection is devoted to the derivation of the multiscale model in 1D
space dimension.

The drift-diffusion Eqs. (1.45-1.46) in 1D read

$\displaystyle\displaystyle\frac{\partial c}{\partial t}+\frac{\partial
J}{\partial x}=0\qquad{\rm in}\,\Omega$ (1.47)

where

$\displaystyle J=-D\left(\frac{\partial c}{\partial
x}+\frac{1}{k_{B}T}cV^{\prime}\right),$ (1.48)

and the full domain is $\Omega=[-\varepsilon,1]$, with
$\Omega^{\varepsilon}_{b}=[-\varepsilon,L\varepsilon]$ the bubble domain,
inside which the ions feel the effect of the potential, and the fluid domain
is $\Omega^{\varepsilon}_{f}=[L\varepsilon,1],$ which is not affected by the
bubble, see Fig. 1.10. Therefore, the potential has effect only in
$\Omega^{\varepsilon}_{b}$, while it is constant for
$x\in\Omega^{\varepsilon}_{f}$, say $V(x)=0$.

A suitable potential is the Lennard-Jones (LJ) type, which describes
attraction at long distances and repulsion at short distances, due to Van der
Waals and Pauli terms, respectively. A typical shape of the potential $V(x)$
is shown in Fig. 1.10 and it takes the form

$\displaystyle
V(x)=E\left(\left(\frac{x+\varepsilon}{\varepsilon}\right)^{-12}-2\left(\frac{x+\varepsilon}{\varepsilon}\right)^{-6}\right)$
(1.49)

where $\varepsilon$ denotes the range of the potential and $E$ represents the
depth of the well (see Fig. 1.10, left panel). It is convenient to adopt a non
dimensional form of the potential, expressed as a function of the rescaled
variable $\xi=1+x/\varepsilon\in[0,L+1]$:

$U(\xi)=\phi\left(\xi^{-12}-2\xi^{-6}\right)$ (1.50)

where $\displaystyle\phi={E}/{k_{B}T}$ represents the ratio between the depth,
$E$, of the potential well of the bubble and $k_{B}T$. $L$ is a scaled
distance beyond which the potential is negligible.

The assumption made is that the trap is located at $x=0$ and its ”thickness”
is of order of $\varepsilon$. The original model requires to solve Eq. (1.47)
in $\Omega^{\varepsilon}=\Omega^{\varepsilon}_{b}\cup\Omega^{\varepsilon}_{f}$
with initial condition $c=c_{0}(x)$ when $t=0$ and with boundary conditions

$J=0\quad{\rm at}\>x=-\varepsilon,{\rm and}\>x=1$ (1.51)

to simulate wall effects.

It is clear that the solution $c$ and the flux $J$ depend on $\varepsilon$ as
well, and then the original problem can be written as

$\displaystyle\displaystyle\frac{\partial c_{\varepsilon}}{\partial
t}+\frac{\partial J_{\varepsilon}}{\partial x}$ $\displaystyle=$
$\displaystyle 0\quad\rm{in}\,\Omega^{\varepsilon}$ (1.52)
$\displaystyle\displaystyle J_{\varepsilon}$ $\displaystyle=$
$\displaystyle-D\left(\frac{\partial c_{\varepsilon}}{\partial
x}+\frac{1}{k_{B}T}\,c_{\varepsilon}\,V^{\prime}_{\varepsilon}\right).$ (1.53)

It is important to notice that the space derivative can be expressed ad
$\partial/\partial x=\partial/\partial\xi\cdot\partial\xi/\partial
x=\partial/\partial\xi\cdot 1/\varepsilon$, and, with this notation, Eq.
(1.53) becomes

$J_{\varepsilon}=-D\frac{1}{\varepsilon}\left(\frac{\partial
c_{\varepsilon}}{\partial\xi}+c_{\varepsilon}\,U^{\prime}\right).$ (1.54)

Rewriting 1.52 in the rescaled variable and using Eq. (1.54), it becomes

$\frac{\partial c_{\varepsilon}}{\partial
t}=-\frac{1}{\varepsilon}\frac{\partial
J_{\varepsilon}}{\partial\xi}=\frac{D}{\varepsilon^{2}}\frac{\partial}{\partial\xi}\left(\frac{\partial
c_{\varepsilon}}{\partial\xi}+c_{\varepsilon}U^{\prime}(\xi)\right).$ (1.55)

The range of the scaled variable $\xi$ does not depend on $\varepsilon$, then
the following expansion in $\Omega^{\varepsilon}_{b}$ is valid for the
concentration $c_{\varepsilon}(\xi,t)$ :

${c_{\varepsilon}(\xi,t)=c^{(0)}(\xi,t)+\varepsilon^{2}c^{(1)}(\xi,t)+\varepsilon^{4}c^{(2)}(\xi,t)+\cdots.}$
(1.56)

Inserting the expansion 1.56 into Eq. (1.54), the following expansion for the
flux hold true

$\varepsilon
J_{\varepsilon}=J^{(0)}+\varepsilon^{2}J^{(1)}+\varepsilon^{4}J^{(2)}+\cdots,$
(1.57)

with

$J^{(k)}=-D\left(\frac{\partial
c^{(k)}}{\partial\xi}+c^{(k)}U^{\prime}(\xi)\right),\quad k\geq 0.$

Using the expansion 1.57 in 1.55, the Eqs. (1.52-1.53) become, to the various
orders in $\varepsilon$,

$\displaystyle O(\varepsilon^{-1}):$
$\displaystyle\displaystyle\phantom{\quad\frac{\partial c^{(k)}}{\partial
t}+}\frac{\partial J^{(0)}}{\partial\xi}$ $\displaystyle=0,$ (1.58)
$\displaystyle O(\varepsilon^{k}):$
$\displaystyle\displaystyle\quad\frac{\partial c^{(k)}}{\partial
t}+\frac{\partial J^{(k)}}{\partial\xi}$ $\displaystyle=0,\quad k\geq 0.$
(1.59)

where Eq. (1.58) states that the lowest order flux $J^{(0)}$ is constant.
Condition $\lim_{x\to-\varepsilon}J(x,t)=0$, which is a consequence of the
singularity of the potential, gives:

$\frac{\partial c^{(0)}}{\partial\xi}+U^{\prime}(\xi)c^{(0)}=0.$ (1.60)

Therefore, it is possible to integrate the equation

$\frac{1}{c^{(0)}}\frac{\partial c^{(0)}}{\partial\xi}=-U^{\prime}(\xi)$
(1.61)

whose solution is

$c^{(0)}(\xi,t)=c^{(0)}(L+1,t)\exp\left(-U(\xi)+U(L+1)\right)=c^{(0)}(L+1,t)\exp\left(-U(\xi)\right)$
(1.62)

where $U(L+1)=0$ because the potential has effect only inside the bubble
domain (that is $\Omega^{\varepsilon}_{b}$ for $V(x)$ and $[0,L+1]$ for
$U(\xi)$), while it is constant outside, say $U(\xi)=0$ for $\xi\geq L+1$.

The following assumption is made

$V_{x}(\varepsilon L)=0,$

which is compatible with having a smooth potential with compact support.
Notice that, strictly speaking, LJ potential is not with compact support,
however it is negligible small, together with its derivative, at distances
which are sufficiently larger than its range $\varepsilon$.

Integrating 1.47 in the range of the potential, using the approximation
$c(x,t)=c^{(0)}(\xi,t),\,x\in\varepsilon[-1,L]$, the zero boundary condition
1.51 at $x=-\varepsilon$ and that $V_{\varepsilon}(\varepsilon
L)=V^{\prime}_{\varepsilon}(\varepsilon L)=0$, we obtain

$\frac{d}{dt}\int_{-\varepsilon}^{\varepsilon
L}c(x,t)\,dx+J_{\varepsilon}(\varepsilon L)=0,$

for which it follows

$\displaystyle\frac{\partial c(\varepsilon L,t)}{\partial
t}\;\varepsilon\int_{0}^{L+1}\exp\left(-U(\xi)\right)d\xi-D\frac{\partial
c(\varepsilon L,t)}{\partial x}$ $\displaystyle=$ $\displaystyle 0$

that represents a boundary condition at $x=\varepsilon L$. Notice that, out of
the range of the potential, the flux $J$ is given by just the diffusion term,
thus Eq. (1.59) will not be used since for sufficiently small values of
$\varepsilon$ it is $c(x,t)\approx c^{(0)}(x,t)$, and Eq. (1.58) is used to
provide an effective boundary condition for the diffusion equation.

Summarizing, the multiscale model can be stated as a limit model for
$\varepsilon\to 0$ as follows:

$\displaystyle\frac{\partial c}{\partial t}$ $\displaystyle=$ $\displaystyle
D\frac{\partial^{2}c}{\partial x^{2}}\quad{\rm{in}}\,x\in[0,1]$ (1.63)
$\displaystyle\frac{\partial c}{\partial x}$ $\displaystyle=$ $\displaystyle
0\quad{\rm{at}}\,x=1$ (1.64) $\displaystyle M\frac{\partial c}{\partial t}$
$\displaystyle=$ $\displaystyle D\frac{\partial c}{\partial
x}\quad{\rm{at}}\,x=0$ (1.65)

where

$M=\varepsilon\int_{0}^{L+1}\exp\left(-U(\xi)\right)d\xi.$ (1.66)

As $\varepsilon\to 0$ one could let $L\to\infty$ in such a way that
$L\varepsilon\to 0$, so in the expression for $M$ the integral can be written
from $0$ to $\infty$. In practice, $L$ is finite, and it is chosen $L=2$.

####  1.3.2  Extension to more dimensions

Figure 1.11: Scheme of the boundary condition in 2D where $\hat{n}$ is the
outgoing normal vector to the domain $\Omega$, and
$\varepsilon=\varepsilon(S)$ is the thickness parameter of the bubble surface
$S$. The internal line near the bubble surface represents the boundary
$\Gamma_{B}$.

In this section the condition in Eq. (1.65) is extended to 2D and 3D. The
bubble $\mathcal{B}$ is defined by a level set function (see e.g. [27, 28, 29,
30])

$\phi(\vec{x})=r-\sqrt{x^{2}+y^{2}},$ (1.67)

such that $\phi(\vec{x})<0$ inside the bubble and $\phi(\vec{x})>0$ outside.
The unit normal $\hat{n}$ can be defined as

$\hat{n}(\vec{x})=\frac{\nabla\phi(\vec{x})}{|\nabla\phi(\vec{x})|}\qquad\forall\vec{x}\in\Gamma_{B}$
(1.68)

where $\Gamma_{B}=\partial\mathcal{B}$.

Near the boundary, the flux can be decomposed into the normal and the
tangential components, such that

$\vec{J}=J_{n}\hat{n}+J_{\tau}\hat{\tau}$ (1.69)

where $J_{n}$ and $J_{\tau}$ denote the components of $\vec{J}$ parallel and
orthogonal to $\nabla V$, respectively. The equation for the concentration $c$
becomes:

$\frac{\partial c}{\partial t}=-\left(\frac{\partial
J_{n}}{\partial\hat{n}}+\frac{\partial J_{\tau}}{\partial\hat{\tau}}\right).$
(1.70)

Let $\vec{x}$ be a point in $\Gamma_{B}$ and, integrating along the normal
direction, it becomes

$\displaystyle\displaystyle\frac{\partial}{\partial
t}\int_{-\varepsilon}^{\varepsilon L}c(\vec{x}+r\hat{n})dr$
$\displaystyle=-J_{n}(\vec{x}-L\varepsilon\hat{n})-\int_{-\varepsilon}^{\varepsilon
L}\frac{\partial J_{\tau}}{\partial\tau}(\vec{x}+r\hat{n})dr$ (1.71)

since $J_{n}(\vec{x}+\varepsilon\hat{n})=0$. Now, we consider a change of
variable from $\vec{x}$ to $\xi$ in the integration intervals and, making use
of $\displaystyle J_{\tau}=-D\frac{\partial c}{\partial\tau}$, the Eq. (1.71)
becomes

$\varepsilon\frac{\partial}{\partial
t}\int_{0}^{L+1}c(\vec{x}+\varepsilon\xi\hat{n})d\xi=-J_{n}(\vec{x})+\varepsilon\int_{0}^{L+1}D\frac{\partial^{2}c}{\partial\tau^{2}}(\vec{x}+\varepsilon\xi\hat{n})d\xi.$
(1.72)

Let us consider the term $J_{n}$:

$J_{n}=-D\left(\frac{\partial
c}{\partial\hat{n}}+\frac{1}{k_{B}T}\frac{\partial
V}{\partial\hat{n}}c\right),$ (1.73)

since $\displaystyle V(x)/k_{B}T=U\left(\xi\right)$, it follows that

$\displaystyle\frac{\partial(V/k_{B}T)}{\partial\hat{n}}=\frac{1}{\varepsilon}U^{\prime}(\xi)$

with $\xi=(r-R+\varepsilon)/\varepsilon$. For this reason, there is a
dependence on $\xi$ of the solution inside the layer,
$c=c(\vec{x}+\varepsilon\xi\hat{n},t)=\tilde{c}(\xi,t)$, thus there is a
common factor $1/\varepsilon$:

$J_{n}=-D\frac{1}{\varepsilon}\left(\frac{\partial
c}{\partial\xi}+\frac{\partial U}{\partial\xi}c\right)$ (1.74)

where we omitted the $\tilde{\cdot}$. Following the same argument used in 1D
to derive Eq. (1.56), the expansion in $\varepsilon$ of the solution, to the
lowest order is

$\frac{\partial c^{(0)}}{\partial\xi}+\frac{\partial
U}{\partial\xi}c^{(0)}=0.$ (1.75)

Integrating in $d\xi$, the solution of the Eq.(1.75) is

$c^{(0)}(\vec{x}+\varepsilon\xi\hat{n})=c^{(0)}(\vec{x})\exp\left(-U(\xi)\right).$
(1.76)

Using this expression we can compute the line density of entrapped ions:

$\mathcal{C}(\vec{x})\varepsilon\int_{0}^{L+1}c(\vec{x}+\varepsilon\xi\hat{n})d\xi\simeq
Mc(\vec{x})$ (1.77)

with $\displaystyle M=\varepsilon\int_{0}^{1}\exp\left(-U(\xi)\right)d\xi$ and
where we assumed $U(L+1)=U^{\prime}(L+1)=0$.

At the end, in 2D (and also in 3D), the system and the new boundary condition
for the concentration becomes:

$\displaystyle\displaystyle\frac{\partial c}{\partial t}$ $\displaystyle=$
$\displaystyle D\Delta c\quad\rm{in}\>\Omega$ (1.78)
$\displaystyle\displaystyle\frac{\partial c}{\partial n}$ $\displaystyle=$
$\displaystyle 0\quad\rm{on}\,\partial\Omega\backslash\Sigma_{B}$
$\displaystyle\displaystyle M\frac{\partial c}{\partial t}$ $\displaystyle=$
$\displaystyle MD\Delta_{\perp}c+D\frac{\partial c}{\partial
n}\quad\rm{on}\,\Sigma_{B}$ (1.79)

because $J_{n}=-D\partial c/\partial n$ just out of the range of the
potential.

Notice that in 2D the Laplace-Beltrami operator reduces to the second
derivative with respect to the arclength of the boundary:333Here by
$\partial/\partial\tau$ we denote the derivative on $\Gamma_{B}$, i.e. the
derivative along the arclength that parametrizes the curve, likewise
$\partial^{2}/\partial\tau^{2}$ denotes the second derivative along
$\Gamma_{B}$, not the second derivative along the tangent direction. See Fig.
1.11 and Appendix 1.6.2.

$\Delta_{\perp}c=\frac{\partial^{2}c}{\partial\tau^{2}}.$ (1.80)

####  1.3.3 Validation of the model and results

In this section we validate the model in 1D and 2D.

First, a 1D domain is considered as intermediate step. Then, a more realistic
2D domain is considered, where a circular hole, at the center of a squared
domain, describes the attractive bubble surface.

##### Validation 1D

The initial domain in 1D is
$\Omega^{\varepsilon}=[-\varepsilon,1]=\Omega^{\varepsilon}_{b}\cup\Omega^{\varepsilon}_{f},\,\Omega^{\varepsilon}_{b}=[-\varepsilon,2\varepsilon],\,\Omega^{\varepsilon}_{f}=[2\varepsilon,1]$,
and the space derivatives are discretized with Finite-Difference (FD) schemes,
see Chapter ?? for more details. In this case, we refer to the system
(1.47-1.49).

For the multiscale model, the domain is $\Omega^{0}=[0,1]$, and the equations
are in (1.63-1.66). Internal equations are discretized by finite differences
by the classical three points approximation of the second derivative, while
the Eq. (1.65) is discretized to second order as

$M\frac{\partial}{\partial t}\frac{c_{0}+c_{1}}{2}=D\frac{c_{1}-c_{0}}{h}.$
(1.81)

Let $c^{\varepsilon}$ and $c^{0}$ be the solutions of the full and multiscale
models, respectively, and their comparison is showed in Fig. 1.12.

Figure 1.12: Here there are two aspects of the same plot. On the left it is
shown the effect of the attractive potential that becomes higher decreasing
the value of $\varepsilon$. On the right there is the comparison between
$c^{\varepsilon}$, and the dashed red line that represents the solution of the
multiscale model, $c^{0}$.

##### Validation 2D in presence of a circular bubble

The arbitrary boundary of the bubble is described by a level set function. The
set of equations is then solved by a finite difference scheme on a regular
Cartesian mesh, which makes use of regular grid points (inside the domain
$\Omega^{0}$) and ghost points (neighbors of regular grid points which are out
of $\Omega^{0}$). The equations on the regular points are obtained by finite
difference discretization of the Laplacian operator, while the equations on
each ghost point are given by interpolating the numerical solution near the
ghost point and imposing the boundary condition on the interpolant at the
orthogonal projection of the point on $\partial\Omega^{0}$ (see [31, 32] for
more details).

Figure 1.13: The effect of the attractive potential in a 2D domain. Here the
bubble has a reasonable circular shape. As before, in the first three panels
(a,b,c) there are the solution $c^{\varepsilon}$ for three different values of
$\varepsilon$. In panel (d) there is a 1D plot, that is a section of the
previous ones for $y=0$ with an additional line. The black dashed line
represents the solution of the multiscale model, $c^{0}$.

###  1.4 Coupling conditions for cardiovascular system

Talking about the circulatory system, it is very important to introduce the
concept of geometric multiscale. In this case, it is related to the coupling
of different dimensions and scales, in the same system of equations.
Considering the relationships between different scales and dimensions, it is
possible to gain a more comprehensive understanding of the system.

The cardiovascular system is very complex for different reasons. Apart from
the geometric complexity, it is very important to mention the difficulty in
obtaining data that are used to deduce meaningful initial and boundary
conditions, when an artificial ’chopped’ domain is considered for
computational purposes (see, e.g. Fig. 1.14). Local artificial domains may
have global consequences, and to take into consideration these properties,
suitable boundary conditions need to be applied between the free surfaces of
the artificial domain.

\begin{overpic}[abs,width=260.17464pt,unit=1mm,scale={.25}]{Figures_new/chopped_domains.pdf}
\put(10.0,22.0){{$\Sigma^{t}$}} \put(18.0,65.0){{$\Omega^{t}_{f}$}}
\put(48.0,65.0){{$\Omega_{s}$}} \put(40.0,55.0){{$\Sigma$}} \end{overpic}
Figure 1.14: An example of a chopped domain. On the left, the purple object
represents the time dependent fluid domain, while the red one represents the
structure domain.

####  1.4.1 The 3D model

Figure 1.15: Classification of Newtonian and Non-Newtonian fluids

In order to obtain a realistic description of a vessel flow, it is convenient
to refer to a 3D Navier-Stokes formulation for incompressible homogeneous
Newtonian fluids. Since the diameter of blood cells is about $10^{-5}m$, while
the smallest diameter of arteries and veins is of the order of $10^{-3}m$,
linear equations govern the relationships between internal forces and velocity
gradient, and blood can be treated as a called Newtonian fluid (see, e.g. [21,
22]). In Fig. 1.15 we see few examples of Newtonian (water, Helium and Oxygen)
and Non-Newtonian fluids (tooth paste, oil, blood). As it is shown in that
figure, the blood belongs to the Non-Newtonian fluids, but with the
approximations mentioned above, it is possible to consider it as a Newtonian
fluid.

An Eulerian formulation is more convenient for the fluid equations and the
frame of reference is denoted by $\Omega_{\rm f}\subset\mathbb{R}^{3}$ (see
Fig. 1.14 (left panel)). While, for deformable solid structure (the arteries)
a Lagrangian framework is more suitable (see Fig. 1.14 (right panel)) and the
reference domain is $\tilde{\Omega}_{\rm s}\subset\mathbb{R}^{3}$. Therefore,
there exists a map $\mathcal{L}:\tilde{\Omega}_{\rm s}\to{\Omega}_{\rm s}$,
where $\Omega_{\rm s}$ is the spatial domain. $\Sigma$ will be the interface
between the fluid and the structure domains.

Now it is possible to define the system for the velocity field $\vec{u}$, the
fluid pressure $p$ and the density $\rho$, both for the fluid and the
structure parts, such that,

$\partial_{t}{\bf V}_{\rm l}+M_{\rm l}({\bf V}_{\rm l})\,\nabla\cdot{\bf
V}_{\rm l}=0,\quad M_{\rm l}({\bf V}_{\rm l})\in\mathbb{R}^{3,3},\quad{\rm
l=f,s}$ (1.82)

where $l=f$ stands for the fluid system and variables and $l=s$ for the
structure ones. The vectors of unknowns are

${\bf V}_{\rm l}^{\text{T}}=(\rho_{\rm l},u_{\rm l},p_{\rm l}),\quad{\rm
l=f,s}.$

If the matrix $M_{\rm l},\,{\rm l=f,s}$ admits real eigenvalues, and it is
diagonalizable (i.e., the eigenvectors form a bases) the system is called
hyperbolic, and it is well-posed.

###### Example 1

Let us consider Euler’s system, in the case of ideal gases. The equation for
the conservation of mass is the following

$\frac{\partial\rho}{\partial
t}+u_{i}\partial_{i}\rho+\rho\partial_{i}u_{i}=0\quad\Rightarrow\quad\frac{D\rho}{Dt}+\rho\nabla\cdot\vec{u}=0,$

the one for the conservation of the momentum reads

$\rho\frac{Du}{Dt}+\nabla p=0,$

while the relation for the energy can be reduced as

$\frac{Dp}{Dt}=\biggl{(}\frac{\partial
p}{\partial\rho}\biggr{)}_{S}\frac{D\rho}{Dt}=a^{2}(-\rho\nabla\cdot\vec{u}),\quad
a^{2}:=\biggl{(}\frac{\partial p}{\partial\rho}\biggr{)}_{S}.$ (1.83)

In this way we obtain a system for the variables $\rho$, $p$ and $\vec{u}$:

$\left\\{\begin{aligned} &\frac{D\rho}{Dt}+\rho\nabla\vec{u}=0\\\
&\frac{Du}{Dt}+\frac{1}{\rho}\nabla p=0\\\
&\frac{Dp}{Dt}+a^{2}\rho\nabla\cdot\vec{u}=0\\\ \end{aligned}\right.$ (1.84)

For simplicity we consider the system (1.84), in 1D, that becomes

$\left\\{\begin{aligned}
\partial_{t}\rho&+u\,\partial_{x}\rho+\rho\,\partial_{x}u&=0\\\
\partial_{t}u&+u\,\partial_{x}u+\frac{1}{\rho}\,\partial_{x}p&=0\\\
\partial_{t}p&+u\,\partial_{x}p+a^{2}\rho\,\partial_{x}u&=0\\\
\end{aligned}\right.$ (1.85)

Now we define the following vector

$V^{\text{T}}=(\rho,u,p)$

and the system (1.85) can be rewritten as

$\partial_{t}V+M(V)\,\partial_{x}V=0,\quad M(V)=\begin{pmatrix}u&\rho&0\\\
0&u&1/\rho\\\ 0&a^{2}\rho&u\end{pmatrix}$ (1.86)

Now we calculate the three eigenvalues of the matrix $M$:

$\text{det}(\lambda\mathbb{I}-M)=(\lambda-u)\bigl{[}(\lambda-u)^{2}-a^{2}\bigr{]}=0$
(1.87)

that are $\lambda_{0}=u$ e $\lambda_{\pm}=u\pm a$, and they belong to
$\mathbb{R}$ if $a\in\mathbb{R}$.

If $a\in\mathbb{C}$ the system is not hyperbolic!

####  1.4.2 The 1D and 0D model

Modern supercomputers struggle with the complexity of creating a complete 3D
model of the cardiovascular system. The computational process for such a model
takes hours or even days, whereas 1D or 0D models achieve good accuracy in a
much shorter time frame. Interestingly, in certain situations, the accuracy
provided by the 3D models exceeds the requirements of the problem at hand.

The focus should be on obtaining a detailed depiction of blood flow in
specific critical areas where a pathology is identified or suspected. For this
purpose, a highly efficient approach involves utilizing a comprehensive 3D
model solely in the designated region, while employing simplified 1D or 0D
models elsewhere [9, 11].The 1D-branches are then connected to the rest of the
circulatory system by 0D equations, i.e. by a system governed by ODE’s.

At the end, the entire domain will be a combination of several domains of 3D,
1D and 0D models. What is important is the continuity relations at the
junctions between them. The remaining domain is represented by attaching a
structured tree, to simulate the spreading of certain tissues, see Fig. 1.16.

Figure 1.16: Systemic arterial tree consisting of larger arteries originates
at the heart (A) and terminates at the points marked with (B) and structured
trees representing smaller arteries originate at these terminals and provide
main tree with outflow boundary conditions.

To deduce a 1D model, some assumptions are needed:

  1. 1.

the axis $z$ of the cylinder is assumed rectilinear and the component $u_{z}$
of the velocity is the dominant one

  2. 2.

all the cross sections of the cylinders are circles $\mathcal{S}(t,z)$ with
radius $R(t,z)$, and the pressure is considered constant over each section
$\mathcal{S}(t,z)$

  3. 3.

the membranes have constant thickness

To derive the 1D model, it is enough to integrate the Eqs. (1.82) over the
section $\mathcal{S}$ in the $z$-direction. An additional order reduction is
obtained integrating the system along the axis direction, obtaining the 0D
model. In this kind of models, only the time dependence survives, and they are
common to describe phenomena in capillaries.

####  1.4.3 Coupling of 3D-1D, 3D-0D and 1D-0D models

To accurately describe the entire cardiovascular system, it is crucial to
select the appropriate spatial dimensions that adequately represent each
component, depending on the desired level of accuracy. In certain regions,
where precision is of utmost importance, a 3D model is employed. However, in
other areas, computational complexity can be reduced by utilizing simplified
models such as 1D or 0D.

When combining both 3D and 1D (or 0D) models, a combination of these models is
utilized in specific regions, requiring an approximation to ensure continuity
of the quantities involved. To facilitate this process, the first step
involves a technique called domain splitting. It implies the division of a
cylindrical vessel into two subdomains, as depicted in the accompanying Fig.
1.17. Subsequently, the 3D-1D (or 3D-0D) model can be formulated by
incorporating appropriate interface conditions.

\begin{overpic}[abs,width=151.76964pt,unit=1mm,scale={.25}]{Figures_new/domain_splitting.pdf}
\put(-7.0,56.5){(1)} \put(22.0,51.5){$\Sigma$} \put(17.0,43.0){$\Omega$}
\put(-7.0,38.0){(2)} \put(12.0,27.0){$\Omega_{\rm f}$}
\put(28.0,23.0){$\Omega_{\rm s}$} \put(-7.0,18.0){(3)}
\put(23.0,16.0){$\Gamma=\\{\Gamma_{\rm f},\Gamma_{\rm s}\\}$}
\put(25.0,11.0){{$z=0$}} \put(42.0,6.5){{$z=L$}} \end{overpic} Figure 1.17:
Scheme of domain splitting (the subscript ${\rm f}$ stands for fluid and ${\rm
s}$ for solid): (1) $\Omega$ is the initial considered domain and $\Sigma$ the
related boundary. (2) $\Omega_{\rm f}\cup\Omega_{\rm s}$ is a partition of
$\Omega$ and $\Gamma$ is the common interface. (3) The 3D domain $\Omega_{\rm
s}$ is replaced with a 1D domain $[0,L]$ in $z$-direction.

Regarding the 3D-0D coupling, there are three properties that mathematical
modeling would face: (i) the non linear terms are neglected; (ii) the boundary
conditions will be integrated in the equations; (iii) domain splitting (see
Fig. 1.17) is required, and some continuity property can be lost (the
subscript ${\rm f}$ stands for fluid and ${\rm s}$ for solid).

At the end, we show the strategy of the 1D-0D coupling. Suppose there is a 1D
model

$\mathcal{M}_{1D}({\bf V}_{1D},\Omega)=0$ (1.88)

where ${\bf V}_{1D}$ is the vector variables and $\Omega=[-L,L]$ is the
domain. The problem can be split in two sub-problems $\mathcal{M}^{1}_{1D}$
and $\mathcal{M}^{2}_{1D}$ and the domain in two subdomains
$\Omega_{1}=[-L,0]$ and $\Omega_{2}=[0,L]$:

$\displaystyle\mathcal{M}^{1}_{1D}({\bf V}^{1}_{1D},[-L,0])=$ $\displaystyle
0,\quad{\rm in}\,\Omega_{1}$ (1.89) $\displaystyle\mathcal{M}^{2}_{1D}({\bf
V}^{2}_{1D},[0,L])=$ $\displaystyle 0,\quad{\rm in}\,\Omega_{2}$ (1.90)
$\displaystyle{\bf V}^{1}_{1D}(z=0)=$ $\displaystyle{\bf V}^{2}_{1D}(z=0)$
(1.91)

where Eqs. (1.89),(1.90) are the equations for the domains $\Omega_{1}$ and
$\Omega_{2}$, respectively, while Eqs. (1.91) is the continuity conditions at
the intersection between $\Omega_{1}$ and $\Omega_{2}$, $z=0$. The next step
is to replace the 1D problem in $\Omega_{2}$, with a suitable 0D model, with
appropriate boundary conditions. To do this, average quantities are needed in
0D, leading to the following problem

$\displaystyle\mathcal{M}^{1}_{1D}({\bf V}^{1}_{1D},[-L,0])=$ $\displaystyle
0,\quad{\rm in}\,\Omega_{1}$ (1.92) $\displaystyle{\bf V}^{1}_{1D}(z=0)=$
$\displaystyle{\bf V}_{0D}$ (1.93) $\displaystyle\mathcal{M}_{0D}({\bf
V}_{0D})=$ $\displaystyle 0.$ (1.94)

###  1.5 Coupling conditions for shallow water equations at channel junctions

The shallow water equations is a hyperbolic system that describes surface
water flows, such as in artificial channels, rivers, lakes and oceans. The
system consists of a mass and a momentum balance equation, and, in the most
common scenario, it is considered in a two dimensional domain.

This section is devoted to the description of an efficient way to consider one
dimensional problems. In Fig. 1.6 we see the schematic representation of the
intersection of three channels. As a reference example, we consider the first
channel the one that is parallel to the $x$-axis, and the other two create two
angles with the first one. The figure shows how to represent with a 1D domain
a 2D problem.

As we said in the Introduction, the tricky part is to define the boundary
condition of each channel at the junction, and to couple them with the SW
equations of the three channels. The main idea is to see the junctions as in
the Fig. 1.6, focusing at the intersections with a more detailed 2D domain.

Choosing the coupling conditions it can be seen as a Riemann problem defined
by imposing conservation of mass and continuity of the water height through
junctions. To explain better the concept, we can say that the following
problem is solved at the junctions:

$\left\\{\begin{array}[]{c}\displaystyle\partial_{t}u_{l}+\partial_{x}f(u_{l})=0,\quad
t\in\mathbb{R}^{+},\quad l=1,2,3\\\ u_{l}(0,x)=\bar{u}_{l},\quad
x\in\mathbb{R}^{+},\,u_{l}\in\Omega,\quad l=1,2,3.\end{array}\right.$ (1.95)

where $\bar{u}_{l},l=1,2,3$ are the constant states in $\Omega$. $u_{l}$ has
two components, i.e., $u_{l}=(u_{l,L},u_{l,R}),l=1,2,3,$ that denote the
densities of the conserved quantities in the $l-th$ tube.

####  1.5.1 Mathematical model

As we said before, the main strategy is to consider 1D models along the
straight channels, and 2D models at the junctions, in order to obtain detailed
boundary conditions.

The shallow water equations in 2D read:

$\partial_{t}{\bf Q}_{2D}+\partial_{x}{\bf F}_{2D}({\bf
Q}_{2D})+\partial_{y}{\bf G}_{2D}({\bf Q}_{2D})={\bf S}({\bf Q}_{2D})$ (1.96)

where ${\bf Q}_{2D}=[h,hu,hv]^{T}$ is the vector of conserved variables, ${\bf
F}_{2D}({\bf Q}_{2D})=[hu,hu^{2}+\frac{1}{2}gh^{2},huv]^{T}$ is the flux in
$x$-direction and ${\bf G}_{2D}({\bf
Q}_{2D})=[hv,hvu,hv^{2}+\frac{1}{2}gh^{2}]^{T}$ is the flux in $y$-direction,
$g$ is the gravity and ${\bf S}_{2D}({\bf Q}_{2D})=[0,-ghS_{x},-ghS_{y}]^{T}$
is the variation of the bottom topography, with $S_{x},S_{y}$ known functions.

The same model, in 1D, reads

$\partial_{t}{\bf Q}_{1D}+\partial_{s}{\bf F}_{1D}({\bf Q}_{1D})={\bf
S}_{1D}({\bf Q}_{1D})$ (1.97)

where the second component of the velocity $v$ is dropped, the $y$-direction
the flux G disappears and $u_{1D}$ is the velocity along the $s$-direction.

It is possible to combine the 1D system, in Eq. (1.97), defined in a generic
direction $s$, with the 2D system defined in Eq. (1.96). We start integrating
Eq. (1.96), that becomes

$\displaystyle\frac{\partial}{\partial t}\int_{V}{\bf
Q}_{2D}dV+\int_{\Omega}\left(\cos(\theta){\bf F}_{2D}({\bf
Q}_{2D})+\sin(\theta){\bf G}_{2D}({\bf Q}_{2D})\right)d\Omega={\bf 0}.$ (1.98)

where $[\cos\theta,\sin\theta]$ are the two components of the outward unit
normal vector, such that, ${\bf n}=[n_{1},n_{2}]=[\cos\theta,\sin\theta]$, and
$V$ is the control volume in the Cartesian plane.

If ${\bf T}={\bf T}(\theta)$ is a rotation matrix, and ${\bf T}^{-1}(\theta)$
its inverse, such that

$T=\left[\begin{array}[]{ccc}1&0&0\\\ 0&\cos(\theta)&\sin(\theta)\\\
0&-\sin(\theta)&\cos(\theta)\\\ \end{array}\right],\qquad
T^{-1}=\left[\begin{array}[]{ccc}1&0&0\\\ 0&\cos(\theta)&-\sin(\theta)\\\
0&\sin(\theta)&\cos(\theta)\\\ \end{array}\right]$ (1.99)

it is easy to prove that

${\bf H}={\bf n}\cdot[{\bf F}({\bf Q}),{\bf G}({\bf Q})]=\cos(\theta){\bf
F}({\bf Q})+\sin(\theta){\bf G}({\bf Q})={\bf T}^{-1}{\bf F}({\bf T}({\bf
Q}))$ (1.100)

for all vectors ${\bf Q}$ and all normal directions of $\Omega$

The final steps to describe are the following:

  * •

1D system (Eq. (1.97)) where the direction $s$ is the normal direction, with
rotated Riemann initial condition

  * •

calculate the angle $\theta_{e}$ between the outward unit normal to the edge
$e$ and the reference $x$ direction, being positive in the counter-clock
direction

  * •

calculate the corresponding rotation matrix ${\bf T_{e}}$ and rotate left and
right initial data

  * •

solve the 1D Riemann problem and compute an approximation of the flux in the
edge $e$

  * •

rotate back the flux (see Eq. (1.100))

\begin{overpic}[abs,width=238.49231pt,unit=1mm,scale={.25}]{Figures_new/junctions2_1D_2D.pdf}
\put(18.5,30.0){1D} \put(52.0,45.0){1D} \put(57.0,16.0){1D}
\put(38.0,26.0){2D} \end{overpic} Figure 1.18: Coupling 2D and 1D domain.

####  1.5.2 Methodology

The basic idea of the first method is to consider a 1D formulation for the
straight channels and a single 2D element at the junctions, as shown in Fig.
1.18. It is important to explain how to switch from 1D to 2D model and vice-
versa.

The 2D flux is calculated by solving the rotate 1D Riemann problem (1.101) in
local coordinates, that reads as follows:

$\left\\{\begin{array}[]{c}\displaystyle\partial_{t}{\bf
Q}_{1D}+\partial_{s}{\bf F}_{1D}({\bf Q}_{1D})={\bf 0}.\\\ \displaystyle{\bf
Q}_{1D}(s,0)=\left\\{\begin{array}[]{c}\displaystyle{\bf Q}_{1D,L}={\bf
T}_{e}({\bf Q}_{2D,L})\quad{\rm if}\,s<0\\\ \displaystyle{\bf Q}_{1D,R}={\bf
T}_{e}({\bf Q}_{2D,R})\quad{\rm if}\,s>0\end{array}\right.\end{array}\right.$
(1.101)

where, the first line is the PDEs in normal direction, while the other
condition are the rotated initial conditions.

Near to the junctions, a 2D system of equations is required. Referring to Fig.
1.6(a) the result is a 1D system, and it is possible to obtain a non-zero
transversal components for the 1D element, as follows

$\displaystyle u_{1D}={\rm
sign}\left(u_{1D}\right)\sqrt{u_{2D}^{2}+v_{2D}^{2}}.$ (1.102)

###  1.6 Appendix

####  1.6.1 Shallow water equations

Shallow water theory is an approximation where the depth of the water is
assumed to be sufficiently small compared to some significant length, like the
radius of curvature of the water surface. This theory applies to many
occurrences in nature, such as flows in open channels, like roll waves, and
flood waves in rivers. It is commonly used in hydraulics by engineers dealing
with these applications. It can be shown that under such conditions, the
pressure in the fluid region is essentially the hydrostatic pressure.

The formal way to derive SW equations, starting from the exact hydrodynamical
system, is obtained as the lowest order approximation in the asymptotic
expansion procedure. The goal is to construct and justify the asymptotic
expansion of the solution in powers of the characteristic quantity defined as
the ratio between the typical length in $y$-direction and the one in $x$\-
(and $z$-) direction. In [18] Section 2.4, there is an explicit derivation of
the $0^{th}$ and $1^{st}$ orders of the asymptotic expansion, and, after some
integration, one is able to re-obtain the hydrostatic pressure relation that
is the analogue of the assumption made for the SW equations.

Here we provide a simple derivation of the shallow water equation under the
assumption of irrotational flow and hydrostatic approximation. We restrict to
consider a two dimensional problem in the $x$-$y$ vertical plane.

The $x$-axis is represented by the undisturbed free surface of the water, and
the $y$-axis is taken vertically upwards. Thus, the bottom expression is given
by $y=-h(x)$, where $h$ is the variable depth of the undisturbed water. The
variable $y=\eta(x,t)$ represents the surface displacement, and the velocity
components are $u(x,y,t)$ and $v(x,y,t)$.

Incompressible Euler equations in a vertical $x$-$y$ plane, in the domain
occupied by the fluid (see Figure 1.5) are given by

$\displaystyle u_{t}+uu_{x}+vu_{y}$ $\displaystyle=$ $\displaystyle-
p_{x}/\rho$ (1.103) $\displaystyle v_{t}+uv_{x}+vv_{y}$ $\displaystyle=$
$\displaystyle-p_{y}/\rho-g$ (1.104)

The continuity equation is

$u_{x}+v_{y}=0.$ (1.105)

It can be shown that, because of the different scaling in $x$ and $y$, the
vertical component $v$ of the velocity is much smaller than the horizontal one
$u$. Furthermore, in incompressible inviscid flow, vorticity
$\omega=v_{y}-u_{x}$ is carried unchanged aling material lines. As a
consequence, if the vorticity is zero at initial time, it will remain zero for
all time. As a consequence of these to facts, one can assume that $u_{y}$ is
negligible, therefore in the shallow water approximation the horizontal
velocity is independent of $y$. As a consequence, Eq. (1.103) becomes

$u_{t}+uu_{x}=g\eta_{x}(x,t)$

with $u=u(x,t)$.

At the free surface, the kinematic condition to be satisfied is

$(\eta_{t}+u\eta_{x}-v)\Big{|}_{y=\eta}=0$ (1.106)

and the dynamic condition is

$p\Big{|}_{y=\eta}=0.$ (1.107)

At the bottom, the only condition to be satisfied is

$(uh_{x}+v)\Big{|}_{y=-h}=0.$ (1.108)

Integrating Eq. (1.105) with respect to $y$, it becomes

$\int_{-h}^{\eta}(u_{x})dy+v\Big{|}_{-h}^{\eta}=0.$ (1.109)

Substituting the conditions 1.106 and 1.108 in Eq. (1.109), it follows that

$\int_{-h}^{\eta}u_{x}dy+\eta_{t}+u|_{\eta}\cdot\eta_{x}+u\Big{|}_{-h}\cdot
h_{x}=0.$ (1.110)

Now, this relation is valid

$\frac{\partial}{\partial
x}\int_{-h(x)}^{\eta(x)}udy=u\Big{|}_{y=\eta}\cdot\eta_{x}+u\Big{|}_{y=-h}\cdot
h_{x}+\int_{-h}^{\eta}u_{x}dy$ (1.111)

and combining it with Eq. (1.109), it becomes

$\frac{\partial}{\partial x}\int_{-h}^{\eta}udy=-\eta_{t}.$ (1.112)

Until here, the formulation is exact because no approximation has been
applied.

The shallow water approximation results from the assumption that the pressure
$p$ is given by

$p=g\rho(\eta-y)$ (1.113)

as in hydrostatics, where $\rho$ is the water density. In this way $p$ is
assumed to be independent from the $y$-component of the acceleration of the
water particles.

Deriving the pressure with respect to $x$ gives

$p_{x}=g\rho\eta_{x}$ (1.114)

and thus $p_{x}$ is independent of $y$. Therefore, it follows that the
acceleration of the water particles in the $x$ direction is also independent
of $y$, and so is $u$, the $x$-component of the velocity. So, from Eq. 1.114
the equation of motion in the $x$-direction becomes

$u_{t}+uu_{x}=-g\eta_{x}.$ (1.115)

This equation makes use of $u_{y}=0$, and Eq. 1.112 can be written as

$[u(\eta+h)]_{x}=-\eta_{t}$ (1.116)

having $\displaystyle\int^{\eta}_{-h}udy=u\int^{\eta}_{-h}dy$. Equations 1.115
and 1.116 are the first order differential equations that represent the non-
linear shallow water theory. If the initial conditions for $u$ and $\eta$ are
specified, the solution of these two equations describes the motion of the
fluid.

####  1.6.2 Laplace-Beltrami operator

Here we prove that

$\partial_{i}\tilde{\partial}_{i}c=\tilde{\partial}_{i}\tilde{\partial}_{i}c=:\Delta_{\perp}c$
(1.117)

where $\Delta_{\perp}$ is the Laplace-Beltrami operator. We also provide a
representation of $\Delta_{\perp}$ in Cartesian coordinates.

First, observe that

$\partial_{i}=\tilde{\partial}_{i}+n_{i}n_{j}\partial_{j},$

therefore

$\partial_{i}\tilde{\partial}_{i}c=\tilde{\partial}_{i}\tilde{\partial}_{i}c+n_{i}n_{j}\partial_{j}\tilde{\partial}_{i}c.$

Now consider the term

$n_{i}n_{j}\partial_{j}\tilde{\partial}_{i}c=n_{i}\frac{\partial}{\partial
n}\tilde{\partial}_{i}c.$

Since $n_{i}\tilde{\partial}_{i}c\equiv 0$, it is

$0=\frac{\partial}{\partial n}(n_{i}\tilde{\partial}_{i}c)=\frac{\partial
n_{i}}{\partial n}\tilde{\partial}_{i}c+n_{i}\frac{\partial}{\partial
n}\tilde{\partial}_{i}c.$

Now $\partial n_{i}/\partial n=n_{j}\partial_{j}n_{i}=0$, therefore

$n_{i}\frac{\partial}{\partial n}\tilde{\partial}_{i}c=0$

from which 1.117 follows.

Now we give a representation for $\Delta_{\perp}$. We start observing that
$\tilde{\partial}_{i}c=h_{ij}\partial_{j}c$ is orthogonal to $\Sigma_{B}$,
since $n_{i}h_{ij}=0$. It is

$\displaystyle\Delta_{\perp}c=\tilde{\partial}_{i}\tilde{\partial}_{i}c$
$\displaystyle=$
$\displaystyle\tilde{\partial}_{i}(h_{ij}\partial_{j}c)=\tilde{\partial}_{i}(\partial_{i}c-n_{i}n_{j}\partial_{j}c)$
(1.118) $\displaystyle=$ $\displaystyle
h_{ij}\partial^{2}_{ij}c-\tilde{\partial}_{i}n_{i}\frac{\partial c}{\partial
n}-n_{i}(\tilde{\partial}_{i}n_{j})\partial_{j}c-n_{i}n_{j}\tilde{\partial}_{i}\partial_{j}c$
(1.119)

Observe that

$\tilde{\partial}_{i}n_{i}=h_{ij}\chi_{ji}=\nabla\cdot\hat{n}$

where $\chi_{ij}\equiv\partial_{i}n_{j}$ is the second fundamental form of the
surface in Cartesian coordinates $\Sigma_{B}$ and is orthogonal to $\hat{n}$.
The last two terms are zero, because $n_{i}\chi_{ij}=0$ and
$n_{i}\tilde{\partial}_{i}\partial_{j}c=0$.

Therefore it follows:

$\Delta_{\perp}c=\tilde{\partial}_{i}\tilde{\partial}_{i}c=h_{ij}\partial^{2}_{ij}c-\chi_{ii}n_{j}\partial_{j}c$

We recall that the trace of the second fundamental form of a 2D surface in 3D
is the sum of the two principal curvatures, i.e. twice the mean curvature of
the surface.

##### Laplace-Beltrami operator in two dimensions

In two space dimensions (see Fig. 1.11) expressions considerably simplify.
Here we prove Eq. (1.80), i.e.

$\frac{\partial^{2}c}{\partial\tau^{2}}\equiv\tau_{j}\partial_{j}(\tau_{i}\partial_{i}c)=h_{ij}\partial^{2}_{ij}c-\chi_{ii}\frac{\partial
c}{\partial n}\equiv\Delta_{\perp}c.$ (1.120)

Indeed it is

$\frac{\partial^{2}c}{\partial\tau^{2}}=\tau_{j}\partial_{j}(\tau_{i}\partial_{i}c)=\tau_{i}\tau_{j}\partial^{2}_{ij}c+\tau_{j}(\partial_{j}\tau_{i})\partial_{i}c=\tau_{i}\tau_{j}\partial^{2}_{ij}c-\frac{1}{R}n_{i}\partial_{i}c$
(1.121)

the latter equality is a consequence of the relation

$\frac{\partial\hat{\tau}}{\partial\tau}=-\frac{1}{R}\hat{n}$

where $R$ denotes the local radius of curvature. Finally we observe that in 2D
it is

$h_{ij}\equiv\delta_{ij}-n_{i}n_{j}=\tau_{i}\tau_{j}$ (1.122)

as it can be easily shown by inspection checking that

$n_{i}n_{j}+\tau_{i}\tau_{j}=\delta_{ij}$

Using Eq. (1.121) and 1.122 in 1.120 completes the proof of relation 1.80.

###### Exercises 1

Implement the discretization of the Eqs. (1.63-1.65) in $\Omega^{0}=[0,1]$
using second order accurate Finite Differences scheme (see Chapter … for more
details). Note that Eq. (1.65) is already discretized in Eq. (1.81).

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures_new/exercise_V_new.png}
\put(1.0,40.0){(a)} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures_new/exercise_V1_V2.png}
\put(1.0,40.0){(b)} \end{overpic}

Figure 1.19: (a): Plot of the functional defined in Eq. (1.129). (b): Plot of
the functional defined in Eqs.(1.133-1.135).

Now consider a 2D domain, in which the potential is defined as a function
$W:=W(x,y)$ and it is constant in $y$-direction, i.e., $W(x,y)=V(x)$ (see Eq.
(1.49)). In this case the domain is defined as
$\Omega^{\varepsilon}=[-\varepsilon,1]\times[0,1]$. Discretize in space the
following problem with Finite Differences scheme

$\displaystyle\frac{\partial c}{\partial t}$
$\displaystyle=D\nabla\cdot\vec{J}\quad{\rm in}\,\Omega^{\varepsilon}$ (1.123)
$\displaystyle\vec{J}$ $\displaystyle=\nabla c+\frac{1}{k_{B}T}c\nabla W$
(1.124) $\displaystyle\vec{J}\cdot\vec{n}$ $\displaystyle=0\quad{\rm
on}\,\partial\Omega^{\varepsilon}$ (1.125)

where $\vec{n}$ is the outgoing normal vector to $\Omega^{\varepsilon}$.

Starting from the system (1.63-1.65)

$\displaystyle\displaystyle\frac{\partial c}{\partial t}+\frac{\partial
J}{\partial x}$ $\displaystyle=0\qquad{\rm in}\,\Omega^{\varepsilon}$ (1.126)
$\displaystyle\displaystyle J$ $\displaystyle=-D\left(\frac{\partial
c}{\partial x}+\frac{1}{k_{B}T}cV^{\prime}\right)$ (1.127) $\displaystyle J$
$\displaystyle=0\qquad{\rm on}\,\partial\Omega^{\varepsilon}$ (1.128)

where, in this case, the full domain is $\Omega=[0.5,1]$, with
$\Omega^{\varepsilon}_{b}=[0.5,0.5+L\varepsilon]$ the potential domain,
$\Omega^{\varepsilon}_{f}=[0.5+L\varepsilon,1]$ the fluid domain, and $L=2$.

Following the strategy in Section 1.3.1, find the new expression for $M$ (see
Eq. (1.66)) considering the following expression for the potential (see Fig.
1.19 (a))

$V(x)=a_{1}\exp\left(-b_{1}\left(\frac{x-x_{0}}{\varepsilon}\right)^{2}\right)-a_{2}\exp\left(-b_{2}\left(\frac{x-x_{0}}{\varepsilon}\right)^{2}\right)$
(1.129)

with $a_{1}=6,b_{1}=30,a_{2}=3,b_{2}=10,x_{0}=0.5$.

Starting from the system (1.63-1.65)

$\displaystyle\displaystyle\frac{\partial c}{\partial t}+\frac{\partial
J}{\partial x}$ $\displaystyle=0\qquad{\rm in}\,\Omega^{\varepsilon}$ (1.130)
$\displaystyle\displaystyle J$ $\displaystyle=-D\left(\frac{\partial
c}{\partial x}+\frac{1}{k_{B}T}cV^{\prime}\right)$ (1.131) $\displaystyle J$
$\displaystyle=0\qquad{\rm on}\,\partial\Omega^{\varepsilon}$ (1.132)

where, in this case, the potential $V$ is defined as follows (see Fig. 1.19
(b))

$\displaystyle V(x)=V_{1}(x)+V_{2}(x)$ (1.133) $\displaystyle
V_{1}(x)=f(x,x_{1}),\quad V_{2}(x)=f(x,x_{2}),\quad x_{1}=0,\quad x_{2}=1,$
(1.134)

and

$f(x,x_{l})=a_{1}\exp\left(-b_{1}\left(\frac{x-x_{l}}{\varepsilon}\right)^{2}\right)-a_{2}\exp\left(-b_{2}\left(\frac{x-x_{l}}{\varepsilon}\right)^{2}\right),$
(1.135)

with $l=1,2$, $a_{1}=6,b_{1}=30,a_{2}=3,b_{2}=10,x_{0}=0.5$. The full domain
is $\Omega=[-\varepsilon,1+\varepsilon]$, with
$\Omega^{\varepsilon}_{1,b}=\varepsilon[-1,L]$ and
$\Omega^{\varepsilon}_{2,b}=[1-L\varepsilon,1+\varepsilon]$ the potential
domains, $\Omega^{\varepsilon}_{f}=[L\varepsilon,1-L\varepsilon]$ the fluid
domain, and $L=2$.

Following the strategy in Section 1.3.1, find the new boundary conditions (see
Eqs. (1.65-1.66)).

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures_new/domain_level_set}
\put(1.0,40.0){(a)} \put(15.0,30.0){$\Omega^{0}$} \put(30.0,33.0){$\varphi<0$}
\put(29.5,17.0){$\Big{\uparrow}$} \put(27.0,12.0){$\varphi>0$}
\put(27.0,26.0){$R_{b}$} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures_new/exercise_V_squared.png}
\put(1.0,40.0){(b)} \end{overpic}

Figure 1.20: (a): Scheme of the level-set function. (b): Scheme of a potential
V(x) that represents a squared sink.

Let us consider the multiscale model in Eqs.(1.78-1.79), with
$\Omega^{0}=[-1,1]^{2}\setminus\mathcal{B}$ and a ball $\mathcal{B}$ centered
in $C=(0,0)$, with radius $R_{b}=0.25$. Using a level-set approach, it is
possible to implicitly know the boundaries of the ball (that are marked in red
in Fig. 1.20 (a)) defining a function
$\varphi=R_{b}-\sqrt{x^{2}+y^{2}},\,(x,y)\in[-1,1]^{2}$. With this technique,
the values of the level-set function are negative inside the domain and
positive outside of the domain (i.e., the points that are inside the ball; see
Fig. 1.20 (a)).

Given two balls in the domain $\Omega^{0}$, s.t.
$\Omega^{0}=[-1,1]^{2}\setminus(\mathcal{B}_{1}\cup\mathcal{B}_{2})$, with
$\mathcal{B}_{1}$ centered in $C_{1}=(-0.5,-0.5)$ and radius $R_{b,1}=0.2$,
and $\mathcal{B}_{2}$ centered in $C_{2}=(0.5,0.5)$ and radius $R_{b,2}=0.25$.
Determine a possible level set function $\varphi$, s.t. $\varphi>0$ in
$\mathcal{B}_{1}\cup\mathcal{B}_{2}$, and $\varphi<0$ in $\Omega^{0}$.

Given the initial conditions in Eq. (1.17), determine the constants $a,b$ in
Eq. (1.15).

Given the initial conditions in Eq. (1.17), determine the constants
$C_{+},C_{-}$ in Eq. (1.20).

Find the expression of the solution of Eq. (1.22) in the three case
$\Delta>0,\Delta=0,\Delta<0$.

Deduce the equation of motion of the pendulum starting from equations (1.43)
and (1.40).

Prove that Eq. (1.100) is valid.

###### Project 1

Find an analytical solution in 1D of the multiscale model in Eqs. (1.63-1.65),
with a squared sink in $x=0$, as it shows in Fig. 1.20 (b).

###### Project 2

In a 2D squared domain $\Omega^{0}=[-1,1]^{2}\setminus\mathcal{B}$, solve
numerically the system (1.78-1.79) with a squared hole $\mathcal{B}$ of
coordinates: $A=(-0.5,-0.5),\,B=(-0.5,0.5),\,C=(0.5,0.5),\,D=(0.5,-0.5)$.

## References

  * [1] B. Blonder, S. Both, M. Jodra, H. Xu, M. Fricker, I. Matos, N. Majalap, D. Burslem, Y. Teh, Yit and Y. Malhi, Linking functional traits to multiscale statistics of leaf venation networks, New Phytologist, 2020
  * [2] C. Astuto, D. Boffi, J. Haskovec, P. Markowich, and G. Russo, _Comparison of Two Aspects of a PDE Model for Biological Network Formation_ , Mathematical and Computational Applications, 2022, 10.3390/mca27050087
  * [3] C. Astuto, D. Boffi, J. Haskovec, P. Markowich, and G. Russo, _Asymmetry and condition number of an elliptic-parabolic system for biological network formation_ , Communications on Applied Mathematics and Computation, 2023, arXiv:2301.12926
  * [4] C. Astuto, D. Boffi, and F. Credali,  _Finite element discretization of a biological network formation system: a preliminary study_ , 2023, arXiv:2303.10625
  * [5] Saint-Venant, A.J.C. Barré de , Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et a l’introduction de marées dans leurs lits, Comptes Rendus de l’Académie des Sciences, 73: 147–154 and 237–240, (1871).
  * [6] Saint-Venant, B de Memories des Savants Etrangers, Vol. 14, Paris Academie des Sciences, Paris, France (in French), 1855
  * [7] Synge, JL , The problem of Saint Venant for a cylinder with free sides, Quarterly of Applied Mathematics, 1945
  * [8] Boffi, Daniele and Credali, Fabio and Gastaldi, Lucia, _On the interface matrix for fluid–structure interaction problems with fictitious domain approach_ , Computer Methods in Applied Mechanics and Engineering, 2022
  * [9] Quarteroni, Alfio and Veneziani, Alessandro and Vergara, Christian, Geometric multiscale modeling of the cardiovascular system, between theory and practice Computer Methods in Applied Mechanics and Engineering, Elsevier, 2016
  * [10] Townsley, Mary I, Structure and composition of pulmonary arteries, capillaries and veins, Comprehensive Physiology, 2012
  * [11] Quarteroni, Alfio and Manzoni, Andrea and Vergara, Christian, The cardiovascular system: mathematical modelling, numerical algorithms and clinical applications, Acta Numerica, 2017
  * [12] Bellamoli, Francesca and Müller, Lucas O and Toro, Eleuterio F A numerical method for junctions in networks of shallow-water channels, Applied Mathematics and Computation, 2018
  * [13] Brochard, F and Lennon, JF, Frequency spectrum of the flicker phenomenon in erythrocytes, Journal de Physique, 1975
  * [14] Raudino, A. and Raciti, D. and Grassi, A. and Pannuzzo, M. and Corti, M., Oscillations of Bubble Shape Cause Anomalous Surfactant Diffusion: Experiments, Theory, and Simulations, Langmuir, 2016
  * [15] C. Astuto, A. Raudino and G. Russo, Multiscale modeling of sorption kinetics, Multiscale Modeling & Simulation (SIAM), 2023, https://doi.org/10.1137/21M1463872
  * [16] Dukhin, Stanislav Samuilovich and Kretzschmar, Günter and Miller, Reinhard, Dynamics of adsorption at liquid interfaces: theory, experiment, application, Elsevier, 1995
  * [17] Nguyen, Phuc Nghia and Waton, Gilles and Vandamme, Thierry and Krafft, Marie Pierre, Behavior of an adsorbed phospholipid monolayer submitted to prolonged periodical surface density variations, Wiley Online Library, 2013
  * [18] Stoker, James Johnston, Water waves: The mathematical theory with applications, John Wiley & Sons, 1992
  * [19] Vignoli, G. and Titarev, V.A. and Toro, E.F., ADER schemes for the shallow water equations in channel with irregular bottom elevation, Journal of Computational Physics, 2008
  * [20] Ata, R. and Pavan, S. and Khelladi, S. and Toro, E.F ,A Weighted Average Flux (WAF) scheme applied to shallow water equations for real-life applications, Advances in Water Resources, 2013
  * [21] Perktold, K and Hilbert, D, Numerical simulation of pulsatile flow in a carotid bifurcation model, Journal of biomedical engineering, 1986
  * [22] Formaggia, Luca and Quarteroni, Alfio and Veneziani, Alessandro, Multiscale models of the vascular system, Springer, 2009
  * [23] Euler, Leonhard, Principes généraux du mouvement des fluides, Mémoires de l’académie des sciences de Berlin, 1757
  * [24] Euler, Leonard, Principia pro motu sanguinis per arterias determinando. Opera Postuma, 2, 814-823, 1862.
  * [25] Coco, Armando, A multigrid ghost-point level-set method for incompressible Navier-Stokes equations on moving domains with curved boundaries, Journal of Computational Physics, 2020
  * [26] C. Astuto, A. Coco and Giovanni Russo, A finite–difference ghost–point multigrid method for multi–scale modelling of sorption kinetics of a surfactant past an oscillating bubble, Journal of Computational Physics, 2023, https://doi.org/10.1016/j.jcp.2022.111880
  * [27] Sussman, Mark and Smereka, Peter and Osher, Stanley, A level set approach for computing solutions to incompressible two-phase flow, Journal of Computational physics, 1994
  * [28] Sethian, James Albert, Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science, Cambridge university press, 1999
  * [29] Russo, Giovanni and Smereka, Peter, A remark on computing distance functions, Journal of computational physics, 2000
  * [30] Osher, Stanley and Fedkiw, Ronald and Piechor, K, Level set methods and dynamic implicit surfaces, Appl. Mech. Rev., 2004
  * [31] Armando Coco and Giovanni Russo , Finite-difference ghost-point multigrid methods on Cartesian grids for elliptic problems in arbitrary domains, Journal of Computational Physics, 2013
  * [32] Armando Coco and Giovanni Russo, Second order finite-difference ghost-point multigrid methods for elliptic problems with discontinuous coefficients on an arbitrary interface, Journal of Computational Physics, 2018
  * [33] A. Raudino, A. Grassi, G. Lombardo, G. Russo , Giovanni C. Astuto and M. Corti, Anomalous Sorption Kinetics of Self-Interacting Particles by a Spherical Trap, Communications in Computational Physics, 2022, https://doi.org/10.4208/cicp.OA-2021-0101
  * [34] Uri M. Ascher and Linda R. Petzold, Computer methods for ordinary differential equations and differential-algebraic equations, Volume 61, SIAM, 1998.

## Index

  * Coupling of 3D-1D, 3D-0D and 1D-0D models §1.4.3
  * Multiscale
    * Modeling Chapter 1
  * The 1D and 0D model §1.4.2
  * The 3D model §1.4.1