# Device-independent Quantum Fingerprinting for Large Scale Localization

Ahmed Shokry  Dept. of Computer Science and Engineering
American University in Cairo
Cairo, Egypt
<EMAIL_ADDRESS>Moustafa Youssef  Dept. of Computer Science and
Engineering
American University in Cairo & Alexandria University
Cairo, Egypt
<EMAIL_ADDRESS>
###### Abstract

Although RF fingerprinting is one of the most commonly used techniques for
localization, deploying it in a ubiquitous manner requires addressing the
challenge of supporting a large number of heterogeneous devices and their
variations. We present QHFP, a device-independent quantum fingerprint matching
algorithm that addresses two of the issues for realizing worldwide ubiquitous
large-scale location tracking systems: storage space and running time as well
as devices heterogeneity. In particular, we present a quantum algorithm with a
complexity that is exponentially better than the classical techniques, both in
space and running time. QHFP also has provisions for handling the inherent
localization error due to building the large-scale fingerprint using
heterogeneous devices. We give the details of the entire system starting from
extracting device-independent features from the raw RSS, mapping the classical
feature vectors to their quantum counterparts, and showing a quantum cosine
similarity algorithm for fingerprint matching.

We have implemented our quantum algorithm and deployed it in a real testbed
using the IBM Quantum machine simulator. Results confirm the ability of QHFP
to obtain the correct estimated location with an exponential improvement in
space and running time compared to the traditional classical counterparts. In
addition, the proposed device-independent features lead to more than 20%
better accuracy in median error. This highlights the promise of our algorithm
for future ubiquitous large-scale worldwide device-independent fingerprinting
localization systems.

###### Index Terms:

quantum computing, device-independent location determination systems,
practical quantum algorithms, quantum location determination, next generation
location tracking systems, quantum pervasive algorithms and systems.

##  I Introduction

Recent years have witnessed the advent of the field of quantum computing (QC),
with a variety of large corporations and startups investing in quantum
computing [15, 14, 27]. QC algorithms can solve problems that are infeasible
to be solved by classical computers [17], providing exponential gains in some
cases [25]. This opens the door for investigating quantum algorithms
advantages in new fields such as worldwide location tracking.

Large-scale worldwide fingerprinting localization are commonly used indoor [7,
29, 2, 6] and outdoor [22, 20, 5, 1, 3] due to their accuracy. The
fingerprinting techniques work in two phases: the offline RF fingerprint
building phase and the online tracking phase. In the offline phase, the
received signal strength (RSS) coming from the different reference points
(RPs111These reference points can be, e.g., WiFi access points; cellular cell
towers; or Bluetooth beacons.) in the environment are recorded at the
different discrete locations. In the online phase, the RSS coming from RPs at
unknown locations are matched aganist the fingerprint and the closest location
in the RSS space becomes the estimated location. However, when the training
devices that are used to build the fingerprint and the test devices that are
used for localization during the online phase are different, the system
accuracy severely degrades. Hence, deploying such systems on different phone
types is not a straightforward task as the RSS readings vary for the different
kinds of phones, even at the same location and time. This problem, which is
known as the device heterogeneity problem [13, 11, 12, 18, 21], can prevent
ubiquitous large-scale fingerprinting-based localization. Moreover, all the
traditional fingerprinting techniques need to match the online RSS
measurements from the different RPs to the RSS measurements at each
fingerpirnt location, making their time and space complexity $o(MN)$, where
$M$ is the number of fingerprint locations and $N$ is the number of RPs. This
can hinder deploying the traditional outdoor/indoor localization in a the
large-scale worldwide settings, especially for IoT environments, where the
number of RPs in an environment can be significant.

In this paper, we present QHFP: a practical device-independent quantum
fingerprint-based location tracking algorithm that requires space and runs in
$o(M\log(N))$. We show the details of how to extract device-independent
features from the received signal strength measurements, how to construct the
quantum fingerprint, how to encode the extracted features in quantum states,
and finally; how to calculate the quantum similarity between the online
features and the offline ones stored in the fingerprint.

We validate our quantum algorithm in a real testbed. In addition, we quantify
its performance using simulations on an IBM quantum computer. The results show
the ability of QHFP to correctly obtain the estimated location with the same
accuracy as its classical counterparts. This comes with an exponential
enhancement compared to the traditional classical fingerprinting techniques in
both space and time. In addition, QHFP provides 20% enhancement in median
localization error in addressing the heterogeneity problem.

The remaining sections are organized as follow: we begin with a brief
background on quantum computing in Section II. Section III provides the
details of our device-independent quantum fingerprint localization algorithm.
Details of the implementation and evaluation results are presented and
discussed in Sections IV. Finally, Section V concludes the paper.

##  II Background on Quantum Computing

In this section, we give a brief background on the basic concepts of quantum
computing that we will build on in the rest of the paper [23, 24].

A quantum bit (qubit) is the basic unit of information and is analogue to the
classical bit. Contrary to classical bits, a qubit can exist in a
superposition of the zero and one states. This superposition is what allows
quantum computations to work on both states at the same time. This is often
referred to as quantum parallelism. Qubits can have various physical
implementations, e.g. the polarization of photons.

Formally, the Dirac notation is commonly used to describe the state of a qubit
as $\ket{\psi}=\alpha\ket{0}+\beta\ket{1}$, where $\alpha$ and $\beta$ are
complex numbers called the amplitudes of classical states $\ket{0}$ and
$\ket{1}$, respectively. The state of the qubit is normalized, i.e.
$\alpha^{2}+\beta^{2}=1$. When the state $\ket{\psi}$ is measured, only one of
$\ket{0}$ or $\ket{1}$ is observed, with probability $\alpha^{2}$ and
$\beta^{2}$, respectively. The measurement process is destructive, in the
sense that the state collapses to the value $\ket{0}$ or $\ket{1}$ that has
been observed, losing the original amplitudes $\alpha$ and $\beta$ [17].

Operations on qubits are usually represented by gates, similar to a classical
circuit. An example of a common quantum gate is the NOT gate (also called
Pauli-X gate) that is analogous to the not gate in classical circuits. In
particular, when we apply the NOT gate to the state
$\ket{\psi_{0}}=\alpha\ket{0}+\beta\ket{1}$, we get the state
$\ket{\psi_{1}}=\beta\ket{0}+\alpha\ket{1}$. Gates are usually represented by
unitary matrices while states are represented by column vectors222The ket
notation $\ket{.}$ is used for column vectors while the bra notation $\bra{.}$
is used for row vectors.. The matrix for the NOT gate is
$\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}$ and the above operation can be
written as $\ket{\psi_{1}}=NOT(\ket{\psi_{0}})=\begin{bmatrix}0&1\\\
1&0\end{bmatrix}\begin{bmatrix}\alpha\\\ \beta\end{bmatrix}.$

Another important gate is the Walsh–Hadamard gate, $H$, that maps $\ket{0}$ to
$\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$, i.e. a superposition state with equal
probability for $\ket{0}$ and $\ket{1}$; and maps $\ket{1}$ to
$\frac{1}{\sqrt{2}}(\ket{0}-\ket{1})$. Figure 1 shows a simple quantum
circuit. Single lines carry quantum information while double lines carry
classical information (typically after measurement). The simple circuit
applies an $H$ gate to state $\ket{0}$, which produces the state
$\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$ at the output of the gate. The
measurement step produces either 0 or 1 with equal probability (the squared
amplitude of the measured state). The state collapses to the observed
classical bit value.

It is important to note that the concept of quantum interference is at the
core of quantum computing. Using quantum interference, one uses gates to
cleverly and intentionally bias the content of the qubits towards the needed
state, hence achieving a specific computation result.

The notion of qubit can be extended to higher dimensions using a quantum
register. A quantum register $\ket{\psi}$, consisting of $n$ qubits, lives in
a $2^{n}$-dimensional complex Hilbert space $\mathcal{H}$. Register
$\ket{\psi}=\sum_{0}^{2^{n}-1}\alpha_{i}\ket{i}$ is specified by complex
numbers $\alpha_{0},...,\alpha_{2^{n}-1}$, where $\sum|\alpha_{i}|^{2}=1$.
Basis state $\ket{i}$ denotes the binary encoding of integer $i$. We use the
tensor product $\otimes$ to compose two quantum systems. For example, we can
compose the two quantum states $\ket{\psi}=\alpha\ket{0}+\beta\ket{1}$ and
$\ket{\phi}=\gamma\ket{0}+\delta\ket{1}$ as
$\ket{\omega}=\ket{\psi}\otimes\ket{\phi}=\alpha\gamma\ket{00}+\alpha\delta\ket{01}+\beta\gamma\ket{10}+\beta\delta\ket{11}$.

Figure 1: An example of quantum circuit. Single lines represent quantum
information and double lines represent classical information.

(a) The CNOT gate.

(b) The CSWAP gate.

Figure 2: Simple controlled gates. (a) The NOT gate is applied to the target
qubit, if and only if, the source qubit is $\ket{1}$. (b) The two target
qubits are swapped, if and only if, the source qubit is $\ket{1}$.

Controlled gates act on multiple qubits, where one or more qubits act as a
control for some operation on the other qubits (Figure 2). Figure 2(a)
illustrates the controlled NOT gate (CNOT). When the source qubit is
$\ket{1}$, the NOT operation will be applied to the target qubit. Figure 2(b)
further shows the controlled SWAP gate (CSWAP) with three qubits as input. In
the CSWAP gate, the swap operation is performed on the target wires, if and
only if, the source line is $\ket{1}$. This can be used to “entangle” qubits
together. Entangled qubits are correlated with one another, in the sense that
information on one qubit will reveal information about the other unknown
qubit, even if they are separated by large distance [17].

A common way to describe a quantum algorithm is to use a quantum circuit,
which is a combination of the quantum gates (e.g as in Figure 1). The input to
the circuit is a number of qubits (in quantum registers) and the gates act on
them to change the combined circuit state using superposition, entanglement,
and interference to reach a desired output state that is a function of the
algorithm output. The final step is to measure the output state(s), which
reveals the required information.

##  III QHFP: A Device-independent Quantum Fingerprinting Localization System

In this section, we present the detailed description of our QHFP quantum
device-independent fingerprinting localization algorithm.

The fingerprinting techniques work in two phases: the offline RF fingerprint
building phase and the online tracking phase. In the offline phase, the RSS
coming from the different RPs in the environment are recorded at the different
discrete locations. In the online phase, the RSS coming from RPs at unknown
locations are matched aganist the fingerprint and the closest location in the
RSS space becomes the estimated location. Hence, fingerprint localization can
capture the relation between RSS coming from the different RPs in the
environment and user location [30, 31, 16].

There are different metrics that can be used to match the online RSS with the
fingerprint. In our paper, we proposed a quantum algorithm that uses cosine
similarity, which is one of the popular approaches usually used to mitigate
device heterogeneity effects [10, 28].

We start the section by explaining how to construct the device-independent
fingerprint and how to obtain the quantum fingerprint through encoding the
offline RSS vectors from the device-independent fingerprint and the online RSS
vector in qubits (i.e. state preparation phase). Then, we present the details
of the quantum fingerprinting matching algorithm. Finally, we give an example
on how the algorithm can work in detail.

###  III-A Device-independent Fingerprint Transformation

In this section, we describe two different techniques to handle the devices
heterogeneity problem: the power ratio and the power difference. The basic
idea is to transform the traditional RSS fingerprint to a device-independent
fingerprint.

####  III-A1 Power ratio

This technique assumes that the ratio between the RSS values coming from
different cell towers remains the same on the different phones. Hence, it
transforms the RSS values (i.e. power) to a relative power. That is, instead
of using the raw RSS from individual RPs in the fingerprint building phase
(offline phase) and the raw RSS in the online matching, it uses the ratio of
the RP RSS readings from each pair of RPs.

More formally, given a RSS features vector $X$= $(f_{1},f_{2},..,f_{N})$ from
$N$ RPs in the environment, where $f_{i}$ is the RSS coming from RP $i$, this
technique transforms $X$ to $X_{r}$= $(r_{1,2},r_{1,3},..,r_{N-1,N})$, where
$r_{i,j}=\frac{f_{i}}{f_{j}}$ and $|X_{r}|={N\choose 2}$.

####  III-A2 Power difference

This approach is similar to the previous approach but takes the power
difference instead of the power ratio. It assumes that the difference between
the RSS values coming from different RPs remains the same on the different
phones.

More formally, given a RSS features vector $X$= $(f_{1},f_{2},..,f_{N})$ from
$N$ cell-towers in the environment, where $f_{i}$ is the RSS coming from cell
tower $i$, this technique transforms $X$ to $X_{d}$=
$(d_{1,2},d_{1,3},..,d_{N-1,N})$, where $d_{i,j}=f_{i}-f_{j}$ and
$|X_{d}|={N\choose 2}$.

In the next subsection, we show how to map these relative classical feature
vectors into a quantum register.

###  III-B State Preparation

The state preparation step aims to encode the relative classical RSS vectors,
coming from the power ratio/difference, in quantum registers. Next, the
quantum fingerprinting matching algorithm is applied between the online RSS
register and the offline ones. Assume there are $N$ RPs in the environment.
Therefore, the ${N\choose 2}$-dimensional normalized RSS vector from the $N$
RPs $v=(\beta_{0},\beta_{1},...,\beta_{{N\choose 2}-1}),\sum_{i=0}^{{N\choose
2}-1}\beta_{i}^{2}=1$, can be encoded using a quantum register $\ket{\delta}$
with $n=log({N\choose 2})$ qubits. Where $\ket{\delta}=\sum_{i=0}^{{N\choose
2}-1}\alpha_{i}\ket{i}$, and the basis state $\ket{i}$ represents the binary
encoding of integer $i$ [9, 26]. Note that the ${N\choose 2}$-dimensional RSS
vector can be encoded in $log({N\choose 2})$ qubits, which is an exponential
saving in the space.

To perform the state preparation, we used the quantum random access memory
(QRAM) which is a popular technique in which the binary code of the RSS
measurements vector (i.e. $\beta_{i}$) is loaded into a qubit register in
parallel and conditional rotations are performed to the qubits in the register
in order to encode the RSS measurements as amplitudes in the quantum register
[9, 19, 32, 26].

We give an example of how to prepare this state from classical vectors later
in this section.

Figure 3: Quantum fingerprint matching circuit between a single fingerprint
RSS vector (encoded in qubit $\gamma$) and the online RSS vector (encoded in
qubit $\delta$). Figure 4: A detailed example of the quantum fingerprint
matching circuit using two-values RSS vectors. The circuit shows the state
preparation step, i.e. how to map the testing RSS vector $(0.43,0.9)$ and
training RSS vector $(0.24,0.97)$ to a quantum state, starting from $\ket{0}$.

###  III-C The Quantum Cosine Similarity Algorithm

Figure 3 shows the quantum circuit for calculating the cosine similarity
between two normalized RSS vectors coming from applying the RSS power
ratio/difference encoded in the quantum registers $\ket{\delta}$ (e.g. a test
RSS vector during the online location tracking phase) and $\ket{\gamma}$ (e.g.
a single fingerprint RSS vector) based on the CSWAP gate [4, 23, 24]. In
particular, the circuit calculates:

$\textrm{sim}(\ket{\delta},\ket{\gamma})=\cos[2](\delta,\gamma)=|\bra{\delta}\ket{\gamma}|^{2}$
(1)

where $\cos(\delta,\gamma)$ is cosine the angle between the two normalized
vectors $\delta$ and $\gamma$. The circuit in Figure 3 takes its input which
is the ancilla qubit at state $\ket{1}$ and the two quantum RSS vectors
encoding $\ket{\delta}$ and $\ket{\gamma}$. Then, it applies a series of gates
to transform the input to the following joint state,

$\frac{1}{2}\left(\sqrt{2-2|\bra{\delta}\ket{\gamma}|^{2}}\ket{0}\frac{\ket{\delta}\ket{\gamma}-\ket{\gamma}\ket{\delta}}{\sqrt{2-2|\bra{\delta}\ket{\gamma}|^{2}}}\right.\\\
+\left.\sqrt{2+2|\bra{\delta}\ket{\gamma}|^{2}}\ket{1}\frac{\ket{\delta}\ket{\gamma}+\ket{\gamma}\ket{\delta}}{\sqrt{2+2|\bra{\delta}\ket{\gamma}|^{2}}}\right)$
(2)

Finally, the probability of measuring the top (ancilla) qubit to be 1 is
$\frac{1}{2}(1+|\bra{\delta}\ket{\gamma}|^{2})$, which is a function of the
required similarity measure between the two vectors. We repeat this circuit
$K$ times to estimate the cosine similarity as $2\times\\#\ket{1}/K-1$.
Algorithm 1 explains our device-independent quantum fingerprinting matching
algorithm. The $H$-gate is applied to the ancilla qubit in order to put it in
a superposition state. The CSWAP gate is applied to entangle the ancilla qubit
with the training and testing quantum registers. Finally, the $H$-gate is
applied to generate the desired computation where the probability of receiving
$\ket{1}$ for the ancilla qubit is a function of the required similarity.

###  III-D Example

In this section, we illustrate the quantum fingerprint matching algorithm
described in the previous section using a simple example with two-values RSS
vector. The two normalized RSS vectors to be matched are $v_{1}=(0.43,0.9)$
and $v_{2}=(0.24,0.97)$. Figure 4 shows the complete circuit for calculating
the cosine similarity between $v_{1}$ and $v_{2}$.

The circuit starts by the state preparation stage, i.e. mapping the RSS
vectors $v_{1}$ and $v_{2}$ to the quantum equivalent
$\ket{\delta}=0.43\ket{0}+0.9\ket{1}$ and
$\ket{\gamma}=0.24\ket{0}+0.97\ket{1}$, respectively. This is achieved by
using the $U$ gate, which is represented as,

$U(\theta)=\begin{bmatrix}\cos(\theta/2)&-\sin(\theta/2)\\\
\sin(\theta/2)&\cos(\theta/2)\end{bmatrix}$ (3)

Where $\theta$ is double the angle between $\ket{0}$ and the quantum
representation of the normalized RSS vector $\ket{v}$.

The fingerprint matching part is the same as the one described in Algorithm 1.
Since we have only two values in the RSS vectors, the quantum registers
contain only one qubit. For the given example, the probability of measuring
the first qubit to be in state 1 is 0.9765 and hence the similarity score is
0.9529.

Note that the proposed quantum circuit needs to be repeated for each of the
fingerprint locations to determine its similarity score to the test RSS
vector. The fingerprint location with the highest score becomes the estimated
user location.

Algorithm 1 QHFP Fingerprint Matching

1:

2:1- Two n-qubits quantum registers $\ket{\delta}$ and $\ket{\gamma}$, storing
RSS vectors of the FP and test locations, coming from the power
ratio/difference techniques, to be compared. $n=log({N\choose 2})$, where $N$
is the number of RPs.

3:2- An ancilla qubit $\ket{a}=\ket{1}$

4:3- Number of iterations $K$.

5:

6:Compute an estimate of the similarity between $\ket{\delta}$ and
$\ket{\gamma}$ as $|\bra{\delta}\ket{\gamma}|^{2}$

7:for $k\leftarrow 1$ to $K$ do

8: Apply $H(\ket{a})$

9: for $i\leftarrow 1$ to $n$ do

10: Apply CSWAP($\ket{a}$,$\ket{\delta_{i}}$,$\ket{\gamma_{i}}$)

11: end for

12: Apply $H(\ket{a})$

13: $\eta_{k}$ $\leftarrow$ measurement of $\ket{a}$.

14:end for

15:return $\frac{2}{K}\sum_{k=1}^{K}\eta_{k}-1$

##  IV Implementation and Evaluation

In this section, we implement the proposed quantum localization algorithm and
evaluate its performance in a real testbed.

###  IV-A Experiment setup

Figure 5 shows our real cellular testbed that spans an 0.2Km2 outdoor urban
area. The area is covered by 8 different cell-towers. Data is collected by
war-driving uniformly over the entire area of interest. We also collected an
independent test set. Both the fingerprint and test locations are uniformly
distributed over the entire area of interest. We use different Android devices
for data collection including HTC Nexus One, Prestigio Multipad Wize 3037 3G,
HTC One X9 and Motorola Moto G5 plus phones among others. The deployed
collector software collects GPS ground-truth locations, received signal
strengths, and timestamps.

The 8 cell-towers RSS can be encoded in 5 qubits after applying RSS power
ratio or difference techniques. Therefore, we require a total of 11 qubits for
running the circuit in Figure 3, 10 qubits for encoding the training and
testing RSS vectors and one for the Ancilla qubit. We use the IBM quantum
machine simulator, which supports up to 28 qubits.

Figure 5: The outdoor testbed.

###  IV-B Comparison Between the Different Heterogeneity Handling Techniques

Figure 6 and Table I show a comparison between the traditional raw RSS, the
quantum power ratio, and the quantum power difference fingerprinting
techniques. Specifically, the results show that QHFP, through its
heterogeneity handling techniques, can achieve a better localization accuracy
than traditional raw RSS localization by more than 20% in the median error.
The results also show that the power difference technique is superior to the
power ratio technique. This can be explained by noting that the RSS reported
by the phone API is in log-scale.

Figure 6: Effect of device heterogeneity techniques. TABLE I: Comparison between quantum and classical localization systems error quantiles. Numbers between parenthesis indicate percentage of enhancement over raw RSS. System |  Q1 |  Median |  Q3
---|---|---|---
Raw RSS - Traditional |  40 |  103 |  121
QHFP \- Power Ratio |  39 (2.6%) |  98 (5.1%) |  121 (0%)
QHFP \- Power Difference |  36 (11.1%) |  85.5 (20.5%) |  119 (1.7%)

###  IV-C Comparison with Classical Localization

Figure 7 shows the CDF of distance error for the quantum and classical
localization system using the same device-independent features (i.e. after
applying the power difference heterogeneity technique). The figure shows that
QHFP has the same accuracy as the classical one, but with the potential
exponential saving in time and space.

Figure 7: CDF of localization error for quantum and classical localization
using the power difference technique.

###  IV-D Effect of Number of Shots

Figure 8 shows the impact of increasing the number of algorithm iterations,
i.e. re-running the system (parameter $K$ in Algorithm 1), on the power
difference quantum localization accuracy. It is evident from the figure that
increasing the number of iterations leads to a better localization accuracy
for both techniques till it starts saturating around 4096 shots.

Figure 8: Effect of number of shots on the quantum power difference
localization error.

###  IV-E Discussion

Last section shows that the device-independent fingerprinting localization
techniques have accuracy that is better than the traditional raw RSS one in
the real testbed used in our experiments.

The device-independent classical version of the localization system requires
$o(N^{2}M)$ space and matching runs in $o(N^{2}M)$, for $N$ RPs and $M$
fingerprint locations. On the other hand, QHFP requires $o(M\log N)$ for both
space and time. This is an exponential enhancement in both space and running
time. Figure 9 explains how the running time complexity for the quantum and
the classical localization systems scale with the increase of the number of
RPs for a fixed number of locations in the fingerprint. The figure confirms
the exponential saving in running time of the proposed quantum algorithm.

Figure 9: Running time complexity for the quantum and the classical
algorithms.

This is important in many scenarios such as worldwide localization and IoT
applications, where there are alot of RPs with different technologies
(cellular, WiFi, BT, etc). Moreover, in order to get a better accuracy, one
needs to fuze the signal coming from the significant number of RPs in the
environment. Hence, the exponential saving in time and space becomes a must.

QHFP is not only useful for the online fingerprint matching process but also
for offline state preparation (quantum fingerprint construction) as it can
reduce the size of the fingerprint RSS vectors from $o(N^{2})$ to o(log(N)).
This is important not only to speed up the matching process but also to save
the space required for downloading and storing the large-scale worldwide
fingerprint from the server to mobile phones if the fingerprint matching is
performed on the mobile phones. In addition, quantum co-processors (similar to
GPUs) have started to appear [8] which can enable quantum localization on
mobile devices.

Finally, the space and time complexity of QHFP can be reduced to $o(\log(NM))$
by encoding all the fingerprint locations data in one circuit. This might have
a significant implications to several localization and spatial algorithms.

##  V Conclusion

In this paper, we have presented, QHFP, a practical device-independent quantum
fingerprinting algorithm. QHFP can provide high accuracy localization with
different phone types while leading to an exponential saving in the space and
the running time of the current fingerprinting localization systems. Results
from deploying QHFP in a real testbed confirm its advantages compared to the
traditional classical techniques. Moreover, the proposed device-independent
features achieves a median distance error of 85.5m, which is 20% better than
the traditional fingerprinting methods.

Currently, we are exploring different quantum similarity metrics, obtaining
theoretical quantum bounds on performance, and optimizing other spatial
algorithms using quantum computing.

## References

  * [1] Heba Aly, Anas Basalamah, and Moustafa Youssef.  Lanequest: An accurate and energy-efficient lane detection system.  In 2015 IEEE International Conference on Pervasive Computing and Communications (PerCom), pages 163–171. IEEE, 2015.
  * [2] Heba Aly and Moustafa Youssef.  New insights into wifi-based device-free localization.  In Proceedings of the 2013 ACM conference on Pervasive and ubiquitous computing adjunct publication, pages 541–548, 2013.
  * [3] Heba Aly and Moustafa Youssef.  semMatch: Road semantics-based accurate map matching for challenging positioning data.  In Proceedings of the 23rd SIGSPATIAL International Conference on Advances in Geographic Information Systems, pages 1–10, 2015.
  * [4] Harry Buhrman, Richard Cleve, John Watrous, and Ronald De Wolf.  Quantum Fingerprinting.  Physical Review Letters, 87(16):167902, 2001.
  * [5] Gu Chen, Ahmed Shokry, and Moustafa Youssef.  The Effect of Ground Truth Accuracy on the Evaluation of Localization Systems.  IEEE INFOCOM, 2021.
  * [6] Rizanne Elbakly and Moustafa Youssef.  A robust zero-calibration RF-based localization system for realistic environments.  In 2016 13th Annual IEEE International Conference on Sensing, Communication, and Networking (SECON), pages 1–9. IEEE, 2016.
  * [7] Moustafa Elhamshary and Moustafa Youssef.  Towards Ubiquitous Indoor Spatial Awareness on a Worldwide Scale.  SIGSPATIAL Special, 9(2):36–43, 2017.
  * [8] George J Frangou, Stephane Chretien, and Ivan Rungger.  The first quantum co-processor hybrid for processing quantum point cloud multimodal sensor data.  Proceedings of the Future Technologies Conference, pages 411–426, 2019.
  * [9] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Quantum Random Access Memory.  Phys. Rev. Lett., 100:160501, Apr 2008.
  * [10] S. Han, C. Zhao, W. Meng, and C. Li.  Cosine Similarity based Fingerprinting Algorithm in WLAN Indoor Positioning against Device Diversity.  2015 IEEE International Conference on Communications (ICC), pages 2710–2714, 2015.
  * [11] Mohamed Ibrahim and Moustafa Youssef.  CellSense: A probabilistic RSSI-based GSM positioning system.  In 2010 IEEE Global Telecommunications Conference GLOBECOM 2010, pages 1–5. IEEE, 2010.
  * [12] Mohamed Ibrahim and Moustafa Youssef.  CellSense: An accurate energy-efficient GSM positioning system.  IEEE Transactions on Vehicular Technology, 61(1):286–296, 2011\.
  * [13] Mohamed Ibrahim and Moustafa Youssef.  Enabling wide deployment of GSM localization over heterogeneous phones.  In 2013 IEEE International Conference on Communications (ICC), pages 6396–6400. IEEE, 2013.
  * [14] Julian Kelly.  A preview of Bristlecone, Google’s New Quantum Processor.  Google Research Blog, 5, 2018.
  * [15] Will Knight.  IBM raises the Bar with a 50-qubit Quantum Computer.  Sighted at MIT Review Technology, 2017.
  * [16] Nesma Mohssen, Rana Momtaz, Heba Aly, and Moustafa Youssef.  It’s the human that matters: Accurate user orientation estimation for mobile computing applications.  Mobiquteous, 2014.
  * [17] Michael A. Nielsen and Isaac L. Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, USA, 10th edition, 2011.
  * [18] Jun-geun Park, Dorothy Curtis, Seth Teller, and Jonathan Ledlie.  Implications of Device Diversity for Organic Localization.  pages 3182–3190. IEEE, 2011.
  * [19] Patrick Rebentrost, Masoud Mohseni, and Seth Lloyd.  Quantum Support Vector Machine for Big Data Classification.  Physical review letters, 113(13):130503, 2014.
  * [20] Ahmed Shokry, Mostafa Mahmoud Elhamshary, and Moustafa Youssef.  DynamicSLAM: Leveraging Human Anchors for Ubiquitous Low-overhead Indoor Localization.  IEEE Transactions on Mobile Computing, 2020.
  * [21] Ahmed Shokry, Moustafa Elhamshary, and Moustafa Youssef.  The Tale of Two Localization Technologies: Enabling Accurate Low-overhead WiFi-based Localization for Low-end Phones.  ACM SIGSPATIAL, pages 1–10, 2017.
  * [22] Ahmed Shokry, Marwan Torki, and Moustafa Youssef.  DeepLoc: a Ubiquitous Accurate and Low-overhead Outdoor Cellular Localization System.  Proceedings of the 26th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems, pages 339–348, 2018.
  * [23] Ahmed Shokry and Moustafa Youssef.  Quantum Computing for Location Determination.  arXiv preprint arXiv:2106.11751, 2021.
  * [24] Ahmed Shokry and Moustafa Youssef.  Towards Quantum Computing for Location Tracking and Spatial Systems.  pages 278–281, 2021.
  * [25] Peter W Shor.  Polynomial-time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer.  SIAM review, 41(2):303–332, 1999.
  * [26] Andrei N. Soklakov and Rüdiger Schack.  Efficient State Preparation for a Register of Quantum Bits.  Phys. Rev. A, 73:012307, Jan 2006.
  * [27] Swamit S Tannu and Moinuddin K Qureshi.  A Case for Variability-aware Policies for Nisq-era Quantum Computers.  arXiv preprint arXiv:1805.10224, 2018.
  * [28] Qing Yang, Shijue Zheng, Ming Liu, and Yawen Zhang.  Research on Wi-Fi Indoor Positioning in a Smart Exhibition Hall based on Received Signal Strength Indication.  EURASIP Journal on Wireless Communications and Networking, 2019(1):275, Dec 2019.
  * [29] Moustafa Youssef.  Towards Truly Ubiquitous Indoor Localization on a Worldwide Scale.  Proceedings of the 23rd SIGSPATIAL International Conference on Advances in Geographic Information Systems, pages 1–4, 2015.
  * [30] Moustafa Youssef and Ashok Agrawala.  The Horus WLAN Location Determination System.  pages 205–218, 2005.
  * [31] Moustafa Youssef and Ashok Agrawala.  Location-clustering techniques for WLAN location determination systems.  International Journal of Computers and Applications, 28(3):278–284, 2006.
  * [32] Zhikuan Zhao, Jack K Fitzsimons, and Joseph F Fitzsimons.  Quantum-assisted Gaussian Process Regression.  Physical Review A, 99(5):052331, 2019.


# Finite volume based film flow and ice accretion models on aircraft wings

Tong Liu  Nuclear Power Institute of China, Chengdu 610213, PR China
Department of Fluid Mechanics, Northwestern Polytechnical University, Xi’an,
710072, PR China  Jinsheng Cai  Department of Fluid Mechanics, Northwestern
Polytechnical University, Xi’an, 710072, PR China  Kun Qu  Department of Fluid
Mechanics, Northwestern Polytechnical University, Xi’an, 710072, PR China
Shucheng Pan<EMAIL_ADDRESS>Department of Fluid Mechanics,
Northwestern Polytechnical University, Xi’an, 710072, PR China

###### Abstract

The thin runback water films driven by the gas flow, the pressure gradient and
the gravity on the iced aircraft surface are investigated in this paper. A
three-dimensional film flow model based on Finite Volume Method (FVM) and the
lubrication theory is proposed to describe the flow. The depth-averaged
velocity of the film is stored in Cartesian coordinates to avoid the
appearance of the metric tensors. The governing equations are discretized in
the first layer structured grid cell which is selected as the grids for film
flow. In order to verify this method, comparisons between numerical results
and experimental results of ice shapes on NACA 0012 airfoil and GLC-305 swept
wing are presented, both showing a good agreement for rime and glaze ice
condition. Overall, this model shows great potential to model ice accretion
reasonably under different icing conditions. Besides, the present method
doesn’t require analytic metric terms, and can be easily coupled to existing
finite volume solvers for logically Cartesian meshes.

††preprint: AIP/123-QED

##  I Introduction

The ice and water film may occur at the windward side of an aircraft when it
flies through clouds in which supercooled droplets are suspended. With the
driven force of the shear stress of the air, pressure gradient and gravity,
unfrozen water film may run back downward on the curved aircraft surface,
freeze, and form what called is glaze ice. Glaze ice is a typical glossy ice
shape with a single or double horns at the leading edge. Such severe ice
accretion would modify aircraft’s geometry, degenerate its aerodynamic
characteristics, and pose a serious threat to the safety of the aircraft
flight. Numerical simulation provides a low-cost way to study the ice
accretion and runback water film flow on the aircraft surface.

Typically, the film thickness on the iced surface is distinctly smaller than
its lateral dimension, and there are many other aspects dealing with the thin
water film evolution on curved three-dimensional surface, such as shallow
water problem, draining, coating and gravity driven film flow. Generally,
complex mathematical model and a mass of computational time are needed to
numerically solve the film flow. Lubrication (or long-wave, also long-scale)
theoryBenney (1966); Oron, Davis, and Bankoff (1997) is usually used to
address this problem. With respect to the full Navier–Stokes equations,
lubrication model simplifies the film flow problem reasonably and saves the
computational cost by solving only a single degenerate nonlinear equation for
the film thickness. To simulate the evolution of the thickness of a film on a
general curved substrate, RoyRoy, Roberts, and Simpson (2002) derived a
lubrication model expressed in terms of the film thickness and in a coordinate
system fitted to the curved substrate. The effects of the curvature of the
substrate, gravity and inertia are included to accurately describe the film
flow. Based on Roy’s work, RobertsRoberts and Li (2006) derived a more
comprehensive model of the dynamics of the film, and the model is expressed in
terms of the film thickness and the averaged lateral velocity. The model
resolved wide range of physical interactions between the various physical
process of inertia, surface tension, gravity and substrate curvature. Flow on
different substrate shapes, including flat, cylindrical, channel and
spherical, were simulated to illustrate its wide application. Jean-
LucThiffeault and Kamhawi (2006a) studied the steady gravity-driven flow of a
thin layer of viscous fluid over a curved substrate with topographical
variations. Different from Roy’s modelRoy, Roberts, and Simpson (2002), the
controlling equation is expressed in nonorthogonal coordinates. Besides, a
correction terms is introduced to the mass-conservation equation and vertical
velocity to ensure that the kinematic boundary condition at the free surface
is satisfied exactly. Using lubrication theory, HowellHowell (2003) derived
the general leading-order equations governing the flow of a thin liquid film
over a moving, curved substrate. Furthermore, the effects of the curvature of
the substrate are investigated, and three possible distinguished limits are
identified.

Aircraft ice accretion is a much more complex situation involving not only
film flow but also phase change, heat transfer and interface propagation.
Based on the Lubrication theoryBenney (1966) and Stefan problemHill (1987),
MyersMyers, Charpin, and Thompson (2002) first proposed the mathematical model
for ice accretion and water flow on the flat plate. The glaze icing rate is
determined by Stefan condition, and the film flow is driven by air shear,
gravity, pressure gradient and surface tension. For film flow on arbitrary
three-dimensional surface, a governing equation is derived based on the
curvilinear orthogonal coordinates, in which the first and second fundamental
metric tensors are used to carry information about the geometry of the curved
surfaceMyers, Charpin, and Chapman (2002); Myers and Charpin (2004). Cao and
Hou (2016) extended Myers’ model for curvilinear nonorthogonal grids systems.
As a result, besides the first and second fundamental metric tensors, second
cross fundamental forms of the surface were included to address the influence
of the metrics. Similarly, to simulate the ice accretion on helicopter rotors,
ChenChen  _et al._ (2018) developed a three-dimensional ice accretion model on
body-fitted nonorthogonal curvilinear coordinates. In this model, besides the
gravity force, the centrifugal caused by the rotation of the coordinate was
also accounted as body force.

However, film flow models aforementioned share the same disadvantage that
metric tensors should be evaluated to derive the governing equations on
curvilinear coordinate systems. The velocity field is expressed in curvilinear
coordinate systems, either. While in a general finite volume method flow
solver, the velocity field is usually expressed in Cartesian coordinate
system. Hence it would take a lot effort to implement the film flow module
with airflow solver. To combat this problem, CalhounCalhoun, Helzel, and
LeVeque (2008); Calhoun and Helzel (2009) proposed a finite volume method for
solving parabolic equations on logically Cartesian curved surface, and
commented that this method didn’t require analytic metric terms, showed second
order accuracy and can be easily coupled to existing finite volume solvers for
logically Cartesian meshes. Tukovi${\rm\acute{c}}$Tuković and Jasak (2012) and
RauterRauter  _et al._ (2018); Rauter and Tuković (2018) developed the finite
area method (FAM), and derived the controlling equation of a depth-integrated
shallow flow model for granular materials on three-dimensional mildly curved
topographies. The governing equatuons are expressed in three-dimensional
Cartesian coordinates, and share the advantage of easy coupling with three-
dimensional ambient flow equations.

In this paper, the governing equations for film flow and ice accretion process
on arbitrary curved surfaces are derived based on the finite area method, the
Lubrication theory, and the Stefan problem. The governing equations of film
flow are expressed in film thickness and depth-averaged velocity. The velocity
of the film flow is expressed in Cartesian coordinate system, therefore the
equations are easy to be coupled with the existing FVM based flow solver
ExstreamCai, Tsai, and Liu (2006); Xu, Liu, and Cai (2013); Xu  _et al._
(2014). Besides, the film flow flow is dominated by shear stress, pressure
gradient and gravity. The velocity profile in the film is approximated as a
polynomial function with respect to the film thickness, and the depth-averaged
velocity of the film can be expressed as a function of the film thickness and
the driving forces acting on the film.

The article is constructed as follows. Section II derives the film flow
governing equations based on the lubrication theory and finite volume method,
briefly introduces the ice accretion model, and presents the discretization
method. Section III is devoted to the simulated ice accretion results.
Finally, a conclusion will be given in section IV.

##  II Mathematical model

The problem considered hereby is shown schematically in Fig. 1. Ice
accumulates on the curved aircraft surface, and the unfrozen runback water
flows over the iced surface driven by the shear stress, pressure gradient and
gravity force. With reference to this figure, $B$ denotes the thickness of the
ice in the normal direction to the wall, and $h$ denotes the thickness of the
film. $T$ and $\theta$ are the temperature in the ice and water layers,
respectively. The substrate is denoted as $\mathscr{S}$, and the top surface
of the film is denoted as $\mathscr{F}$. With this article, our aim is to
develop a low-dimensional and easy established model for the film flow and
icing process on three-dimensional aircraft surface, and hence assumptions are
made as follows reasonably:

  * •

The film flow is incompressible, namely the density of the film flow is
assumed to be constant.

  * •

The aspect ratio and the Reynolds number of the thin layer film flow are
sufficiently small, which allows the lubrication theory to be used. The the
normal velocity component is negligible compared to tangential one, which
gives the constraint condition of the film velocity,
$\bm{u}\cdot\bm{n}_{wall}=0$ or $\bm{u}=\bm{u}_{t}$, where $\bm{n}_{wall}$ is
the unit normal vector of the wall pointing to the fluid, and $\bm{u}_{t}$ is
the tangential component of the film velocity.

  * •

The energy transfer across the film is driven by conduction rather than
advectionMyers, Charpin, and Chapman (2002).

Figure 1: Diagram of the two-dimensional ice and water layers on the
substrate.

###  II.1 Film flow model based on the finite volume methods

The control volume of the film flow above a patch of the substrate is depicted
in Fig. 2, which is extending across the film layer from $\zeta=0$ to
$\zeta=h$. Noting that the substrate refers to the ice surface rather than the
wall if ice accumulates on the wall. $\bm{n}_{b}$, $\bm{n}_{fs}$ and
$\bm{n}_{io}$ are the outward-pointing unit normal vectors of the faces of the
control volume, in which subscript fs, b and io denote the top surface, the
bottom surface, and the side surface, respectively. Because the thickness of
the film is very thin, the top surface is assumed to be parallel to the bottom
surface, namely $\bm{n}_{wall}=\bm{n}_{fs}=-\bm{n}_{b}$.

Figure 2: Sketch of the control volume of the film flow.

####  II.1.1 Mass conservation

To model the advection of the film flow on the aircraft surface, we employ
integral equations for mass balances written in the conservative forms

$\frac{\partial}{\partial t}\int_{V}\rho_{w}{\rm
d}V+\int_{V}\nabla\cdot\left[\rho_{w}(\bm{u}-\bm{u}_{S})\right]{\rm
d}V=\int_{S_{fs}}\dot{m}_{imp}{\rm d}S-\int_{S_{b}}\rho_{i}\frac{\partial
B}{\partial t}{\rm d}S$ (1)

Where $V$ is the arbitrary control volume for film, and $S$ denotes the
surfaces enclosing the control volume $V$, respectively. $\rho_{w}$,
$\rho_{i}$, $\bm{u}$, $\bm{u}_{S}$, $\dot{m}_{imp}$ and $B$ are bulk density
of water, bulk density of ice, film velocity, velocity of boundary surface,
incoming mass rate due to supercooled droplets impingement and ice thickness
respectively. The second term of right-hand-side of Eq. (1) denotes mass loss
due to ice accretion at the water-ice interface. Following the assumption that
the fluid is incompressible, the density of the film is constant. Applying the
Gauss Theorem and extracting the constant density we obtain

$\frac{\partial}{\partial t}\int_{V}{\rm
d}V+\oint_{S}(\bm{u}-\bm{u}_{S})\cdot\bm{n}{\rm
d}S=\int_{S_{fs}}\frac{\dot{m}_{imp}}{\rho_{w}}{\rm
d}S-\int_{S_{b}}\frac{\rho_{i}}{\rho_{w}}\frac{\partial B}{\partial t}{\rm
d}S$ (2)

where $\bm{n}$ is the outward unit normal vector of $S$. The first term of Eq.
(2) can be transformed into a surface-aligned curvilinear coordinate system:

$\int_{V}{\rm d}V=\int_{S_{b}}\int_{0}^{h}{\rm det}(\bm{J}){\rm d}\zeta{\rm
dS}$ (3)

where $\bm{J}$ is the Jacobian matrix induced by coordinate transformation
from Cartesian coordinates to curvilinear coordinates. According to
Ref.Thiffeault and Kamhawi (2006b), the determinant of matrix $\bm{J}$

${\rm det(\bm{J})}=1-\kappa\zeta+\mathcal{G}\zeta^{2}$ (4)

where $\kappa$ and $\mathcal{G}$ are the mean curvature and the Gaussian
curvature of the surface, respectively. $\kappa\zeta$ can be estimated as the
ratio of flow thickness to curvature radius. The film flow thickness during
ice accretion is on the scale of $1.0\times 10^{-5}{\rm m}$. Sharp convex and
concave are avoided during grid evolution and the maximum grid length along
$\xi$ and $\eta$ is usually great than $1.0\times 10^{-3}{\rm m}$. Therefore
$\kappa\zeta\ll 1$ and hence Eq. (3) can be written as

$\int_{V}{\rm d}V=\int_{S_{b}}\int_{0}^{h}{\rm d}\zeta{\rm
dS}=\int_{S_{b}}h{\rm dS}.$ (5)

Similarly, we could obtain

$\int_{S_{io}}{\rm d}S=\oint_{L_{io}}\int_{0}^{h}{\rm det}(\bm{J}){\rm
d}\zeta{\rm d}L\approx\oint_{L_{io}}\int_{0}^{h}{\rm d}\zeta{\rm d}L$ (6)

where $L_{io}$ is the side length of the bottom surface.

The second term of Eq. (2) can be split into the integral on each surface:

$\begin{split}&\oint_{S}(\bm{u}-\bm{u}_{S})\cdot\bm{n}{\rm
d}S=\int_{S_{b}}(\bm{u}-\bm{u}_{b})\cdot\bm{n}_{b}{\rm d}S\\\
&+\int_{S_{fs}}(\bm{u}-\bm{u}_{fs})\cdot\bm{n}_{fs}{\rm
d}S+\int_{S_{io}}(\bm{u}-\bm{u}_{io})\cdot\bm{n}_{io}{\rm d}S\end{split}$ (7)

The velocity of each surface of the control volume is negligible due to ice
grows slowly, namely $\bm{u}_{S}=\bm{u}_{b}=\bm{u}_{fs}=\bm{u}_{io}=\bm{0}$.
Besides, the assumptions mentioned before, $\bm{u}\cdot\bm{n}_{wall}=0$ and
$\bm{n}_{wall}=\bm{n}_{fs}=-\bm{n}_{b}$, indicate that the first and the
second terms of right-hand-side of Eq. (7) are 0. Combining with Eq. (6), the
third term of Eq. (2) cam be written as

$\int_{S_{io}}(\bm{u}-\bm{u}_{io})\cdot\bm{n}_{io}{\rm
d}S=\oint_{L_{io}}\int_{0}^{h}\bm{u}\cdot\bm{n}_{io}{\rm d}\zeta{\rm
d}L=\oint_{L_{io}}\left[\left(\int_{0}^{h}\bm{u}{\rm
d}\zeta\right)\cdot\bm{n}_{io}\right]{\rm d}L$ (8)

The depth-averaged velocity is defined as

$\bar{\bm{u}}=\frac{1}{h}\int_{0}^{h}\bm{u}{\rm d}\zeta$ (9)

and then Eq. (2) can be simplified as

$\oint_{S}(\bm{u}-\bm{u}_{S})\cdot\bm{n}{\rm
d}S=\oint_{L_{io}}h\bm{\bar{u}}\cdot\bm{n}_{io}{\rm d}L$ (10)

combining with Eq. (5), Eq. (10) and the relation between top surface and
bottom surface of the control volume, which reads $S_{fs}=S_{b}{{\rm
det}(\bm{J})}\approx S_{b}$, we get the mass conservation equation of the film
flow in conservation form

$\frac{\partial}{\partial t}\int_{S_{b}}h{\rm
d}S+\oint_{L_{io}}h\bm{\bar{u}}\cdot\bm{n}_{io}{\rm
d}L=\frac{1}{\rho_{w}}\int_{S_{b}}\dot{m}_{imp}{\rm
d}S-\frac{\rho_{i}}{\rho_{w}}\frac{\partial}{\partial t}\int_{S_{b}}B{\rm d}S$
(11)

The control volume $V$ shrinks to the control area as depicted in Fig. 3.
$\bm{n}_{e}$ is the outward normal vector of the side of the control area, and
is equal to $\bm{n}_{io}$. The equation above can be rewritten in the
differential form for any control area $S_{b}$Rauter and Tuković (2018) as

$\frac{\partial h}{\partial
t}+\nabla_{s}\cdot\left(h\bm{\bar{u}}\right)=\frac{\dot{m}_{imp}}{\rho_{w}}-\frac{\rho_{i}}{\rho_{w}}\frac{\partial
B}{\partial t}$ (12)

where $\nabla_{s}$ is the Nabla operator on the curves surface.

Figure 3: Sketch of the control area of the film flow.

####  II.1.2 Moment conservation

Before deriving the moment conservation equation for the film flow, we
introduce the boundary condition related to the film flow first. The fluid
immediately in contact with the bottom substrate $\mathscr{S}$ does not slip
along the stationary substrate, namely at the wall the no-slip boundary
condition reads $\left.\bm{u}\right|_{\zeta=0}=0$. At the film surface
$\mathscr{F}$, where $\eta=h$, forces acting on the film surface should be in
equilibrium. The tangential force balance at the air-water interface yields a
relation between normal derivative of tangential velocity on the two sides of
the water-air interfaceTuković and Jasak (2012):

$\left.\mu_{w}\nabla\bm{u}\cdot\bm{n}_{fs}\right|_{\zeta=h}=\left.\mu_{a}\nabla\bm{u}_{a,t}\cdot\bm{n}_{fs}\right|_{\zeta=h},$
(13)

where $\mu_{w}$ is the dynamic viscosity of water, and $\bm{u}_{a,t}$ is the
velocity component of air tangent to the wall. This equation is equivalent to

$\left.\mu_{w}\frac{\partial\bm{u}}{\partial
n}\right|_{\zeta=h}=\bm{\tau}_{a},$ (14)

where $\tau_{a}$ is the air shear stress acting on the film. The normal force
balance at the air-water interface gives

$\left.p\right|_{\zeta=h}=p_{fs}=p_{a}+\sigma_{film}\kappa$ (15)

where $p_{fs}$ is the pressure of the film at the top surface of the film,
$p_{a}$ is the pressure of the air at the air-water interface, and
$\sigma_{film}$ is the surface tension of the film. Surface tension has little
effect on the film velocity, therefore it is negligible during the simulation.

According to the assumptions mentioned before and the lubrication theory, the
film flow is treated as an incompressible laminar flow and the convection term
can be neglected, therefore the conservation form of the moment equation take
the form

$-\int_{S}p\bm{n}{\rm d}S+\int_{V}\rho_{w}\bm{g}{\rm
d}V+\int_{V}\mu_{w}\nabla^{2}\bm{u}{\rm d}V=0,$ (16)

where $p$ is the pressure in the film, and $\bm{g}$ is the gravitational
acceleration vector. Applying the approximation in the Appendix A, the moment
conservation equation can be expressed as integral on the control area

$\begin{split}&-\oint_{L_{io}}\bm{n}_{io}h\bar{p}{\rm
d}L-\int_{S_{b}}(p_{b}-p_{a})\bm{n}_{wall}{\rm d}S\\\
&+\int_{S_{b}}h\rho_{w}\bm{g}{\rm d}S+\int_{S_{b}}\bm{\tau}_{a}{\rm
d}S-\int_{S_{b}}\bm{\tau}_{b}{\rm d}S=0,\end{split}$ (17)

Rewritting this equation in differential form, we obtain

$-\nabla_{s}\left(h\bar{p}\right)-(p_{b}-p_{a})\bm{n}_{wall}+h\rho_{w}\bm{g}+(\bm{\tau}_{a}-\bm{\tau}_{b})=0$
(18)

Applying surface normal and surface tangential projection to Eq. (18) yields

$p_{b}=p_{a}+h\rho_{w}g_{n}$ (19)

and

$-\nabla_{s}\left(h\bar{p}\right)+h\rho_{w}\bm{g}_{s}+(\bm{\tau}_{a}-\bm{\tau}_{b})=0,$
(20)

where $\bm{g}_{n}=\bm{n}_{wall}\cdot\left(\bm{n}_{wall}\cdot\bm{g}\right)$ and
$\bm{g}_{s}=\left(\bm{I}-\bm{n}_{wall}\bigotimes\bm{n}_{wall}\right)\cdot\bm{g}$.
$\bar{p}$ is the depth-averaged pressure of the film, which reads

$\bar{p}=\dfrac{p_{a}+p_{b}}{2}=p_{a}+\frac{1}{2}h\rho_{w}g_{n}.$ (21)

The velocity profile in the water film usually follows a linearBourgault,
Beaugendre, and Habashi (2000) or polynomial functionMyers (2001); Myers,
Charpin, and Chapman (2002); Wang and Rothmayer (2009); Cao and Hou (2016);
Chen  _et al._ (2018) along the film thickness $\zeta$, and the latter leads
to

$\bm{u}(\zeta)=\bm{a}\zeta^{2}+\bm{b}\zeta+\bm{c},$ (22)

where $\bm{a}$, $\bm{b}$ and $\bm{c}$ are the coefficient vectors. Applying
the no-slip boundary condition at the water-substrate interface and the shear
stress boundary condition (14) and integrating the velocity along the film
thickness finally yields

$\bar{\bm{u}}=\frac{\left(\bm{\tau}_{a}+2\bm{\tau}_{b}\right)h}{6\mu_{w}}$
(23)

The vectors $\bm{\tau}_{a}$ and $\bm{\tau}_{b}$ are parallel to the substrate,
which makes the film velocity satisfy the constraint condition
$\bm{u}\cdot\bm{n}_{b}=0$ automatically. Moreover, combining Eq. (20) and Eq.
(23) gives another form of the depth-averaged velocity:

$\bar{\bm{u}}=\frac{h}{2\mu_{w}}\bm{\tau}_{a}+\frac{h}{3\mu_{w}}\left[-\nabla_{s}\left(h\bar{p}\right)+h\rho_{w}\bm{g}_{s}\right]$
(24)

It indicates that air shear play a key role in the driving forces of the film
flow. Furthermore, the film flux in Eq. (12) can be expressed as

$\bm{F}=h\bar{\bm{u}}=\frac{h^{2}}{2\mu_{w}}\bm{\tau}_{a}+\frac{h^{2}}{3\mu_{w}}\left[-\nabla_{s}\left(h\bar{p}\right)+h\rho_{w}\bm{g}_{s}\right]$
(25)

####  II.1.3 Energy conservation

The temperature at the ice-water interface is freezing temperature of water
$T_{f}$, and the film is very thin. Therefore the depth-averaged temperature
of the film is close to $T_{f}$, hence the assumption is made that there is no
heat exchange between adjacent control areas along the wall surface. Applying
the lubrication theory and the approximation aforementioned in a similar way,
the energy conservation equation for the film flow on the arbitrary curved
surface can be reduced toMyers, Charpin, and Thompson (2002); Myers, Charpin,
and Chapman (2002)

$\int_{V}\nabla^{2}\theta{\rm d}V=0.$ (26)

The equation above indicates that the temperature gradient in the film normal
to the wall is constant, i.e., the temperature profile in the film in linear.

###  II.2 Rime and glaze ice accretion model

Under rime ice condition, all impinging droplets freeze and there is no film
flow on the surface, hence all terms on the left-hand side of Eq. (11) are
identically zero. The icing rate is simply proportional to the droplet mass
rate impinging on the control area, and the control equation can be written in
differential form as

$\frac{\partial B}{\partial t}=\frac{\dot{m}_{imp}}{\rho_{i}}.$ (27)

Under glaze ice condition instead, by applying the energy balance on the ice-
water interfaceMyers, Charpin, and Thompson (2002); Myers, Charpin, and
Chapman (2002), the icing rate in differential form can be derived as

$\rho_{i}L_{f}\frac{\partial B}{\partial t}=k_{i}\left.\frac{\partial
T}{\partial n}\right|_{\zeta=0}-k_{w}\left.\frac{\partial\theta}{\partial
n}\right|_{\zeta=0},$ (28)

where $L_{f}$, $k_{i}$ and $k_{w}$ are the latent heat of ice accretion, the
thermal conductivities of ice and water, respectively. Eq. (28) implies that
the energy released during the solidification of impinging water is conducted
away through the ice and water layers, which is also known as the Stefan
condition or phase change conditionHill (1987). The temperature distribution
in the ice layer is assumed to be linearMyers, Charpin, and Chapman (2002);
Myers, Charpin, and Thompson (2002). Under this assumption, the energy balance
(28) can be expressed in its final form

$\rho_{i}L_{f}\frac{\partial B}{\partial
t}=k_{i}\frac{T_{f}-T_{s}}{B}-k_{w}\frac{Q_{glaze}^{+}-q_{glaze}^{-}\left(T_{f}-T_{a}\right)}{k_{w}+hq_{glaze}^{-}},$
(29)

where $Q_{glaze}^{+}$ and $q_{glaze}^{-}$ are the heat energy gained and the
rate of the heat energy lost during ice accretion, respectively. Details of
these two terms are referred to Liu  _et al._ (2019).

###  II.3 Discretization and solution methods

The ice accretion on the aircraft surface is governed by Eq. (12), Eq. (20),
Eq. (27) and Eq. (29). The ice thickness $B$, the film thickness $h$ and the
depth-averaged velocity of the film $\bar{\bm{u}}$ are unknown, and now we
seek the discretization and solution methods of the governing equations to
solve this problem. Under rime ice condition (dry accretion), the ice
thickness increment is solved according to the discretized form of Eq. (27)

$\Delta B^{n}=\frac{\dot{m}_{imp}\Delta t}{\rho_{i}},$ (30)

where $\Delta t$ is the time step employed in the calculation. Under glaze ice
condition (wet accretion), the discretized form of Eq. (29) gives

$\Delta B^{n}=\frac{\Delta
t}{\rho_{i}L_{F}}\left[k_{i}\frac{T_{f}-T_{s}}{B^{n}}-k_{w}\frac{Q_{glaze}^{+}-q_{glaze}^{-}\left(T_{f}-T_{a}\right)}{k_{w}+h^{n}q_{glaze}^{-}}\right].$
(31)

The pressure gradient of the film should be solved first to evaluate the
depth-averaged velocity and the flux of the film flow explicitly, and applying
the Gauss Theorem yields

$\nabla_{s}\left(h\bar{p}\right)=\frac{\sum_{e}\bm{n}_{e}h_{e}^{n}\bar{p}_{e}L_{e}}{S_{b}}$
(32)

where $L_{e}$ is the edge length and $\bm{n}_{e}$ is the unit normal vectors
on lateral surface. A high-quality body-fitted structured grids with
orthogonality of grids ensured near the wall are generated hence $\bm{n}_{e}$
is equal to $\bm{n}_{io}$ approximatively. The mass conservation equation (11)
can be discretized as follows to obtain the film thickness increment

$\Delta S\frac{\Delta h^{n}}{\Delta
t}+\sum_{e}\bm{F}^{n+1}\cdot\bm{n}_{e}L_{e}=\Delta
S\left(\frac{\dot{m}_{imp}}{\rho_{w}}-\frac{\rho_{i}}{\rho_{w}}\frac{\Delta
B^{n}}{\Delta t}\right)$ (33)

The discretization of the convective term is performed with the first-order
upwind scheme to smooth the oscillations of the film thickness and the Line
Successive Over Relaxation (LSOR) is used to solve the equations on the iced
surface.

A solution strategy presented by Myers, Charpin, and Chapman (2002) is adopted
to solve the ice accretion problem. The type of the ice accretion at each
control area is assumed to be wet accretion therefore the ice and film
thickness are determined by Eq. (31) and Eq. (33), which requires that
$B^{n}\neq 0$. Hence a precursor ice thickness $B_{p}=1.0\times
10^{-12}\mathrm{m}$ is specified initially to avoid the problem mentioned
above. Besides, a precursor film $h_{p}$ is introduced to avoid the difficulty
associated with the advancing contact line. In this paper, the aircraft
surface is set to be covered with a thin film whose thickness is $1.0\times
10^{-12}\mathrm{m}$ initially. The wet accretion assumption holds in the
control area if the new film thickness $h^{n+1}$ is greater than the precursor
film thickness $h_{p}$. While the wet assumption is invalid if the new film
thickness is smaller than the precursor film thickness. In this case, the new
ice thickness increment is calculated by Eq. (30) and the film thickness is
set to $h_{p}$.

###  II.4 Numerical methods

As shown in Fig. 4, Current ice accretion code is composed of several modules,
i.e., a flow solver, a droplets flow solver, an ice accretion module with
thermodynamic and film flow model, and a mesh module. The air flow solver
called ExstreamXu, Liu, and Cai (2013); Xu  _et al._ (2014); Cai, Tsai, and
Liu (2006) solves the Reynolds-averaged Navier–Stokes (RANS) equations to
evaluate the driving force acting on the film. The influence of the surface
roughness on the convective heat transfer characteristics is taken into
account by extending the turbulence modelKnopp, Eisfeld, and Calvo (2009). The
droplet flow solver based on the Eulerian method is utilized to obtain the
water collection efficiencies on the wall. The mesh module adopts hierarchical
overset grid strategyCai, Tsai, and Liu (2006) and parabolic grid generation
approach to generate high-quality curvilinear body-fitted structured grids on
the clean and complex iced geometry. A quasi-steady multi-step
algorithmVerdin, Charpin, and Thompson (2009) is employed to simulate the
unsteady icing process. During the unsteady icing process in each single step,
the geometry of the iced airfoil changes slowly, implying that the unsteady
effects of the air flow can be neglected. Thus, the parameters of the air flow
and droplet flow are kept as constant during the simulation of ice accretion
and film flow. For further information on these modules mentioned above, the
reader is referred to Liu  _et al._ (2019).

Figure 4: Flowchart of ice accretion.

##  III Numerical results and Analysis

The verification of the film flow and ice accretion models is verified on both
two and three dimensional cases where numerical results obtained by current
method are compared with available experimental and numerical results. In
addition, the film thickness distribution under different icing conditions are
studied.

###  III.1 Film flow on NACA 0012 airfoil

A simulation performed on NACA0012 airfoil presented by Lavoie (2017) is
carried out to verify current film flow model. The test parameters are the
same with that of Run 308 in Table 1, except that the angle of attack is 0∘,
the temperature is 288.15$\mathrm{K}$ and the exposure time is 4$\mathrm{s}$.

Results of the numerical solution are presented in Fig. 5 for the film
thickness evolution, the film depth-averaged velocity evolution, the water
collection efficiency and the air shear stress. The colored solid lines, which
are the time series of streamwise profiles plotted every 0.5 $\mathrm{s}$,
present the distribution of film thickness and the film velocity, while the
black dashed lines present the distribution of the water collection efficiency
and the air shear stress acting on the film. Velocity and air shear stress in
the counter-clockwise direction along the airfoil surface is taken as negative
values. Note that the two subfigures share the same X axis. As we can see, the
impinging limits of the droplets locates around the stagnation point and its
range is 6.6% of the surface arc length. Since the surrounding temperature is
warm enough, impinging droplets do not freeze at all and flow all the way from
the impinging area to the unimpinging area under the influence of the driving
force. The upper subfigure in Fig. 5 shows that the film cover about 70% of
the airfoil surface at $\rm{t=4s}$. The lower subfigure shows that a sharp
increase of the film velocity is observed near the stagnation point, which
indicates that water around the stagnation point is flow downstream rapidly.
Hence the film thickness around the stagnation point is smaller than that at
downstream distinctly even though the water collection efficiency reaches its
maximum value at the stagnation point. Eq. (24) states that the air shear
dominates the depth-averaged velocity of the film, thus a decrease of the air
shear stress leads to the decrease of the film velocity. As a result, the film
thickness at the advancing front increases over time gradually.

Figure 5: The film flow thickness, the film flow depth-averaged velocity, the
water collection efficiency and the air shear stress acting on the film
against the non-dimensional curvilinear distances computed clockwise from the
trailing edge. The stagnation point locates at $\mathrm{s/c}=0$.

###  III.2 Ice accretion on NACA 0012 airfoil

Wright presented a series of icing results for different icing conditions on
different airfoilsWright and Rutkowski (1999). Experimental results of NASA
Lewis Icing Research Tunnel (IRT) and numerical results of LEWICE 2.0 are
provided for each case. In this paper, 9 cases on NACA 0012 airfoil, including
rime and glaze ice accretion, are selected to validate current model
thoroughly. The corresponding parameters for these cases are presented in
Table 1. The chord length of the airfoil is 0.5334${\rm m}$, and the angle of
attack used in these cases is 3.5∘.

Table 1: Ice conditions for NACA0012 airfoil. Case |  $U_{\infty}$, $\rm m/s$ |  $T_{\infty}$, $\rm K$ | LWC, $\rm g/m^{3}$ | MVD, $\rm\mu m$ | Time, $\rm min$
---|---|---|---|---|---
Run 308 | 102.8 | 262.04 | 1.0 | 20 | 3.85
Run 316 | 102.8 | 262.04 | 0.55 | 20 | 3.22
Run 401 | 102.8 | 265.37 | 0.55 | 20 | 7
Run 403 | 102.8 | 262.04 | 0.55 | 20 | 7
Run 405 | 102.8 | 250.37 | 0.55 | 20 | 7
Run 409 | 67.1 | 265.07 | 1.3 | 30 | 6
Run 421 | 67.1 | 268.40 | 1.0 | 20 | 6
Run 422 | 67.1 | 266.74 | 1.0 | 20 | 6
Run 423 | 67.1 | 265.07 | 1.0 | 20 | 6

Fig. 6 illustrates the comparisons between current ice shapes and ice shapes
of IRT experiment and LEWICE 2.0. Note that figures in the same column share
the same X axis, while figures in the same row share the same Y axis. Run 401,
run 403 and run 405 (2nd column of Fig. 6) share the same icing conditions
except for the temperature, which decreases from 265.37$\rm K$ to 250.37$\rm
K$, and the ice shapes present a transition between glaze ice and rime ice.
For run 405, all impinging droplets freeze, thus typical rime ice and smooth
streamwise shape form. For run 403, temperature increases to 262.04$\rm K$.
Unfrozen droplets flow downstream and horns form at the upper and lower
surfaces near the trailing edge. With temperature increasing to 265.37$\rm K$,
unfrozen droplets were transported downstream further. Therefore the upper
horn shifts to downstream, while the lower horn disappear, resulting in a
wider range of ice layer. It is clearly evident that with the increase of the
temperature, the ice layer at the stagnation point grows thicker and the
position of the ice horns shifts more downstream. Run 423, run 422 and run 421
(3rd column of Fig. 6) also share the same icing conditions except for the
temperature, which increases from 265.07$\rm K$ to 267.37$\rm K$, and all the
ice shapes are typical glaze ice. The impinging limits of these cases are the
same, because the trajectories of the droplets are barely affected by the
temperature. As mentioned above, with the increase of the temperature, more
film flow downstream, and ice tends to spread thinner and farther along the
lower surface.

Overall, good agreement among predicted results by current model, experimental
results and that of LEWICE is obtained. For few cases, a small discrepancy
between predicted and experimental ice shape is observed at the upper horn.
Main reason may lie in the imprecise prediction of air shear and the
convective heat transfer coefficient, and further work will be required to
solve this problem.

Figure 6: Ice accretion on NACA 0012 airfoil.

###  III.3 Ice accretion on three-dimensional surface

The experimental results on the GLC-305 swept wing presented by
PapadakisPapadakis  _et al._ (2003) have been widely used to validate ice
accretion models. The wing section in the streamwise direction is the GLC-305
airfoil and is constant from the root to the tip, and there is a geometry
twist of -4∘ from wing root to wing tip. More details of the swept wing
mentioned above refers to Papadakis  _et al._ (2003). In this paper, two glaze
ice cases, CS10 and IS10, are selected to validate current model, and
numerical results are compared to experiment and that of LEWICE. Section
comparisons are made at three spanwise locations A, B and C. The section at
station C is located at the wing root and in the streamwise direction. The
sections at stations B and A are taken normal to the wing leading edge at
36.8% and 73.6% semispan, respectively. The corresponding icing conditions are
listed in Table 2. Similar to the LEWICE3D, the computations are performed
using the single step algorithm.

Table 2: Icing conditions for GLC-305 swept wing. Parameters | CS10 | IS10
---|---|---
MAC, $\rm m$ | 0.4755 | 0.4755
$U_{\infty}$, $\rm m/s$ | 89.99 | 67.06
AOA, deg | 4 | 4
$T_{\infty}$, $\rm K$ | 261.87 | 269.26
$P_{\infty}$, $\rm Pa$ | 101300 | 101300
LWC, $\rm g/m^{3}$ | 0.68 | 0.65
MVD, $\rm\mu m$ | 20 | 20
Time, $\rm min$ | 10 | 10

Fig. 7 presents the sectional ice accretion shapes of case CS10 at sections A,
B and C. Overall, good agreement in ice accretion predictions is observed for
three sections between the numerical results and the experimental results. The
ice thickness around the stagnation point and the ice orientation for all the
three sections are in accordance with the experimental data, whereas LEWICE
trends to underpredict the ice thickness at the stagnation point for the lack
of consideration of the heat conduction through the ice and water layers. The
ice growth direction and the ice thickness at the upper and lower horns are
well captured at section A (Fig. 7(a)) and section B (Fig. 7(b)). Among the
three sections, section C (Fig. 7(c)) presents the poorest agreement with the
experimental data both for the icing limits and the ice thickness in
downstream regions of the stagnation point. Possible explanations for such
differences may be attributed to the fact that in the experiment a fuselage-
like body is used to mount the swept wing in the wind tunnel and away from the
near-wall flow, while in the numerical simulation a symmetry boundary
condition is imposed at the wing rootPapadakis  _et al._ (2003); Cao and Hou
(2016); Liu  _et al._ (2019). Different treatments leads to different flow
characteristic near the wing root, which results in the discrepancy in the
water collection efficiency and the ice shapes. Another possible reason may
lie in the fact that the ice shape is computed in a single-step method.
Results show that compared to multi-step calculations, single-step
calculations significantly overpredicts the ice thickness around the icing
limitsVerdin, Charpin, and Thompson (2009). Fig. 8 presents the sectional ice
accretion shapes of case IS10 at sections A, B and C, and the predicted
results compared favorably with the experimental results. Compared to current
model, LEWICE still underestimates the ice thickness at the stagnation point.
While one significant difference is observed between predicted results and
experiment results. At section C, the ice thickness around the impinging
limits are much thicker than that of experimental results, which is similar
with case CS10.

(a) Section A.

(b) Section B.

(c) Section C.

Figure 7: Ice shapes comparisons of CS10 for GLC-305 swept wing at sections A,
B and C.

(a) Section A.

(b) Section B.

(c) Section C.

Figure 8: Ice shapes comparisons of IS10 for GLC-305 swept wing at sections A,
B and C.

Fig. 9 shows the maximum water collection efficiency $\beta$ and convective
heat transfer coefficient HTC in the streamwise direction across all the
semispan positions, while Fig. 10 and Fig. 11 illustrate the evolution of the
maximum film and ice thickness in the streamwise direction across all the
semispan positions. Note that the spanwise position is normalized by the
semispan length. CS10 and IS10 are typical glaze ice case, therefore similar
trends of $\beta$, HTC, film thickness and ice thickness are observed for
these two cases. Compared to CS10, less flow velocity of IS10 leads to, on the
whole, less water collection efficiency, and finally less ice accretion
volume.

As shown in Fig. 9, $\beta$ and HTC increase from the wing root to the wing
tip, which indicates that the wing tip tends to collect more supercooled water
and the supercooled water are more likely to freeze near the wing tip. This
leads to the decrease of the film thickness in the outward spanwise direction,
as shown in Fig. 10 and Fig. 12. As we can see in these two figures, film
first appears at the wing root and near the wing tip after about
40$\mathrm{s}$. Then the film spreads downstream, and the overall film
thickness increases gradually over time, which means that more unfrozen water
run back over the wing surface and freeze in downstream regions. After about
300$\mathrm{s}$, the overall film thickness trends to reach some constants,
and the film flow becomes stable, which indicates that the mass of water
entering the control cell is nearly equal to that leaving the control cell.
Furthermore, film thickness at the wing root is significantly greater than
that at any other regions because of the smaller HTC at the wing root. As a
result, the predicted ice horn at section C is more distinct than that at
section B and section A.

(a) Case CS10.

(b) Case IS10.

Figure 9: Water collection efficiency and convective heat transfer coefficient
comparisons across the span for case CS10 and case IS10.

(a) Case CS10.

(b) Case IS10.

Figure 10: Film thickness comparison across the span for case CS10 and case
IS10.

(a) Case CS10.

(b) Case IS10.

Figure 11: Ice thickness comparison across the span for case CS10 and case
IS10.

(a) Case CS10.

(b) Case IS10.

Figure 12: Film thickness evolution for case CS10 and case IS10.

##  IV Conclusion

Based on the finite area method, we derive the governing equations for film
flow and ice accretion process on arbitrary curved surfaces. The governing
equations are expressed in Cartesian coordinate system therefore it is easy to
be coupled to existing FVM flow solver Exstream. The film flow flow is
dominated by shear stress, pressure gradient and gravity. The velocity field
is expressed in Cartesian coordinate system. Besides, the velocity profile in
the film is approximated as a polynomial function with respect to the film
thickness and finally the depth-averaged velocity of the film can be expressed
as a function of the film thickness and the shear stress acting on the film.
Numerical simulations are conducted on airfoils and the wing, and current film
and ice accretion model are verified to be reasonable by calculating the flow
behavior of the runback water and the shape feature of the ice shapes.
Numerical results show that, for some glaze ice conditions, as temperature
increases, ice tends to extend higher on the upper surface and spread thinner
and farther along the lower surface.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(No. 11902271 and No. 91952203), the Fundamental Research Funds for the
Central Universities of China (No. G2020KY05101), the Foundation of National
Key Laboratory (No. 6142201190303), and the 111 Project of China (B17037). The
authors acknowledge the computational resources provided by the Tianhe-2
supercomputer of the National Supercomputer Center in Guangzhou (NSCC-GZ).

##  Appendix A Moment conservation

The first term of the left-hand-side of Eq. (16) can be decomposed as

$\int_{S}p\bm{n}{\rm d}S=\int_{S_{io}}p\bm{n}_{io}{\rm
d}S+\int_{S_{b}}p_{b}\bm{n}_{b}{\rm d}S+\int_{S_{fs}}p_{fs}\bm{n}_{fs}{\rm
d}S$ (34)

where $p_{b}$ and $p_{fs}$ are the pressure at the bottom and the surface of
the film respectively. Here a simple linear pressure profile in the film layer
is assumed, then the pressure along the film thickness is
$p(\zeta)=p_{b}+\dfrac{\zeta}{h}(p_{fs}-p_{b})$. Applying the depth-
integration to the first term of the right-hand-side of Eq. (34) we get

$\begin{split}\int_{S_{io}}p\bm{n}_{io}{\rm
d}S&=\oint_{L_{io}}\bm{n}_{io}\int_{0}^{h}p{\rm d}\zeta{\rm d}L\\\
&=\oint_{L_{io}}\bm{n}_{io}h\bar{p}{\rm d}L\end{split}$ (35)

where the depth-averaged pressure
$\bar{p}=\dfrac{1}{2}\left(p_{fs}+p_{b}\right)$. At last Eq. (34) gives

$\int_{S}p\bm{n}{\rm d}S=\oint_{L_{io}}\bm{n}_{io}h\bar{p}{\rm
d}L+\int_{S_{b}}(p_{b}-p_{a})\bm{n}_{b}{\rm d}S$ (36)

For the second term of the left-hand-side of Eq. (16), by performing the
depth-integration the body force induced by the gravity acceleration force can
be easily obtained as

$\int_{V}\rho_{w}\bm{g}{\rm d}V=\int_{S_{b}}h\rho_{w}\bm{g}{\rm d}S$ (37)

The third term of the left-hand-side of Eq. (16) can be split into

$\begin{split}\int_{V}\mu_{w}\nabla^{2}\bm{u}{\rm
d}V&=\oint_{S}\mu_{w}\nabla\bm{u}\cdot\bm{n}{\rm d}S\\\
&=\int_{S_{io}}\mu_{w}\nabla\bm{u}\cdot\bm{n}_{io}{\rm
d}S+\int_{S_{b}}\mu_{w}\nabla\bm{u}_{b}\cdot\bm{n}_{b}{\rm d}S\\\
&+\int_{S_{fs}}\mu_{w}\nabla\bm{u}_{fs}\cdot\bm{n}_{fs}{\rm d}S\end{split}$
(38)

The first term of the right-hand-side of Eq. (38) can be written as

$\int_{S_{io}}\mu_{w}\nabla\bm{u}\cdot\bm{n}_{io}{\rm
d}S=\oint_{L_{io}}\mu_{w}\int_{0}^{h}\nabla\bm{u}\cdot\bm{n}_{io}{\rm
d}\zeta{\rm d}L$ (39)

According to the shear continuity at the water-air interface,
$\mu_{w}\nabla\bm{u}_{fs}\cdot\bm{n}_{b}=\bm{\tau}_{a}$, the third term can be
written as

$\begin{split}\int_{S_{fs}}\mu_{w}\nabla\bm{u}_{fs}\cdot\bm{n}_{fs}{\rm
d}S&\approx-\int_{S_{b}}\mu_{w}\nabla\bm{u}_{fs}\cdot\bm{n}_{b}{\rm d}S\\\
&=-\int_{S_{b}}\bm{\tau}_{a}{\rm d}S\end{split}$ (40)

At the water-substrate/ice interface, we let
$\mu_{w}\nabla\bm{u}_{b}\cdot\bm{n}_{b}=\bm{\tau}_{b}$. Finally Eq. (38) is

$\begin{split}\int_{V}\mu_{w}\nabla^{2}\bm{u}{\rm
d}V&=\oint_{S}\mu_{w}\nabla\bm{u}\cdot\bm{n}{\rm d}S\\\
&=\oint_{L_{io}}\mu_{w}\int_{0}^{h}\nabla\bm{u}\cdot\bm{n}_{io}{\rm
d}\zeta{\rm d}L\\\ &+\int_{S_{b}}\bm{\tau}_{b}{\rm
d}S-\int_{S_{b}}\bm{\tau}_{a}{\rm d}S\end{split}$ (41)

Finally, substituting Eq. (36), Eq. (37) and Eq. (41) into Eq. (16) we obtain
the moment conservation equation Eq. (17).

To solve the pressure at the bottom surface, we multiply Eq. (18) with the
normal vector $\bm{n}_{b}$:

$\begin{split}&-\bm{n}_{b}\cdot\hat{\nabla}\left(h\bar{p}\right)-(p_{b}-p_{a})\bm{n}_{b}\cdot\bm{n}_{b}\\\
&+h\rho_{w}\bm{n}_{b}\cdot\bm{g}+\bm{n}_{b}\cdot(\bm{\tau}_{b}-\bm{\tau}_{a})=0\end{split}$
(42)

The first and the last terms at the left-hand-side of Eq. (42) are
approximately 0, and then the equation along the normal direction can be
deduced to

$-(p_{b}-p_{a})+h\rho_{w}g_{n}=0$ (43)

and thus we obtain

$p_{b}=p_{a}+h\rho_{w}g_{n}$ (44)

and further

$\bar{p}=p_{a}+\frac{1}{2}h\rho_{w}g_{n}$ (45)

To evaluate the surface tangential momentum equation, we multiply Eq. (43)
with $\bm{n}_{b}$ again and it yields

$-(p_{b}-p_{a})\bm{n}_{b}+h\rho_{w}\bm{g}_{n}=0$ (46)

where $\bm{g}_{n}=\bm{n}_{b}\cdot\left(\bm{n}_{b}\cdot\bm{g}\right)$.
Subtracting the equation above with Eq. (18) yields

$-\hat{\nabla}_{s}\left(h\bar{p}\right)+h\rho_{w}\bm{g}_{s}+(\bm{\tau}_{b}-\bm{\tau}_{a})=0$
(47)

The depth-averaged velocity of the film is

$\bar{\bm{u}}=\frac{\bm{a}}{3}h^{2}+\frac{\bm{b}}{2}h+\bm{c}$ (48)

While the gradient of velocity is

$\frac{\partial\bm{u}}{\partial\zeta}=2\bm{a}\zeta+\bm{b}.$ (49)

The no-slip boundary condition at the water-substrate interface gives
$\bm{c}=0$ directly. Applying the shear stress boundary condition (14) yields

$\bm{a}=\frac{\bm{\tau}_{a}-\bm{\tau}_{b}}{2h\mu_{w}}$ (50)

While according to the defination of $\bm{\tau}_{b}$ we can easily obtain

$\bm{b}=\frac{\bm{\tau}_{b}}{\mu_{w}}$ (51)

At last, substituting the coefficients above into Eq. (48) yields

$\bar{\bm{u}}=\frac{\left(\bm{\tau}_{a}+2\bm{\tau}_{b}\right)h}{6\mu_{w}}$
(52)

## References

  * Benney (1966) D. J. Benney, “Long waves on liquid films,” Journal of Mathematics and Physics 45, 150–155 (1966).
  * Oron, Davis, and Bankoff (1997) A. Oron, S. H. Davis, and S. G. Bankoff, “Long-scale evolution of thin liquid films,” Reviews of Modern Physics 69, 931–980 (1997).
  * Roy, Roberts, and Simpson (2002) R. V. Roy, A. J. Roberts, and M. E. Simpson, “A lubrication model of coating flows over a curved substrate in space,” Journal of Fluid Mechanics 454, 235–261 (2002).
  * Roberts and Li (2006) A. J. Roberts and Z. Li, “An accurate and comprehensive model of thin fluid flows with inertia on curved substrates,” Journal of Fluid Mechanics 553, 33–73 (2006).
  * Thiffeault and Kamhawi (2006a) J.-L. Thiffeault and K. Kamhawi, “Transport in thin gravity-driven flow over a curved substrate,” (2006a), arXiv:nlin/0607075 .
  * Howell (2003) P. Howell, “Surface-tension-driven flow on a moving curved surface,” Journal of Engineering Mathematics 45, 283–308 (2003).
  * Hill (1987) J. M. Hill,  _One-dimensional Stefan problems: an introduction_ , Vol. 31 (Longman Sc & Tech, 1987).
  * Myers, Charpin, and Thompson (2002) T. G. Myers, J. P. F. Charpin, and C. P. Thompson, “Slowly accreting ice due to supercooled water impacting on a cold surface,” Physics of Fluids 14, 240–256 (2002).
  * Myers, Charpin, and Chapman (2002) T. G. Myers, J. P. F. Charpin, and S. J. Chapman, “The flow and solidification of a thin fluid film on an arbitrary three-dimensional surface,” Physics of Fluids 14, 2788–2803 (2002).
  * Myers and Charpin (2004) T. Myers and J. Charpin, “A mathematical model for atmospheric ice accretion and water flow on a cold surface,” International Journal of Heat and Mass Transfer 47, 5483–5500 (2004).
  * Cao and Hou (2016) Y. Cao and S. Hou, “Extension to the myers model for calculation of three-dimensional glaze icing,” Journal of Aircraft 53, 106–116 (2016).
  * Chen  _et al._ (2018) N. Chen, H. Ji, G. Cao, and Y. Hu, “A three-dimensional mathematical model for simulating ice accretion on helicopter rotors,” Physics of Fluids 30, 083602 (2018).
  * Calhoun, Helzel, and LeVeque (2008) D. Calhoun, C. Helzel, and R. LeVeque, “Logically Rectangular Grids and Finite Volume Methods for PDEs in Circular and Spherical Domains,” SIAM Review 50, 723–752 (2008).
  * Calhoun and Helzel (2009) D. Calhoun and C. Helzel, “A Finite Volume Method for Solving Parabolic Equations on Logically Cartesian Curved Surface Meshes,” SIAM Journal on Scientific Computing 31, 4066–4099 (2009).
  * Tuković and Jasak (2012) Ž. Tuković and H. Jasak, “A moving mesh finite volume interface tracking method for surface tension dominated interfacial fluid flow,” Computers & Fluids 55, 70–84 (2012).
  * Rauter  _et al._ (2018) M. Rauter, A. Kofler, A. Huber, and W. Fellin, “faSavageHutterFOAM 1.0: depth-integrated simulation of dense snow avalanches on natural terrain with OpenFOAM,” Geoscientific Model Development 11, 2923–2939 (2018).
  * Rauter and Tuković (2018) M. Rauter and Ž. Tuković, “A finite area scheme for shallow granular flows on three-dimensional surfaces,” Computers & Fluids 166, 184–199 (2018).
  * Cai, Tsai, and Liu (2006) J. Cai, H. M. Tsai, and F. Liu, “A parallel viscous flow solver on multi-block overset grids,” Computers & Fluids 35, 1290–1301 (2006).
  * Xu, Liu, and Cai (2013) J. Xu, Q. Liu, and J. Cai, “Numerical Simulations for DLR-F6 Wing/Body/Nacelle/Pylon with Enhanced Implicit Hole Cutting Method,” in  _Parallel Computational Fluid Dynamics_, Communications in Computer and Information Science (Springer, Berlin, Heidelberg, 2013) pp. 185–194.
  * Xu  _et al._ (2014) J. Xu, J. Cai, Q. Liu, and K. Qu, “Flow Simulations by Enhanced Implicit-Hole-Cutting Method on Overset Grids,” Journal of Aircraft 51, 1401–1409 (2014).
  * Thiffeault and Kamhawi (2006b) J.-L. Thiffeault and K. Kamhawi, “Transport in thin gravity-driven flow over a curved substrate,” arXiv preprint nlin/0607075 (2006b).
  * Bourgault, Beaugendre, and Habashi (2000) Y. Bourgault, H. Beaugendre, and W. G. Habashi, “Development of a shallow-water icing model in fensap-ice,” Journal of Aircraft 37, 640–646 (2000).
  * Myers (2001) T. G. Myers, “Extension to the messinger model for aircraft icing,” AIAA Journal 39, 211–218 (2001).
  * Wang and Rothmayer (2009) G. Wang and A. Rothmayer, “Thin water films driven by air shear stress through roughness,” Computers & Fluids 38, 235–246 (2009).
  * Liu  _et al._ (2019) T. Liu, K. Qu, J. Cai, and S. Pan, “A three-dimensional aircraft ice accretion model based on the numerical solution of the unsteady stefan problem,” Aerospace Science and Technology 93, 105328 (2019).
  * Knopp, Eisfeld, and Calvo (2009) T. Knopp, B. Eisfeld, and J. B. Calvo, “A new extension for k- turbulence models to account for wall roughness,” International Journal of Heat and Fluid Flow 30, 54–65 (2009).
  * Verdin, Charpin, and Thompson (2009) P. Verdin, J. P. F. Charpin, and C. P. Thompson, “Multistep Results in ICECREMO2,” Journal of Aircraft 46, 1607–1613 (2009).
  * Lavoie (2017) P. Lavoie,  _Modeling of thin water films on swept wings in icing condition_ , Master’s thesis, École Polytechnique de Montréal (2017).
  * Wright and Rutkowski (1999) W. Wright and A. Rutkowski, “Validation results for lewice 2.0,” Tech. Rep. (CR 1999-208690, NASA, 1999).
  * Papadakis  _et al._ (2003) M. Papadakis, H. W. Yeong, M. Vargas, and M. Potapczuk, “Aerodynamic performance of a swept wing with ice accretions,” in  _41st Aerospace Sciences Meeting and Exhibit_ (American Institute of Aeronautics and Astronautics, 2003).
  * Messinger (1953) B. L. Messinger, “Equilibrium temperature of an unheated icing surface as a function of air speed,” Journal of the aeronautical sciences 20, 29–42 (1953).
  * Craster and Matar (2009) R. V. Craster and O. K. Matar, “Dynamics and stability of thin liquid films,” Reviews of Modern Physics 81, 1131–1198 (2009).
  * Hagemeier, Hartmann, and Thévenin (2011) T. Hagemeier, M. Hartmann, and D. Thévenin, “Practice of vehicle soiling investigations: A review,” International Journal of Multiphase Flow 37, 860–875 (2011).
  * Nicholson  _et al._ (2019) J. M. P. Nicholson, H. Power, O. Tammisola, S. Hibberd, and E. D. Kay, “Fluid dynamics of the slip boundary condition for isothermal rimming flow with moderate inertial effects,” Physics of Fluids 31, 033602 (2019).
  * Dong  _et al._ (2016) W. Dong, M. Zheng, J. Zhu, and G. Lei, “Calculation and analysis of runback water flow on anti-icing airfoil surface,” Journal of Aircraft 53, 1597–1605 (2016).
  * Lavalle  _et al._ (2015) G. Lavalle, J.-P. Vila, G. Blanchard, C. Laurent, and F. Charru, “A numerical reduced model for thin liquid films sheared by a gas flow,” Journal of Computational Physics 301, 119–140 (2015).