# A moving boundary approach of capturing diffusants penetration into rubber:
FEM approximation and comparison with laboratory measurements

S. Nepal 1,∗, R. Meyer2, N. H. Kröger2,3, T. Aiki4, A. Muntean1, Y.
Wondmagegne1, U. Giese2
1 Department of Mathematics and Computer Science, Karlstad University, Sweden
2 Deutsches Institut für Kautschuktechnologie e.V. (DIK), Hanover, Germany
3 material prediction GmbH, Wardenburg, Germany
4 Department of Mathematics, Japan Women’s University, Tokyo, Japan
*<EMAIL_ADDRESS>
###### Abstract

We propose a moving-boundary scenario to model the penetration of diffusants
into dense and foamed rubbers. The presented modelling approach recovers
experimental findings related to the diffusion of cyclohexane and the
resulting swelling in a piece of material made of ethylene propylene diene
monomer rubber (EPDM). The main challenge is to find out relatively simple
model components which can mimic the mechanical behavior of the rubber. Such
special structure is identified here so that the computed penetration depths
of the diffusant concentration are within the range of experimental
measurements. We investigate two cases: a dense rubber and a rubber foam, both
made of the same matrix material.

After a brief discussion of scaling arguments, we present a finite element
approximation of the moving boundary problem. To overcome numerical
difficulties due to the a priori unknown motion of the diffusants penetration
front, we transform the governing model equations from the physical domain
with moving unknown boundary to a fixed fictitious domain. We then solve the
transformed equations by the finite element method and explore the robustness
of our approximations with respect to relevant model parameters. Finally, we
discuss numerical estimations of the expected large-time behavior of the
material.

Keywords: Rubber, absorption, diffusion, swelling, moving boundary problem,
finite element method.

##  1 Introduction

Rubber is one of most usable polymers in our daily life, where it is employed
for multiple purposes ranging from tires to materials for daily use, medical
devices, and toys. The source of the natural rubber is the latex of rubber
trees (esp. hevea brasiliensis). Due to its softness, the use of natural
rubber is limited. Besides natural rubber, synthetic rubbers are widely used.
They are usually produced from petroleum byproducts. To develop the wanted
toughness and strength, crosslinking of the rubber molecules, e.g. by
vulcanization, is required. There are several ways to crosslink rubber, but
the most common way is to add sulfur components that help forming bonds
between rubber molecules. This leads to an increased strength as well as to a
lower solubility in solvents; see e.g. [18] and compare also [24] and
references cited therein. On the other hand, due to the solubility of low
molecular components, depending on their polarity, size of the molecules as
well on the flexibility of the polymer chains, the material is to some extent
permeable. This is precisely where our interest lies. Predicting the
penetration of small molecules in polymeric materials has a great importance
in polymer science and engineering especially if one has in mind applications
in areas such as food packaging [20], controlled drug delivery [17], membrane
separation [27], etc. Especially, the a priori knowledge of critical
diffusants and their distribution in rubber matrices is of great interest,
cf.[8, 9, 22, 23]. To determine the approximate position of the penetration
front, many different type of models involving various degrees of complexity
have been proposed earlier; see e.g. [6, 16, 19, 25, 28]. However, a general
consensus is not yet reached in the literature. Consequently, many fundamental
questions are still open particularly what concerns the large-time behavior of
rubber-based materials when they are exposed to environmental conditions. For
the scenario we have in view, this is due to the fact that one does not
understand yet sufficiently well how the macroscopic rheology of the rubber
changes particularly when diffusants-induced swelling takes place.

In this framework, we deal with the macroscopic modeling and numerical
simulation for the penetration of diffusants via absorption and diffusion
through swellable rubbers. We consider two cases: (i) dense rubber and (ii)
foamed rubber; both being made based on the same polymer chemistry.

We are mainly addressing the following questions:

  * (Q1)

How far does the concentration of a diffusant species111We consider a
population of small-size particles able to get absorbed into the rubber’s
skeleton and diffuse from there inside the material. We refer here to this
diffusing population as ”diffusants”. penetrate into a dense rubber?

  * (Q2)

How far are these diffusants able to penetrate a foamed rubber?

Although these questions may appear fairly simple, they are quite hard to
answer in a full generality. We consider in this paper an experimental setup
(see Section 2) that can be traced numerically by means of one dimensional
evolution model with a freely evolving moving boundary - the penetration front
of the diffusants (see Section 3). We address here (Q1) and (Q2) by comparing
the model output with the measured penetration depths and discuss the expected
large-time behavior of both materials when drastic variations of parameters
are applied.

The experimental setup, detailed in Section 2, is sketched in Figure 4. The
model proposed in Section 3 is meant to fit to this precise setup. As a
consequence of the material’s shape, the modeling domain is a 1D line (located
right in the middle of the rubber sheet), which is supposed to swell as time
elapses, due to a continuous penetration of diffusants up to a location at
time $t>0$ denoted by $s(t)$. Predicting the time evolution of this moving
front $s(t)$ is our main target - the evolution takes place along the
longitudinal line indicated in in Figure 4 so that $0<s(0)\leq s(t)\leq\ell$,
where $\ell$ is a characteristic length of the rubber sample. In Section 3, we
describe the governing equations of our model. They are taken over from [13],
where the model was designed to suit our purpose. We start off with performing
the non-dimensionalization of the model equations to identify the typical
sizes of the characteristic times of diffusion, absorption and swelling, and
consequently also the typical size of the mass transfer Biot number. As next
step, for the ease of handling the numerical simulations and parameter
identification procedure, we transform the non-dimensional model equations
from the a priori unknown domain confined by the target diffusants’s
penetration depth $s(t)$ into a fixed domain. In Section 4, we introduce a
finite element approximation of our model, while in Section 5, we present and
discuss the simulation results. We calibrate our model based on the
experimental data corresponding to the dense rubber case. Afterwards in
Section 6, we use the same model and capture the measured diffusion fronts in
the foamed rubber by changing either the kinetic parameters (entering the
speed law for the diffusion front) or the effective diffusivity), or a
suitable combination of these parameters. We close the paper with a series of
discussions and hints to further work, see Section 7 and Section 8. Of
particular interest for us is the study of the expected large-time behavior of
such materials when they are exposed to less controlled environmental
conditions. This work can be seen as a preliminary step in this direction.

##  2 Experimental investigations

The modelling idea is based on the following experimental results related to
ethylene propylene diene monomer rubber (EPDM) being swollen in cyclohexane.
Especially, the difference in the diffusion respectively the swelling
behaviour between full matrix rubber and its foam is of interest.

Material preparation: The material’s ingredients, see Table 1, were admixed at
a temperature of 50∘C and rotor speed of 30 rotation per minute. Each mixture
was vulcanized as plates for 24 minutes at 160∘C in a press. The full matrix
material was produced in 1 mm plates, whereas plates of the foam were produced
releasing the pressure of the vulcanisation press resulting in thicker plates.
Subsequently, the foam plates were cut in 1 mm slices. The geometries used
were 1 mm thick plates with 17 mm width and 40 mm height.

Ingredients in phr | EPDM | ZnO | Stearic Acid | TAIC (70%) | DCP (40%) | TMQ | Backing powder
---|---|---|---|---|---|---|---
Full rubber matrix | 100 | 2.5 | 2.5 | 0.6 | 8 | 0 | 0
Foam | 100 | 2.5 | 2.5 | 0.6 | 8 | 1.5 | 4.7
Table 1: Recipe of the investigated EPDM mixture (phr $\widehat{=}$ parts per
hundred rubber related to mass parts).

Experimental setup and results: A solution of cyclohexane and extracted black
tea as marker is filled in a beaker. In order to compensate the evaporation of
the solvent a dropping funnel is added, and adjusted accordingly to ensure a
constant volume of solution inside the beaker, see Figure 1(left). The needed
dropping rate is determined beforehand in dependency of the evaporation’s
rate. The test specimens are placed in the solution submerged by an initial
length of 20 mm, see Figure 1 (middle). To characterize the diffusion and
swelling behaviour the diffusion front, the swelling length (submerged
specimen) and the submerged swollen area are evaluated, see Figure 1(right),
Figure 2, Figure 3, and Table 2. For evaluation, the specimens are taken out
of the solution, measured and placed back afterwards.

The position of the diffusion front outside the solution saturates for the
rubber material before ten minutes and for the rubber foam before 30 minutes.
The materials swell into the solution as well as it can be seen in Figure 3
which is indicated by the increasing submerged specimens lengths. In case of
the rubber foam its saturation is reached in a similar time frame as for the
diffusion front, whereas the full rubber only slowly converges. Due to the
porosity of the foam and therefore lower (mechanical) resistance the effects
are pronounced. Especially, the diffusion into the upper part of the specimen
is amplified. Similar but slower effects are observed by evaluating the
swollen, submerged surface areas, see Figure 3.

Figure 1: Experimental setup and evaluation. | Rubber | Foamed rubber
---|---|---
Time | Front | Length | Area | Front | Length | Area
[min] | [mm] | [mm] | [mm2] | [mm] | [mm] | [mm2]
0 | 0 | 20 | 340 | 0 | 20 | 340
3.5 | 1 | 21 | 420 | 5 | 25 | 550
10 | 2 | 21 | 420 | 6 | 28 | 700
30 | 2 | 22 | 462 | 7 | 31 | 837
150 | 2 | 24 | 576 | 7 | 31 | 961
300 | 2 | 25 | 625 | 7 | 31 | 961
Table 2: Determined experimental results related to diffusion front, submerged
length and submerged area of the swollen specimens, cf. Figure 1.
$0$$50$$100$$150$$200$$250$$300$$0$$2$$4$$6$$8$Time in minDiffusion front in
mmRubber (front)Foamed rubber (front)$8$$16$$24$$32$Submerged length in
mmRubber (length)Foamed rubber (length) Figure 2: Position of the diffusion
front and length of the submerged rubber, respectively foam, as a function of
time. $0$$50$$100$$150$$200$$250$$300$$0$$200$$400$$600$$800$$1{,}000$Time in
minSwollen, submerged surface area in mm2RubberFoamed
rubber$0$$200$$400$$600$$800$$1{,}000$ Figure 3: Swollen, submerged surface
area of the rubber, respectively foam, as a function of time.

##  3 Model equations

At time $t=0$, we consider a piece of a rubber denoted by $\Omega$ of vertical
length $\ell>0$ placed in contact with a diffusant reservoir. The initial
configuration of the rubber is divided by the interface $x=s(0)$ into two
parts. The lower part $\Omega_{1}(0)=\Omega\cap\\{x<s(0)\\}$ is in contact
with a solvent rich in diffusant concentration, while the upper part
$\Omega_{2}(0)=\Omega\cap\\{x>s(0)\\}$ is diffusant-free. Borrowing
terminology from phase-change problems [4], we start off from the assumption
that the region $\Omega_{1}(0)$, where the diffusant has already penetrated,
is separated by a sharp interface from the diffusant-free region
$\Omega_{2}(0)$. To keep things simple, we look at a symmetric piece of
material so that we can consider the concentration of diffusants as a function
of the one-dimensional space position $x$ and of the time variable
$t\in(0,T)$. Here $T$ is the observation time. As time elapses, the diffusants
penetrate through the small pores of the rubber driven by absorption and
diffusion processes. As a result, the initial position of the interface
$x=s(0)$ changes continuously in time to new locations inwards the material.
We denote the position of the moving diffusant penetration front by $x=s(t)$.
It is worth observing that its location at time $t$ is a priori unknown and
must be resolved together with determining the diffusant concentration inside
the material. The two regions of interest are the diffusants-penetrated part
and the diffusants-free part, viz.

$\Omega_{1}(t):=\Omega\cap\\{x<s(t):\ t\in(0,T)\\}\;\;\text{(with
diffusants)}$

and

$\Omega_{2}(t):=\Omega\cap\\{x>s(t):\ t\in(0,T)\\}\;\;\text{(without
diffusants)},$

respectively. We illustrate in Figure 4 a sketch of the experimental setup we
have in mind.

$\ell$$x=s(0)$$x=s(T)$Rubber (diffusant free)Modeling domain Figure 4: Sketch
of the experimental setup. A rubber material is put in contact with a
reservoir of diffusants.

As $\Omega_{2}(t)$ is diffusants-free zone for each time $t$, the actual
problem is to find the diffusant concentration profile inside $\Omega_{1}(t)$
and the location of the free boundary $s(t)$. This setting is called one-phase
moving boundary problem; see e.g. [10] concerning modeling with moving
interfaces. In this work, the modeling domain is a 1D slab as shown in Figure
5, which is the longitudinal line where $0<s(0)\leq s(t)\leq\ell$.
For a fixed observation time $T\in(0,\infty)$, the interval $[0,T]$ is the
time span of the process we are considering. Let $x\in[0,s(t)]$ and
$t\in[0,T]$ denote the space and respectively time variable, and let $m(t,x)$
be the concentration of diffusant placed in position $x$ at time $t$. The
function $m(t,x)$ acts in the region $Q_{s}(T)$ defined by

$Q_{s}(T):=\\{(t,x)|t\in(0,T)\;\text{and}\;x\in(0,s(t))\\}.$

The evolution of the diffusant concentration is governed by the diffusion
equation

$\displaystyle\frac{\partial m}{\partial t}-D\frac{\partial^{2}m}{\partial
x^{2}}=0\;\;\;\ \text{in}\;\;\;Q_{s}(T),$ (1)

where $D>0$ is the constant diffusion coefficient. This choice indicates that
we assume that the capillary suction does not contribute to the transport of
the diffusants in the rubber matrix; see [14] for a case where the imbibition
of a liquid into a solid matrix is driven by a suitable combination of viscous
drag and capillary pressure.

$0$$s(0)$$s(T)$$\ell$Zone inside rubber with penetration of diffusants at
$t=0$Diffusant-free rubber at $t=0$ Figure 5: Modeling domain.

What concerns boundary conditions, we prescribe at the left boundary $x=0$ an
inflow. This is sometimes referred to as Newton or Robin-type boundary
condition. This is the place where absorption balances the incoming diffusion
flux via

$-D\frac{\partial m}{\partial
x}(t,0)=\beta(b(t)-Hm(t,0))\;\;\;\text{for}\;\;t\in(0,T).$ (2)

In (2), $H>0$ is Henry’s constant, while $\beta>0$ is a given absorption (or
mass transfer) rate. $\beta$ equals zero corresponds to the case where there
is no inflow of diffusants into rubber from the left boundary. Here we
introduce the threshold $b$ to incorporate in the model an approximate value
of the diffusant concentration that is initially present in contact with the
rubber at the bottom surface as indicated in Figure 4. The scenario in which
diffusants tend to diffuse into rubber corresponds to the case where the value
of $b$ is greater than $Hm(0,0)$.

At the right boundary, i.e. at $x=s(t)$, it holds:

$-D\frac{\partial m}{\partial
x}(t,s(t))=s^{\prime}(t)m(t,s(t))\;\;\;\text{for}\;\;t\in(0,T).$ (3)

The boundary condition (3) describes the mass conservation of diffusant
concentration at the moving boundary. It indicates that the diffusion
mechanism is responsible for pushing the interface.
To have a complete model, we need to prescribe one additional condition on the
right boundary as $s(t)$ is a priori unknown. To this end, we propose an
ordinary differential equation that describes explicitly the speed of the
moving boundary by

$s^{\prime}(t)=a_{0}(m(t,s(t))-\sigma(s(t))\;\;\;\;\text{for
}\;\;\;t\in(0,T),$ (4)

where $a_{0}>0$ is constant. The function $\sigma(s(t))$ incorporates an
eventual contribution of the mechanics of the material. In particular (4)
points out that the mechanical behaviour (here it is about the swelling of the
rubber) also contributes to the motion of the moving penetration front. It is
worth mentioning that when the term $\sigma(s(t))$ vanishes, the model
equations form the so-called one-phase Stefan problem with kinetic condition.
This approach was originally used to model the motion of the ice-water
interfaces (compare e.g. [4]) as well as the motion of sharp carbonation
reaction interfaces in concrete (see [2, 3]).
The initial diffusant concentration and the initial position of the moving
front are specified via

$\displaystyle m(0,x)=m_{0}(x)\;\;\;\text{for}\;\;\;x\in[0,s(0)],$ (5)
$\displaystyle s(0)=s_{0}>0\;\text{with}\;\;0<s_{0}<s(t)<L.$ (6)

The model described in this Section was originally proposed in [12] and
analyzed mathematically in [13]. Now is the moment to verify its range of
validity against the measurements reported in Section 2.

###  3.1 Non-dimensionlization

To quantify the relative sizes of the characteristic time scales for
absorption, diffusion and swelling, it is convenient to formulate the model
equations in dimensionless form. To do so, we introduce the dimensionless
objects:

$\displaystyle z=\frac{x}{x_{ref}},\hskip
19.91692pt\tau=\frac{t}{x_{ref}^{2}/D},\hskip
19.91692ptu=\frac{m}{m_{ref}},\hskip 19.91692pth=\frac{s}{x_{ref}},$ (7)

where $x_{ref}$ is a characteristic length scale, while $m_{ref}$ is a
reference value for diffusant concentration. The dimensionless concentration
$u(\tau,z)$ acts in the region $Q_{h}(T^{*})$ defined by

$Q_{h}(T^{*}):=\left\\{(\tau,z)|\tau\in(0,T^{*})\;\text{and}\;z\in(0,h(\tau))\;\text{with}\;T^{*}=TD/x_{ref}^{2}\right\\}.$

The diffusion equation, now in dimensionless form, reads:

$\frac{\partial u}{\partial\tau}-\frac{\partial^{2}u}{\partial z^{2}}=0\;\;\;\
\text{in}\;\;\;Q_{h}(T^{*}).$ (8)

The boundary conditions (2) and (3) become

$\displaystyle-\frac{\partial u}{\partial
z}(\tau,0)=\text{Bi}\left(\frac{b(\tau)}{m_{ref}}-Hu(\tau,0)\right)\;\;\;\text{for}\;\;\tau\in(0,T^{*}),$
(9) $\displaystyle-\frac{\partial u}{\partial
z}\left(\tau,h(\tau)\right)=h^{\prime}(\tau)u\left(\tau,h(\tau)\right)\;\;\;\text{for}\;\;\tau\in(0,T^{*}),$
(10)

where Bi denotes the standard mass transfer Biot number. In our context, this
is defined as

$\text{Bi}:=\beta\frac{x_{ref}}{D}.$ (11)

Equation (4) takes the form

$\displaystyle
h^{\prime}(\tau)=A_{0}\left(u\left(\tau,h(\tau)\right)-\frac{\sigma(h(\tau))}{m_{ref}}\right)\;\;\;\text{for}\;\;\tau\in(0,T^{*}),$
(12)

where

$A_{0}:=\frac{x_{ref}}{D}m_{ref}a_{0}.$ (13)

It is worth mentioning that $A_{0}$ is a sort of Thiele modulus (also called
$2$nd Damköhler number). Both dimensionless numbers Bi and $A_{0}$ are
relevant for the discussion of our results as we will see in Section 5. The
initial conditions become:

$\displaystyle
u(0,z)=\frac{m_{0}(z)}{m_{ref}}\;\;\;\text{for}\;\;\;z\in[0,h(0)],$ (14)
$\displaystyle h(0)=\frac{s_{0}}{x_{ref}}.$ (15)

Note also the following replacements

$\displaystyle b(t_{ref}\tau)\rightarrow b(\tau),$
$\displaystyle\sigma(x_{ref}h(t_{ref}\tau))\rightarrow\sigma(h(\tau)),$
$\displaystyle m_{0}(x_{ref}z)\rightarrow m_{0}(z).$

To fix ideas, we set as reference values

$\displaystyle x_{ref}:=s(T),\;\;\;\;\;m_{ref}:=m_{0},$ (16)

where $m_{0}$ is the initial amount of the solvent, while the value $s(T)$ is
the length of penetration front when the solvent reached the stationary state.

###  3.2 Fixed-domain transformation

To map the system of equations (8)-(15) into the cylindrical domain
$Q(T^{*}):=\\{(\tau,y)|\tau\in(0,T^{*})\;\text{and}\;y\in(0,1)\\}$, we use the
transformation

$y=\frac{z}{h(\tau)},\;\;\;\tau\in(0,T^{*}).$ (17)

This transformation, originally due to Landau cf. [15], maps the mathematical
model with moving boundary conditions into a mathematical model with fixed
boundary conditions. By using the chain rule of differentiation and the
transformation (17), equation (8) becomes

$\displaystyle\frac{\partial{u}}{\partial{\tau}}-y\frac{h^{\prime}(\tau)}{h(\tau)}\frac{\partial{u}}{\partial{y}}-\frac{1}{(h(\tau))^{2}}\frac{\partial^{2}{u}}{\partial{y^{2}}}=0\;\;\;\
\text{in}\;\;\;Q(T^{*}),\\\ $ (18)

where we have used the prime (′) notation to denote the derivative with
respect to time variable $\tau$. The boundary conditions given in (9) and (10)
take now the form

$\displaystyle-\frac{1}{h(\tau)}\frac{\partial u}{\partial
y}(\tau,0)=\text{Bi}\left(\frac{b(\tau)}{m_{0}}-Hu(\tau,0)\right)\;\;\;\text{for}\;\;\tau\in(0,T^{*}),$
(19) $\displaystyle-\frac{1}{h(\tau)}\frac{\partial u}{\partial
y}(\tau,1)=h^{\prime}(\tau)u(\tau,1)\;\;\;\text{for}\;\;\tau\in(0,T^{*}),$
(20)

where $0<y<1.$ The ordinary differential equation corresponding to the growth
rate of the moving interface given in (12) reads

$h^{\prime}(\tau)=A_{0}\left(u(\tau,1)-\frac{\sigma(h(\tau))}{m_{0}}\right)\;\;\;\text{for}\;\;\tau\in(0,T^{*}).$
(21)

The initial conditions (14) and (15) transform into

$\displaystyle u(0,y)=1\;\;\;\text{for}\;\;y\in[0,1],$ (22) $\displaystyle
h(0)=\frac{s_{0}}{x_{ref}}.$ (23)

##  4 Numerical method

To handle the discretization of $\eqref{a32}$ and (21) we use the method of
lines (see, for instance, [26]). First, the equations are discretized in space
by means of the finite element method. The resulting time-dependent ordinary
differential equations (capturing information at the mesh points used for the
finite element method) are solved by the solver odeint in Python, see [11] for
details on the solver. In this Section, we present the details of our
approximation scheme.
The weak formulation of the model is obtained by multiplying (18) with a test
function $\varphi\in H^{1}(0,1)$ and by integrating the result over the domain
$(0,1)$,

$\displaystyle\int_{0}^{1}\frac{\partial u}{\partial\tau}\varphi
dy-\int_{0}^{1}y\frac{h^{\prime}(\tau)}{h(\tau)}\frac{\partial u}{\partial
y}\varphi
dy-\int_{0}^{1}\frac{1}{(h(\tau))^{2}}\frac{\partial^{2}{u}}{\partial{y^{2}}}\varphi
dy=0.$ (24)

Using integration by parts in (24) and the boundary conditions (19) and (20),
we get the following weak formulation: Find the couple $(u,h)$ satisfying the
equations

$\displaystyle\displaystyle\int_{0}^{1}\frac{\partial u}{\partial\tau}\varphi
dy-\int_{0}^{1}y\frac{h^{\prime}(\tau)}{h(\tau)}\frac{\partial u}{\partial
y}\varphi
dy+\int_{0}^{1}\frac{1}{(h(\tau))^{2}}\frac{\partial{u}}{\partial{y}}\frac{\partial\varphi}{\partial
y}dy$ $\displaystyle\hskip
113.81102pt-\frac{1}{h(\tau)}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-Hu(\tau,0)\right)\varphi(0)+\frac{h^{\prime}(\tau)}{h(\tau)}u(\tau,1)\varphi(1)=0\;\;\text{for
all}\;\varphi\in H^{1}(0,1),$ (25) $\displaystyle
h^{\prime}(\tau)=A_{0}\left(u(\tau,1)-\frac{\sigma(h(\tau))}{m_{0}}\right),\;\;\tau\in(0,T^{*}),$
(26)

subject to the initial conditions (22) and (23).

Existence and uniqueness of weak solution to problem (25)-(26) are discussed
in [13]. We refer the reader to Theorem 3.1 and Theorem 3.2 in [13] for the
details on these matters. Furthermore, it turns out that the weak solution
depends in a continuous way with respect to the model parameters. Now, we
present approximation of this weak solution numerically.

We use here a Galerkin scheme based on piecewise linear polynomials to
discretize the model equations (25) and (26). For $N\in\\{2,3,\cdots\\}$, let
$\\{\phi_{j}\\}_{j=0}^{N-1}$ denote the standard piecewise linear continuous
basis function (i.e., the ”hats basis”) for the finite dimensional subspace
$V_{k}$ of $H^{1}(0,1)$ defined on the uniform mesh
$[0,1/(N-1),2/(N-1),\cdots,1]$. The mesh size is $k:=1/(N-1)$. For
$j\in\\{0,1,\cdots,N-1\\}$, we have

$\phi_{j}(y)=\begin{cases}\displaystyle\;\;\frac{y-y_{j-1}}{k},&\text{if\;
$y_{j-1}\leq y\leq y_{j}$},\\\
\displaystyle\;\;\frac{y_{j+1}-y}{k},&\text{if\; $y_{j}\leq y\leq
y_{j+1}$},\\\ \;\;\;0,&\text{else}.\end{cases}$

$\phi_{0}$$\phi_{1}$$\phi_{2}$$\phi_{3}$$\phi_{4}$$0=y_{0}$$y_{1}$$y_{2}$$y_{3}$$y_{4}=1$$1$
Figure 6: Basis functions for $V_{k}$

We formulate an approximate problem as follows: Find the pair $(u_{k},h)\in
C(0,T^{*};V_{k})\times C^{1}(0,T^{*})$ such that

$\displaystyle\displaystyle\int_{0}^{1}\frac{\partial u_{k}}{\partial\tau}\chi
dy-\int_{0}^{1}y\frac{h^{\prime}(\tau)}{h(\tau)}\frac{\partial u_{k}}{\partial
y}\chi
dy+\int_{0}^{1}\frac{1}{(h(\tau))^{2}}\frac{\partial{u_{k}}}{\partial{y}}\frac{\partial\chi}{\partial
y}dy$ $\displaystyle\hskip
56.9055pt-\frac{1}{h(\tau)}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-Hu_{k}(\tau,0)\right)\chi(0)+\frac{h^{\prime}(\tau)}{h(\tau)}u_{k}(\tau,1)\chi(1)=0,\;\;\text{for
all}\;\;\chi\in V_{k},\;\;\tau\in(0,T^{*})$ (27) $\displaystyle
h^{\prime}(\tau)=A_{0}\left(u_{k}(\tau,1)-\frac{\sigma(h(\tau))}{m_{0}}\right),\tau\in(0,T^{*}),$
(28) $\displaystyle u_{k}(0)=1,$ (29) $\displaystyle
h(0)=\frac{s_{0}}{x_{ref}}.$ (30)

We let

$u_{k}(\tau,y)=\sum_{j=0}^{N-1}\alpha_{j}(\tau)\phi_{j}(y),\;\;\tau\in[0,T^{*}],\;\;y\in[0,1],$
(31)

where

$\alpha(\tau)=(\alpha_{0}(\tau),\alpha_{1}(\tau),\cdots,\alpha_{N-1}(\tau))^{T}\in\mathbb{R}^{N}.$

For $\tau=0$, $\alpha_{j}(0),\;\;j\in\\{0,1,\cdots,N-1\\}$ are the nodal
values given by the initial approximation, i.e. for
$j\in\\{0,1,\cdots,N-1\\}$, it holds:

$\alpha_{j}(0)=1.$ (32)

Substituting (31) in (27) yields

$\displaystyle\sum_{j=0}^{N-1}\alpha_{j}^{\prime}(\tau)(\phi_{j}(y),\varphi)$
$\displaystyle-\frac{h^{\prime}(\tau)}{h(\tau)}\sum_{j=0}^{N-1}\alpha_{j}(\tau)b(\phi_{j}(y),\varphi)+\frac{1}{(h(\tau))^{2}}\sum_{j=0}^{N-1}\alpha_{j}(\tau)a(\phi_{j}(y),\varphi)$
$\displaystyle-\frac{1}{h(\tau)}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-H\sum_{j=0}^{N-1}\alpha_{j}(\tau)\phi_{j}(0)\right)\varphi(0)$
(33)
$\displaystyle+\frac{h^{\prime}(\tau)}{h(\tau)}\sum_{j=0}^{N-1}\alpha_{j}(\tau)\phi_{j}(1)\varphi(1)=0,\;\;\tau\in(0,T^{*}).$
(34)

Taking, in particular, $\chi=\phi_{i}(y)\in V_{k}$, for
$i\in\\{0,1,\cdots,N-1\\}$ in (33), we obtain

$\displaystyle\sum_{j=0}^{N-1}\alpha_{j}^{\prime}(\tau)(\phi_{j}(y),\phi_{i}(y))$
$\displaystyle-\frac{h^{\prime}(\tau)}{h(\tau)}\sum_{j=0}^{N-1}\alpha_{j}(\tau)b(\phi_{j}(y),\phi_{i}(y))$
$\displaystyle+\frac{1}{(h(\tau))^{2}}\sum_{j=0}^{N-1}\alpha_{j}(\tau)a(\phi_{j}(y),\phi_{i}(y))-\frac{1}{h(\tau)}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-H\sum_{j=0}^{N-1}\alpha_{j}(\tau)\phi_{j}(0)\right)\phi_{i}(0)$
$\displaystyle+\frac{h^{\prime}(\tau)}{h(\tau)}\sum_{j=0}^{N-1}\alpha_{j}(\tau)\phi_{j}(1)\phi_{i}(1)=0,\;\;\tau\in(0,T^{*})\;\;\;\text{for
all}\;\;\;i\in\\{0,1,\cdots,N-1\\}.$ (35)

We write this expression in matrix form as

$\displaystyle M\alpha^{\prime}(\tau)$
$\displaystyle-\frac{h^{\prime}(\tau)}{h(\tau)}K\alpha(\tau)+\frac{1}{(h(\tau))^{2}}A\alpha(\tau)$
$\displaystyle-\frac{1}{h(\tau)}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-H\alpha_{0}(\tau)\right)\textbf{e}_{0}+\frac{h^{\prime}(\tau)}{h(\tau)}\alpha_{N-1}(\tau)\textbf{e}_{N-1}=0,\;\;\tau\in(0,T^{*}),$
(36)

where
$\alpha(\tau)=(\alpha_{0}(\tau),\alpha_{1}(\tau),\cdots,\alpha_{N-1}(\tau))^{T}\in\mathbb{R}^{N}$
is the unknown vector. The matrices $M,K,A\in\mathbb{R}^{N\times N}$ are
defined by

$\displaystyle
M:=M_{i,j}=\int_{0}^{1}\phi_{i}(y)\phi_{j}(y)dy,\;\;\;i,j\in\\{0,1,2,\dots,N-1\\},$
$\displaystyle M=\frac{k}{6}\begin{bmatrix}2&1&0&0&0&\cdots&0&0&0&0\\\
1&4&1&0&0&\cdots&0&0&0&0\\\ 0&1&4&1&0&\cdots&0&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&0&\cdots&0&1&4&1&0\\\ 0&0&0&0&\cdots&0&0&1&4&1\\\
0&0&0&0&\cdots&0&0&0&1&2\\\ \end{bmatrix},$ $\displaystyle
K:=K_{i,j}=\int_{0}^{1}y\phi_{i}(y)\phi_{j}^{\prime}(y)dy,\;\;\;i,j\in\\{0,1,2,\dots,N-1\\},$
$\displaystyle
K=\begin{bmatrix}\frac{2k-3y_{1}}{6}&-\left(\frac{2k-3y_{1}}{6}\right)&0&0&\cdots&0&0&0\\\
-\left(\frac{2k-3y_{0}}{6}\right)&\frac{4k-3(y_{0}+y_{2})}{6}&-\left(\frac{2k-3y_{2}}{6}\right)&0&\cdots&0&0&0\\\
0&-\left(\frac{2k-3y_{1}}{6}\right)&\frac{4k-3(y_{1}+y_{3})}{6}&-\left(\frac{2k-3y_{3}}{6}\right)&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&0&\cdots&-\left(\frac{2k-3y_{N-3}}{6}\right)&\frac{4k-3(y_{N-3}+y_{N-1})}{6}&-\left(\frac{2k-3y_{N-1}}{6}\right)\\\
0&0&0&0&\cdots&0&-\left(\frac{2k-3y_{N-2}}{6}\right)&\frac{2k-3y_{N-2}}{6}\\\
\end{bmatrix},$ $\displaystyle
A:=A_{i,j}=\int_{0}^{1}\nabla\phi_{i}(y)\nabla\phi_{j}(y)dy,\;\;\;i,j\in\\{0,1,2,\dots,N-1\\},$
$\displaystyle A=\frac{1}{k}\begin{bmatrix}1&-1&0&0&0&\cdots&0&0&0&0\\\
-1&2&-1&0&0&\cdots&0&0&0&0\\\ 0&-1&2&-1&0&\cdots&0&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&0&\cdots&0&-1&2&-1&0\\\ 0&0&0&0&\cdots&0&0&-1&2&-1\\\
0&0&0&0&\cdots&0&0&0&-1&1\\\ \end{bmatrix}.$

The vectors $\textbf{e}_{0},\textbf{e}_{N-1}\in\mathbb{R}^{N}$ are given by

$\displaystyle\textbf{e}_{0}=(1,0,\cdots,0,0)^{T}\;\text{and}\;\;\textbf{e}_{N-1}=(0,0,\cdots,0,1)^{T}.$

Note that the matrix $M$ is symmetric and positive definite as well as
diagonally dominant, and hence, invertible.

$\displaystyle\alpha^{\prime}(\tau)=\frac{h^{\prime}(\tau)}{h(\tau)}M^{-1}K\alpha(\tau)-\frac{1}{(h(\tau))^{2}}M^{-1}A\alpha(\tau)$
$\displaystyle\hskip
85.35826pt+\frac{1}{h(\tau)}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-H\alpha_{0}(\tau)\right)M^{-1}\textbf{e}_{1}-\frac{h^{\prime}(\tau)}{h(\tau)}\alpha_{N-1}(\tau)M^{-1}\textbf{e}_{N-1},\;\;\tau\in(0,T^{*}).$
(37)

Substituting (31) in (28), we get

$\displaystyle
h^{\prime}(\tau)=A_{0}\left(\alpha_{N-1}(\tau)-\frac{\sigma(h(\tau))}{m_{0}}\right),\;\tau\in(0,T^{*}).$
(38)

Equation (4) together with (38) forms a system of $N+1$ ordinary differential
equations in terms of the unknowns
$(\alpha_{0}(\tau),\alpha_{1}(\tau),\cdots,\alpha_{N-1}(\tau),h(\tau))$. As
the system is not stiff, nor it contains algebraic constraints, it can be
solved numerically by standard solvers such as odeint in Python, relying on
the initial condition
$(\alpha_{0}(0),\alpha_{1}(0),\cdots,\alpha_{N-1}(0),h(0))$ given in (32) and
(30).

##  5 Model calibration. Simulation results for the dense rubber

In this section we present our simulation results for the dense rubber case.
Basically, we calibrate here our model to represent the experimental range.
The situation of the rubber foam is the aim of the next section.

The output consists of the concentration profile of the diffusants and of the
position of the moving front. In this framework, we investigate the parameter
space by exploring eventual effects of the choice of parameters on the overall
diffusants penetration process. Essentially, we deal with (4) and (38) relying
on repeated use of the solver odeint in Python. Our implementation requires
three input arguments which are a vector of first order ordinary differential
equations, a vector of initial conditions, and a finite set of time
discretization points. The output is a matrix where each row is the solution
vector of the ordinary differential equations defined at each discretization
time point, producing the initial conditions vector in the first row of the
output matrix. The simulation results are obtained using the set of reference
parameters given in Table 3.

Parameters |  Dimension |  Typical Values
---|---|---
Diffusion constant for concentration in rubber, $D$ |  $L^{2}T^{-1}$ |  $3.66\times 10^{-4}$ (mm2/min), [18]
Absorption rate, $\beta$ |  $LT^{-1}$ |  $0.564$ (mm/min), [21]
Constant appearing in the speed of the moving boundary $a_{0}$ |  $L^{4}T^{-1}M^{-1}$ |  $500$ (mm4/min/gram)
Initial height of diffusants (experiment), $s_{0}$ |  $L$ |  $0.01$ (mm)
$\sigma(s(t))$ |  $ML^{-3}$ |  $\frac{s(t)}{10}$ (gram/mm3)
Initial diffusant concentration (experiment), $m_{0}$ |  $ML^{-3}$ |  $0.1$ (gram/mm3)
Concentration in lower surface of the rubber, $b$ |  $ML^{-3}$ |  $1$ (gram/mm3)
Henry’s constant, $H$ |  – |  2.50 (dimensionless), [5]
Table 3: Name, dimension and typical values for the model parameters.

Initially, the concentration is uniformly distributed within the rubber up to
$0.01$ mm. We take as observation time $T=40$ minutes for the final time with
time step $\Delta t=1/1000$ minute. We choose the number of space
discretization points $N$ to be $100$. We take the value $10$ mm for
characteristic length scale $x_{ref}$. As reference diffusant concentration
$m_{ref}$, we choose $0.1$ gram/mm3. With our choice of parameters, the
dimensionless numbers Bi and $A_{0}$ defined in (11) and (13), are of order of
$O(10^{4})$ and of $O(10^{6})$, respectively.

The crucial component of the model is (4) which describes the speed of the
moving front. To get insight into the role played by the parameters in (4), we
perform simulation for different values of $\sigma(s(t))$ and $a_{0}$.

###  5.1 Approximated swelling - capturing the effect of $\sigma(s(t))$

We present numerical results for different levels of linearity of
$\sigma(s(t))$. In Figure 7 we show the numerical results for the
concentration profile for different values of $\sigma(s(t))$, while taking
$a_{0}=500$. The first plot in Figure 7 shows a typical large-time behaviour
very similar to the classical Stefan problem (ice-melting problem) where the
concentration of the diffusant at the moving front position approaches zero.
Comparing now the remaining two plots in Figure 7, we see that in this
scenario the function $\sigma(s(t))$ describes a breaking mechanism preventing
diffusants to move further inside the material. As a consequence of this fact,
we notice an accumulation of diffusants at the position of the moving front.

Figure 7: Concentration vs. space with
$\sigma(s(t))=\frac{s(t)}{20},\;\sigma(s(t))=\frac{s(t)}{10},\;\sigma(s(t))=\frac{s(t)}{5}$
(from left to right).

It is worth noting that such an accumulation effect can be used in principle
for material design purposes, since it gives an indication on how to bring a
certain diffusant concentration in a prescribed region inside the material
within a given time span.

Considering all three linear cases for $\sigma(s(t))$, within a short time of
release of diffusant from its initial position, the diffusant quickly enters
the rubber from the left boundary (as we expect from (3)) and then start off
to diffuse further raising the penetration front. In the same Figure 7, we see
that close to the initial time $t=0$, concentration profiles change convexity.
This situation disappears as time elapses. Such feature seems to be due to our
choice of the parameter regime. In this particular case, a transient layer
occurs in the time span of the process close to $t=0$. During this time span,
there is apparently a rather strong competition between diffusion, absorption,
and swelling. The experimental data reported in Table 2 refers to precisely
this transient time during which most effects are visible. As expected, our
numerical simulations indicate that the steady state of the process is not
reached during the transient layer.

In Figure 8 we compare the numerical results and experimental data for the
moving front. The first and third plot in Figure 8 show a large deviation
between numerical results and experimental data whereas the second plot shows
a good agreement between numerical results and experimental data.

Figure 8: Comparison of experimental data and numerical diffusion front with
different choices of $\sigma(s(t))$ for $T=40$ minutes.

###  5.2 Effect of the kinetic parameter $a_{0}$

To get some insight in the role played by the parameter $a_{0}$, we perform
simulation runs for different values of $a_{0}$. As illustrated in Figure 9,
the changes in $a_{0}$ lead to changes in the numerical output. More
precisely, larger values of $a_{0}$ lead to a faster penetration front, as
expected in fact from (4).

Figure 9: Comparison of diffusion front for fixed
$\sigma(s(t))=\frac{s(t)}{10}$ with varying $a_{0}=10,\;100,\;1000$.

To estimate the dependence on time of the moving front $s(t)$, we fit our
numerical results to the form

$s(t)=t^{\gamma},$ (39)

where $\gamma>0.$ In Figure 10, we show log-log plots of penetration fronts
for different values of $\sigma(s(t))$ and $a_{0}$ and compare them against
(39). The slope of a log-log plot gives the value for $\gamma$, and a straight
line is an indication that relationship (39) holds true. The approximated
values for $\gamma$ with varying the levels of linearity of $\sigma(s(t))$ and
$a_{0}$ are listed in Table 4 and Table 6, respectively. We expect that the
actual value of $\gamma$ strongly depends on the choice of structure of
$\sigma(s(t))$.

Figure 10: Short-time behaviour of dense rubber: comparison of
$\log(\sqrt{t})$ vs. $\log(t)$ with $\log(s(t))$ vs. $\log(t)$ for different
values of $\sigma(s(t))$ with $T=40$ minutes and
$a_{0}=100,\;a_{0}=500,\;a_{0}=1000$ (from left to right).

Considering now the information exhibited in Figure 10, we see comparisons of
the position of the front’s penetration depth versus time for various
modifications of the kinetic parameter $a_{0}$ as well as of the levels of
linearity in the swelling term $\sigma(t)$. We compare all of them against the
pure diffusive behavior $\sigma(t)\approx\sqrt{t}$. In other words, any
deviation from a $\sqrt{t}$-like behavior can now be attributed to swelling.
Figure 10 pinpoints what happens at the laboratory scale. Choosing the log-log
scale in these plots allow us to identify the best exponents $\gamma$, which
we list in Table 4. At this stage, we can only say that if $\gamma>0.5$ the
motion of the diffusants penetration front appears to be super-diffusive,
while for $\gamma<0.5$ this is sub-diffusive.

$a_{0}$ $\sigma(s(t))$ |  $\frac{s(t)}{20}$ |  $\frac{s(t)}{10}$ |  $\frac{s(t)}{5}$
---|---|---|---
100 |  0.25684094 |  0.18829908 |  0.11140718
500 |  0.3726316 |  0.26559534 |  0.16501296
1000 |  0.40093566 |  0.28627561 |  0.18102778
Table 4: For the dense rubber: approximated value of $\gamma$ entering
$s(t)=t^{\beta}$ for $T=40$ minutes.

##  6 Simulation results for the foam rubber case

In this section, we perform numerical simulations to find out robust choices
of $\sigma(s(t))$ and $a_{0}$ to mimic the experimental data for the foam
rubber. Except for $\sigma(s(t))$ and $a_{0}$, we use all other parameter
values same as for dense rubber case discussed in Section 5. In Figure 11, we
show the concentration profile for different values of $\sigma(s(t))$, while
taking $a_{0}=2000$. With our choice of parameters, the dimensionless numbers
Bi is in the same range as before, but $A_{0}$ is perhaps $1$ order of
magnitude higher. Comparing the plots in Figure 11 with Figure 7, we see that
the concentration profiles for the dense rubber and foam rubber cases have a
similar shape.

Figure 11: For the foam rubber case: concentration vs. space with
$\sigma(s(t))=\frac{s(t)}{100},\;\sigma(s(t))=\frac{s(t)}{50},\;\sigma(s(t))=\frac{s(t)}{25}$
(from left to right).

In Figure 12, we compare our numerical results with the experimental data for
the moving front. One can observe some deviation between the numerical results
and the experimental data in the case of the solid line ( ) and dashdot ( ),
whereas the dashed line ( ) shows a good agreement.

Figure 12: For the foam rubber case: comparison of the experimental diffusion
front with numerical diffusion front for $a_{0}=2000$ and different choices of
$\sigma(s(t))$.

We also investigate the importance of the parameter $a_{0}$ what concerns the
case of the foam rubber. We compare the experimentally measured diffusion
fronts against simulated diffusion fronts obtained when choosing different
values of $a_{0}$ as shown in Figure 13. Comparing the plots corresponding to
$a_{0}=2500$ and $a_{0}=4000$, the effect of $a_{0}$ on the diffusion front is
mainly seen during the first few minutes of the release of the concentration
from its initial position, but not too much effect can be seen when time
elapses beyond a couple of minutes.

Figure 13: For the foam rubber case: comparison of the experimental diffusion
front with numerical diffusion front for $\sigma(s(t))=\frac{s(t)}{50}$ and
different choices of $a_{0}$.

We wish to see how much the asymptotic behavior of the diffusion fronts
deviates from the $\sqrt{t}$-law. To achieve this task, we fit the simulated
diffusion front $s(t)$ to the form (39). We list in Table 5 and Table 7 the
best fit values of $\gamma$ for different values of $\sigma(s(t))$ and
$a_{0}$. Figure 14 and Figure 16 point out log-log plots of penetration fronts
for different values of $\sigma(s(t))$ and $a_{0}$ compared against (39).

Figure 14: Short-time behaviour of foam rubber : comparison of $\log(\sqrt{t})$ vs. $\log(t)$ with $\log(s(t))$ vs. $\log(t)$ for different values of $\sigma(s(t))$ with $T=40$ minutes and $a_{0}=500,\;a_{0}=1000,\;a_{0}=2000$ (from left to right). $a_{0}$ $\sigma(s(t))$ |  $\frac{s(t)}{100}$ |  $\frac{s(t)}{50}$ |  $\frac{s(t)}{25}$
---|---|---|---
500 |  0.63059691 |  0.52306905 |  0.40861496
1000 |  0.70484987 |  0.56983803 |  0.44029209
2000 |  0.74367275 |  0.59585696 |  0.4592986
Table 5: For the foam rubber case: approximated value of $\gamma$ entering
$s(t)=t^{\gamma}$ for $T=40$ minutes.

##  7 Expected large-time behavior of penetration fronts

We are interested in predicting the behavior of penetration depths of the
diffusants beyond the laboratory timescale, hence the large-time asymptotic is
now of interest. To perform such large-time numerical simulations, we rely on
the penetration depths listed in Table 2 to which we add an additional point
recorded after two days. This additional information was obtained in another
large-time measurement performed for the same type of rubber with a similar
solvent. The role of this last measurement point is to ensure that numerical
results stay within the expected physical range as long as possible.

As we have already shown the concentration profile for large time in Figure 7
(for the dense rubber case) and Figure 11 (for the foam rubber case), we now
study the large time behaviour of diffusion fronts. We refer the reader to see
our results for the dense rubber and foam rubber in Figure 15 and Figure 16,
respectively. We list in Table 6 and in Table 7 various best options for the
exponent $\gamma$ matching the two choices of materials.

Figure 15: For the dense rubber case: comparison of $\log(\sqrt{t})$ vs.
$\log(t)$ with $\log(s(t))$ vs. $\log(t)$ for different values of
$\sigma(s(t))$ with $T=5000$ minutes and $a_{0}=100,\;a_{0}=500,\;a_{0}=1000$
(from left to right).

Figure 16:  For the foam rubber case: Comparison of $\log(\sqrt{t})$ vs. $\log(t)$ with $\log(s(t))$ vs. $\log(t)$ for different values of $\sigma(s(t))$ with $T=5000$ minutes and $a_{0}=500,\;a_{0}=1000,\;a_{0}=2000$ (from left to right). $a_{0}$ $\sigma(s(t))$ |  $\frac{s(t)}{20}$ |  $\frac{s(t)}{10}$ |  $\frac{s(t)}{5}$
---|---|---|---
100 |  0.30145999 |  0.2694511 |  0.23967876
500 |  0.30531516 |  0.26961346 |  0.23571825
1000 |  0.30841906 |  0.27189137 |  0.23634618
Table 6: For the dense rubber case: approximated value of $\gamma$ entering $s(t)=t^{\gamma}$ for $T=5000$ minutes. $a_{0}$ $\sigma(s(t))$ |  $\frac{s(t)}{100}$ |  $\frac{s(t)}{50}$ |  $\frac{s(t)}{25}$
---|---|---|---
500 |  0.39449518 |  0.35425142 |  0.31698023
1000 |  0.39890782 |  0.35803523 |  0.32025865
2000 |  0.40184547 |  0.3605881 |  0.32252626
Table 7: For the foam rubber case: approximated value of $\gamma$ entering
$s(t)=t^{\gamma}$ for $T=5000$ minutes.

Looking for instance at Figure 9 (for the dense rubber case) or at Figure 13
and Figure 14 (for the foam rubber case), we notice that the penetration front
$s(t)$ really behaves like $t^{\gamma}$ as $t$ is sufficiently large with
$\gamma$ some positive number. We also notice that the dimensionless numbers
Bi and $A_{0}$ are invariant with respect to time. This makes us conjecture
that, under suitable conditions, there exist two positive constants $c_{1}$
and $c_{2}$ (independent of the choice of $t$) such that the inequality

$c_{1}t^{\gamma}\leq s(t)\leq c_{2}t^{\gamma}$ (40)

holds true for sufficiently large $t$. In Refs. [2, 3], the authors were able
to prove such estimate on the large time behavior of the moving boundary for a
different set of moving boundary conditions combined with a Dirichlet boundary
condition at $x=0$. We expect that some of the techniques used in [2] are in
principle applicable in our context as well. Hence, we hope to be able to
prove rigorously the conjecture (40) in a follow-up paper. It is worth noting
that if the conjecture were true, then new opportunities would appear
particularly regarding the reduction of complex modeling approaches to the use
of simple (material-dependent) educated guesses for fitted power laws to
capture asymptotically the penetration of the diffusants.

##  8 Conclusion and outlook

Within this framework, we have proposed a one dimensional moving boundary
model capable to capture the motion of sharp diffusion fronts in both dense
and foam rubbers. Including in our modeling work information about the
measured swelling of the area and about the measured elongated height of the
sample can be used further to devise a balance law for the porosity change due
to the penetration of the diffusant. This would implicitly include nonlinear
diffusion effects. Such features need to be accounted for especially if one
wants to capture swelling effects in more than one space dimension. Moreover,
it is of both engineering and mathematical importance to build simple
asymptotic laws that can be verified experimentally in laboratory conditions
and have a mathematically-controlled range of application outside the
laboratory. From this point of view, it is important to understand the nature
of the deviations that usually occur in the position of the penetration front
$s(t)$ from a $\sqrt{t}$-like behavior. As further possible outlet, we plan to
extend the current model to a multiscale framework, possibly based on ideas
from [1]. This approach would allow us to include explicitly in our numerical
investigations various microstructure models for the rubber foam.

## Notation

We list here the notation used in the text.
$A=N\times N$ tridiagonal matrix,
$A_{0}$ = Thiele modulus (dimensionless),
Bi = Biot number (dimensionless),
$C(0,T^{*};V_{k})=$ the space of continuous function from $(0,T^{*})$ to
$V_{k}$ [26],
$C^{1}(0,T^{*})=$ the space of continuously differentiable function defined on
$(0,T^{*})$,
$D$ = diffusion constant for concentration in rubber (mm2/min),
$H$ = Henry’s constant (dimensionless),
$H^{1}(0,1)$ = Sobolev space [26],
$K=N\times N$ tridiagonal matrix,
$M=N\times N$ positive definite matrix,
$V_{k}=$ finite dimensional subspace of $H^{1}(0,1)$ [26],
$a_{0}$ = constant appearing in the speed of the moving front (mm4/min/gram),
$b$ = concentration in the lower surface of the rubber (gram/mm3),
$\textbf{e}_{0},\;\textbf{e}_{N-1}$ = unit vectors,
$h(\tau)$ = position at time $\tau$ of the front separating a diffusant-free
zone from a diffusant-penetrated zone within the rubber between concentration
and rubber, dimensionless form,
$m$= diffusant concentration (gram/mm3),
$m_{0}$ = initial value of diffusant concentration $m$ (gram/mm3),
$s(t)$ = position of the front between concentration and rubber at time $t$
(mm),
$s_{0}=s(0)$ (mm),
$s^{\prime}(t)$ = growth rate of moving boundary (mm/min),
$u$ = dimensionless form of diffusant concentration $m$,
$\beta$ = mass transfer constant from concentration to rubber at the boundary
(mm/min),
$\phi_{j}$ = linear piece-wise continuous basis function.

## Acknowledgements

A.M. thanks K. Kumazaki (Nagasaki, Japan) for very fruitful discussions on
this topic. The work of S.N. and A.M. is financed partly by the VR grant
2018-03648. T.A. and A.M. thank the Knowledge Foundation (project nr. KK
2019-0213) for supporting financially their research. The work of T.A. is
partially supported also by JSPS KAKENHI project nr. JP19K03572.

## References

  * [1] T. Aiki, N. H. Kröger, and A. Muntean.  A macro-micro elasticity-diffusion system modeling absorption-induced swelling in rubber foams – proof of the strong solvability.  arXiv, 2010.03153, 2020.
  * [2] T. Aiki and A. Muntean.  A free-boundary problem for concrete carbonation: Front nucleation and rigorous justification of the $\sqrt{t}$-law of propagation.  Interfaces and Free Boundaries, 15:167–180, 2013.
  * [3] T. Aiki and A. Muntean.  Large-time asymptotics of moving-reaction interfaces involving nonlinear henry’s law and time-dependent dirichlet data.  Nonlinear Analysis, Theory, Methods & Applications, 93:3–14, 2013\.
  * [4] V. Alexiades.  Mathematical Modeling of Melting and Freezing Processes.  CRC Press, 1992.
  * [5] M. Böhm, J. Devinny, F. Jahani, and G. Rosen.  On a moving-boundary system modeling corrosion in sewer pipes.  Applied Mathematics and Computation, 92(2-3):247–269, 1998.
  * [6] A. Fasano, G. Meyer, and M. Primicerio.  On a problem in the polymer industry: theoretical and numerical investigation of swelling.  SIAM J. Appl. Math., 17:945–960, 1986.
  * [7] I. Gawel, R. Stepkowski, and F. Czechowski.  Molecular interactions between rubber and asphalt.  Industrial & Engineering Chemistry Research, 45(9):3044–3049, 2006\.
  * [8] U. Giese, R. Hebisch, K.-H. Pannwitz, and M. Tschickard.  Theoretische Grundlagen der Passivprobenahme.  In A. Hartwig, editor, Luftanalysen, Analytische Methoden zur Prüfung gesundheitsschädlicher Arbeitsstoffe, pages 81–103. 15. Lieferung (DFG-Publikationen), Band 1, 2000.
  * [9] U. Giese, V. Peterseim, and R. H. Schuster.  Charakterisierung der Verteilung von Mineralölweichmacheranteilen in Kautschukverschnitten.  GAK: Gummi, Fasern, Kunststoffe; Fachmagazin für die Polymerindustrie, 51:729–733, 1998.
  * [10] S. C. Gupta.  The Classical Stefan Problem: Basic Concepts, Modelling and Analysis with Quasi-analytical Solutions and Methods, volume 45.  Elsevier, 2017.
  * [11] R. Johansson.  Numerical Python: Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib.  Apress, 2018.
  * [12] K. Kumazaki and A. Muntean.  Local weak solvability of a moving boundary problem describing swelling along a halfline.  Networks & Heterogeneous Media, 14:445–469, 2019.
  * [13] K. Kumazaki and A. Muntean.  Global weak solvability, continuous dependence on data, and large time growth of swelling moving interfaces.  Interfaces and Free Boundaries, 22(1):27–50, 2020.
  * [14] M. Kvick, Martinez M.D., D.R. Hewitt, and N.J. Balmforth.  Imbibition with swelling: Capillary rise in thin deformable porous media.  Physical Review Fluids, 2:074001, 2017.
  * [15] H. G. Landau.  Heat conduction in a melting solid.  Quarterly of Applied Mathematics, 8(1):81–94, 1950.
  * [16] R. Masoodi, H. Tan, and K. M. Pillai.  Numerical simulation of liquid absorption in paper-like swelling porous media.  AIChE Journal, 58(8):2536–2544, 2012.
  * [17] C. Mircioiu, V. Voicu, V. Anuta, A. Tudose, C. Celia, D. Paolino, M. Fresta, R. Sandulovici, and I. Mircioiu.  Mathematical modeling of release kinetics from supramolecular drug delivery systems.  Pharmaceutics, 11(3):140, 2019.
  * [18] M. Morton.  Rubber Technology.  Springer Science & Business Media, 2013.
  * [19] F. Neff, A. Lion, and M. Johlitz.  Modelling diffusion induced swelling behaviour of natural rubber in an organic liquid.  Journal of Applied Mathematics and Mechanics, 99(3):e201700280, 2019\.
  * [20] A. Reynier, P. Dole, and A. Feigenbaum.  Migration of additives from polymers into food simulants: numerical solution of a mathematical model taking into account food and polymer interactions.  Food Additives & Contaminants, 19(1):89–102, 2002.
  * [21] M. G. Rezk and J. Foroozesh.  Determination of mass transfer parameters and swelling factor of co2-oil systems at high pressures.  International Journal of Heat and Mass Transfer, 126:380–390, 2018\.
  * [22] C. Rosca, U. Giese, R. H. Schuster, and T. Früh.  Diffusion of dithiophosphate accelerators in tire compounds.  Kautschuk Gummi Kunststoffe, 57:593–598, 2004.
  * [23] C. Rosca, U. Giese, R. H. Schuster, and T. Früh.  Investigation of diffusion of phthalates in nitrile rubber by means of FT-IR-Spectroscopy.  Kautschuk Gummi Kunststoffe, 59:86–92, 2006.
  * [24] N. Sombatsompop and P. Lertkamolsin.  Effects of chemical blowing agents on swelling properties of expanded elastomers.  Journal of Elastomers & Plastics, 32(4):311–328, 2000.
  * [25] T. Sweijen, C.J. van Duijn, and S. M. Hassanizadeh.  A model for diffusion of water into a swelling particle with a free boundary: Application to a super absorbent polymer particle.  Chemical Engineering Science, 172:407–413, 2017.
  * [26] V. Thomée.  Galerkin Finite Element Methods for Parabolic Problems, volume 1054\.  Springer, 1984.
  * [27] G. B. van den Berg and C. A. Smolders.  Diffusional phenomena in membrane separation processes.  Journal of Membrane Science, 73(2-3):103–118, 1992.
  * [28] J. Wilmers and S. Bargmann.  Simulation of non-classical diffusion in polymers.  Heat Mass Transfer, 50:1543–1552, 2014.


# Label Relation Graphs Enhanced Hierarchical Residual Network for
Hierarchical Multi-Granularity Classification

Jingzhou Chen
Zhejiang University
Peng Wang
Zhejiang University
Jian Liu
Ant Group
Yuntao Qian
Zhejiang University Corresponding author

###### Abstract

Hierarchical multi-granularity classification (HMC) assigns hierarchical
multi-granularity labels to each object and focuses on encoding the label
hierarchy, _e.g_., [“Albatross”, “Laysan Albatross”] from coarse-to-fine
levels. However, the definition of what is fine-grained is subjective, and the
image quality may affect the identification. Thus, samples could be observed
at any level of the hierarchy, _e.g_., [“Albatross”] or [“Albatross”, “Laysan
Albatross”], and examples discerned at coarse categories are often neglected
in the conventional setting of HMC. In this paper, we study the HMC problem in
which objects are labeled at any level of the hierarchy. The essential designs
of the proposed method are derived from two motivations: (1) learning with
objects labeled at various levels should transfer hierarchical knowledge
between levels; (2) lower-level classes should inherit attributes related to
upper-level superclasses. The proposed combinatorial loss maximizes the
marginal probability of the observed ground truth label by aggregating
information from related labels defined in the tree hierarchy. If the observed
label is at the leaf level, the combinatorial loss further imposes the multi-
class cross-entropy loss to increase the weight of fine-grained classification
loss. Considering the hierarchical feature interaction, we propose a
hierarchical residual network (HRN), in which granularity-specific features
from parent levels acting as residual connections are added to features of
children levels. Experiments on three commonly used datasets demonstrate the
effectiveness of our approach compared to the state-of-the-art HMC approaches
and fine-grained visual classification (FGVC) methods exploiting the label
hierarchy.

##  1 Introduction

Traditional single-granularity classification usually assigns a single label
to a given object from a set of mutually exclusive class labels. For instance,
FGVC aims at distinguishing objects from different subordinate-level
categories within a given object category, e.g., subcategories of birds [32],
cars [16], aircraft [20]. However, the definition of what is fine-grained is
subjective, and the image quality may affect the identification, as
illustrated in Fig. 1. A bird can be discerned as Albatross or Laysan
Albatross due to differences in domain knowledge. Moreover, a bird expert
recognizes a bird as Albatross rather than Black-footed Albatross because of
the occlusion of key parts. Airborne or satellite image resolutions often have
large variations, causing objects to be recognized at different levels. These
challenges increase the difficulty of constructing a dataset for single-
granularity classification, while images annotated as coarse categories are
also overlooked.

(a) Differences in domain knowledge and interference from the image occlusion.

(b) Large variations of image resolutions.

Figure 1: Different objects can be discerned at various levels in the label
hierarchy due to differences in domain knowledge or image quality such as
occlusion or resolution.

Compared to single-granularity classification, a more preferable solution is
to employ hierarchical multi-granularity labels to describe an object, which
provides more flexible options for annotators with different knowledge
backgrounds [4]. HMC [15] aims to exploit hierarchical multi-granularity
labels and embeds the label hierarchy in loss function or network
architecture. Whereas conventional HMC usually evaluates each sample with
complete hierarchical labels from the coarsest to the finest granularity. A
more robust HMC model should effectively utilize examples observed at various
levels in the hierarchy, _e.g_., making use of bird images annotated as
[“Albatross”] and [“Albatross”, “Laysan Albatross”].

In this paper, we study the HMC problem in which samples are labeled at any
level of the hierarchy. We factorize this problem into two aspects: (1) how to
effectively use instances labeled at different levels; (2) how to perform
hierarchical feature interaction in the network architecture. For the first
problem, we adopt a tree hierarchy that defines two kinds of semantic
relationships between labels: parent-child correlations between levels and
mutual exclusion at the same level. Inspired by the work of [7], if an
instance is discerned at a label in the hierarchy, we maximize its marginal
probability in the probability space constrained by the tree hierarchy. Such
marginalization enjoys two benefits: learning with the coarse-level label
could impact decisions of fine-grained subclasses while learning with the
fine-level label aids the prediction of coarse-grained superclasses. Moreover,
if the ground truth label is observed at the leaf level, we further impose the
multi-class cross-entropy loss to enhance the discriminative power among fine-
grained categories.

Another critical issue is to design appropriate hierarchical feature
interaction that reflects the label hierarchy. A distinct characteristic of
hierarchical categories is that from coarse-to-fine levels, fine-level classes
not only have unique attributes but also inherit attributes related to coarse-
level superclasses. Based on this property, we propose a hierarchical residual
network (HRN) illustrated in Fig. 2. We first set up granularity-specific
layers to disentangle hierarchical features from the trunk network. Then,
these hierarchical features interact via residual connections [11, 34, 14, 13,
12, 19, 30], _i.e_., features from parent levels acting as skip connections
are added to features of children levels. Experiments on three commonly used
FGVC datasets demonstrate the effectiveness of our approach compared to the
state-of-the-art HMC approaches and FGVC methods exploiting hierarchical
knowledge under two evaluation metrics [31].

Figure 2: The network architecture consists of a trunk network (ResNet-50),
hierarchical feature interaction module, and two parallel output channels:
$O_{Hier}$ and $O_{CE}$ forming the probabilistic classification loss
($\mathcal{L}_{Hier}$) and the cross-entropy loss ($\mathcal{L}_{CE}$),
respectively. We illustrate the network architecture on CUB-200-2011 dataset
that contains three hierarchical levels. Granularity-specific block for each
hierarchical level process feature maps generated from the trunk network, then
these hierarchical features interact via residual connections, _i.e_.,
features from parent levels acting as skip connections are added to features
of children levels. $O_{Hier}$ organizes sigmoid outputs from three
hierarchical levels using the tree hierarchy, and $O_{CE}$ generates softmax
outputs corresponding to the fine-grained leaf categories.

##  2 Related Work

###  2.1 Hierarchical Multi-Granularity Classification

HMC problems naturally arise in many domains, such as text categorization [17,
25, 22] and functional genomics [1, 31, 26]. In text categorization, an
increasing number of works [21, 29, 23, 5] leveraged the label hierarchy to
improve accuracy. In image classification, HMC systems have been used to
annotate medical images [8] and classify diatom images [9]. Based on deep
neural networks (DNNs), the studies usually go along two paths: mapping the
label hierarchy to network architectures [3, 2, 24, 33] or loss functions that
impose the hierarchical constraints [7, 10]. HMC with local multi-layer
perceptrons (HMC-LMLP) [3] proposed to train a chain of multi-layer perceptron
(MLP) networks, each corresponding to a hierarchical level. The input of each
MLP uses the output provided by the previously trained MLP to augment the
feature vector of the instance. This supervised incremental greedy procedure
continues until the last level of the hierarchy is reached. HMC network (HMCN)
[33] comprised multiple local outputs, with one local output layer per
hierarchical level of the class hierarchy plus a global output layer that
captures the cumulative relationships forwarded across the entire network. All
local outputs are then concatenated and pooled with the global output to
generate a final consensual prediction. HMC-LMLP and HMCN embed label
hierarchy in their network architecture. Their loss functions sum over binary
cross-entropy losses from each hierarchical level, which assumes each label is
independent of each other, causing the implicit hierarchical relations between
two semantic labels to be ignored.

Another line of HMC works encodes the label hierarchy in loss functions by
imposing the hierarchical constraints. Coherent HMC neural network (C-HMCNN)
[10] revised the binary cross-entropy loss to satisfy the parent-child
constraint. The revision ensures that no hierarchy violation happens, _i.e_.,
for any threshold, when C-HMCNN predicts a sample belonging to a class, this
sample also belongs to its parent classes. Moreover, C-HMCNN can teach the
network how to better make the prediction on the higher level classes using
the prediction results on the lower level ones. While C-HMCNN only restricts
the parent-child correlation, other kinds of semantic relations between
hierarchical labels can be constructed using graphs. Deng _et al_. [7]
formalized semantic connections between any two labels into a directed acyclic
graph (DAG). They built a modified junction tree algorithm that contains
multiple loops during message passing on the junction tree to compute the
probabilistic classification loss defined on the DAG.

###  2.2 Fine-Grained Visual Classification

Since FGVC inherently forms a hierarchy with different levels of concept
abstraction, many approaches [35, 36, 28, 6, 4] proposed to exploit the
hierarchical label structure of FGVC. Zhang _et al_. [35] generalized the
triplet loss by describing inequalities of the distance between images
belonging to the same fine-grained class, different fine-grained classes but
the same coarse class, and different coarse classes. Shi _et al_. [28]
proposed a generalized large-margin loss that not only reduces between-class
similarity and within-class variance of the learned features but also makes
the subclasses belonging to the same coarse class be more similar than those
belonging to different coarse classes in the feature space. Chen _et al_. [6]
developed a novel hierarchical semantic embedding framework that incorporates
the predicted score vector of the higher level as prior knowledge to learn
finer-grained feature representation at each hierarchical level. During
training, the predicted score vector of the higher level is also employed to
regularize sub-categories prediction by using it as soft targets. Chang _et
al_. [4] leveraged level-specific classification heads to disentangle coarse-
level features with fine-grained ones and allowed fine-grained features to
participate in coarser-grained label predictions but constraining the gradient
flow to only update the parameters within each classification head. Their
method reaches the state-of-the-art results on the traditional single-label
FGVC problem. These approaches refine the feature representation related to
hierarchical levels in the feature space. They developed the loss function
based on the multi-class cross-entropy loss that implies the mutual exclusion
among classes at the same hierarchical level. However, they neglect to encode
other label relations like the parent-child correlation to transfer
hierarchical knowledge between levels using samples observed at different
levels.

##  3 Proposed Methods

###  3.1 Network Architecture

Our network architecture includes a trunk network, a hierarchical feature
interaction module, and two parallel output channels, see Fig. 2. The trunk
network is used to extract features from the input images and any common
network is applicable. Here, we adopt the ResNet-50 as the trunk network since
it is widely used for feature extraction. The hierarchical feature interaction
module contains granularity-specific block and residual connections. These
blocks share the same structure that comprises two convolutional layers and
two fully connected (FC) layers. Each block is designed to extract the
specialized feature for one hierarchical level. The residual connections first
linearly combine features of fine-level subclasses with features of coarse-
level superclasses. Accordingly, subclasses not only have unique attributes
but also inherit the attributes from their superclasses. Then, non-linear
transformation (ReLU) is applied to combined features.

We set up two output channels in our model. The first output channel is
utilized to compute the probabilistic classification loss based on the tree
hierarchy, in which each sigmoid node corresponds to a distinct label in the
hierarchy. We perform the non-linear projection by sigmoid instead of softmax
because sigmoid reflects the independent relations, whereas softmax implies
mutual exclusion. The sigmoid nodes from each hierarchical level are then
organized with the tree hierarchy to comply with the hierarchical constraints.
The second output channel computes the multi-class cross-entropy loss imposed
on the leaf level so that the mutually exclusive fine-grained classes gain
more attention during training. For simplicity, we denote the first and the
second output channel as $O_{Hier}$ and $O_{CE}$, respectively.

###  3.2 Loss Function

The proposed combinatorial loss integrates two forms of losses: the
probabilistic classification loss and the multi-class cross-entropy loss. We
first formalize the tree hierarchy to encode semantic relations between
hierarchical labels. The probabilistic classification loss defined on the tree
hierarchy aims to transfer hierarchical knowledge during training. We
empirically find that if the training samples labeled at the leaf level are
few, the probabilistic classification loss fails to well separate the fine-
grained leaf classes. One simple but feasible solution is to increase the
weight of fine-grained classification loss. Therefore, we further impose the
multi-class cross-entropy loss on the leaf categories, which obeys the
mutually exclusive constraint among fine-grained classes defined in the tree
hierarchy.

####  3.2.1 The Formalism of Tree Hierarchy

The tree hierarchy $G=(V,E_{h},E_{e})$ consists of a set of nodes
$V=\\{v_{1},\ldots,v_{n}\\}$, directed edges $E_{h}\subseteq V\times V$, and
undirected edges $E_{e}\subseteq V\times V$. Each node $v\in V$ corresponds to
a distinct class label. The number of nodes $n$ equals the number of all
labels in the hierarchy. A directed edge $(v_{i},v_{j})\in E_{h}$ is a
subsumption edge, indicating that class $i$ subsumes label $j$, _e.g_.,
Albatross is a parent or superclass of Black-footed Albatross. An undirected
edge $(v_{i},v_{j})\in E_{e}$ is an exclusion edge, denoting that classes
$v_{i}$ and $v_{j}$ are mutually exclusive, _e.g_., a bird cannot be the
Black-footed Albatross and Laysan Albatross simultaneously. Any two nodes
share a subsumption edge or an exclusion edge.

Each class label takes binary values, _i.e_., $v_{i}\in\\{0,1\\}$,
representing whether an object belongs to this class or not. Each edge then
defines a constraint on the binary values that two labels of its incident
nodes can take. An assignment of $(v_{i},v_{j})=(0,1)$ (e.g. a Black-footed
Albatross but not a Albatross) for a subsumption edge $(v_{i},v_{j})\in E_{h}$
is illegal, while $(v_{i},v_{j})=(1,1)$ (it is both Black-footed Albatross and
Laysan Albatross) is also an illegal assignment for an exclusion edge
$(v_{i},v_{j})\in E_{e}$. Defined by these local constraints of individual
edges, a legal global assignment of all labels in the hierarchy is a binary
label vector $\mathbf{y}\in\\{0,1\\}^{n}$ for an object. The set of all legal
global assignments forms the state space $S_{G}\subseteq\\{0,1\\}^{n}$ of tree
$G$. We can infer $S_{G}$ to be a matrix
$\mathbf{S}\in\mathcal{R}^{(n+1)\times n}$, where each row represents a legal
binary label vector $\mathbf{y}$. We traverse all legal assignments by
assigning each label a value of $1$, along with an assignment that is all
zeros.

####  3.2.2 Probabilistic Classification Loss

We calculate the probabilistic classification loss from $O_{Hier}$, and each
sigmoid node in $O_{Hier}$ corresponds to a class label in the tree hierarchy.
Suppose the number of sigmoid nodes is $n$, and $\mathbf{y}\in\\{0,1\\}^{n}$
is the binary label vector representing an assignment of all labels. Given an
input image $\mathbf{x}$, the joint probability of all sigmoid nodes
concerning the assignment $\mathbf{y}$ can be computed as:

$\tilde{P}(\mathbf{y}|\mathbf{x})=\prod^{n}_{i=1}\phi_{i}(\bar{x}_{i},y_{i})\prod_{i,j\in\\{1,\ldots,n\\}}\psi_{i,j}(y_{i},y_{j})$
(1)

where $\bar{x}_{i}$ is the sigmoid output of the $i$-th label node,
$\tilde{P}(\mathbf{y}|\mathbf{x})$ is the unnormalized probability, and
$\phi_{i}(\bar{x}_{i},y_{i})=e^{\bar{x}_{i}[y_{i}=1]}$.
$\psi_{i,j}(y_{i},y_{j})$ is the constraint defined in the tree hierarchy
between any two labels in $\mathbf{y}$:

$\psi_{i,j}(y_{i},y_{j})=\begin{cases}0,&\text{if violates constraints}\\\
1,&\text{otherwise}\end{cases}$ (2)

The joint probability is then normalized by
$Pr(\mathbf{y}|\mathbf{x})=\frac{\tilde{P}(\mathbf{y}|\mathbf{x})}{Z(\mathbf{x})}$,
where $Z(\mathbf{x})$ is the partition function that sums over all legal
assignments $\bar{\mathbf{y}}\in S_{G}$ in the state space of tree $G$:

$Z(\mathbf{x})=\sum_{\bar{\mathbf{y}}\in\\{0,1\\}^{n}}\prod^{n}_{i=1}\phi_{i}(\bar{x}_{i},\bar{y}_{i})\prod_{i,j\in\\{1,\ldots,n\\}}\psi_{i,j}(\bar{y}_{i},\bar{y}_{j})$
(3)

If input image $\mathbf{x}$ is observed at the $i$-th label in the tree
hierarchy, _i.e_., $y_{i}=1$, we can obtain the marginal probability
$Pr(y_{i}=1|\mathbf{x})$ of label $i$ by summing over all legal assignments
$\bar{\mathbf{y}}\in S_{G}$ that include $\bar{y}_{i}=1$:

$Pr(y_{i}=1|\mathbf{x})=\frac{1}{Z(\mathbf{x})}\sum_{\bar{\mathbf{y}}:\bar{y}_{i}=1}\prod_{i}\phi_{i}(\bar{x}_{i},\bar{y}_{i})\prod_{i,j}\psi_{i,j}(\bar{y}_{i},\bar{y}_{j})$
(4)

The marginal probability of a leaf label in tree $G$ relies on the sum of its
ancestors’ scores because all its ancestors must be 1 if the label of this
leaf node takes value 1, which enables the parents’ scores to impact the
descendants’ decisions. On the other hand, the marginal probability of a
parent label is marginalized over all possible states of its descendants,
_i.e_., aggregating the information from all its subclasses.

We propose to compute marginalization via matrix multiplication. Suppose the
network outputs $\mathbf{X}\in\mathcal{R}^{n\times k}$ from $O_{Hier}$, where
$n$ is the number of sigmoid nodes, and $k$ stands for batch size. Each column
in $\mathbf{X}$ is the output vector corresponding to a sample in the batch.
The unnormalized joint probability can be computed as
$\mathbf{J}=\exp(\mathbf{S}\mathbf{X})$, and the partition function
$\mathbf{z}$ can be calculated by summing each column of $\mathbf{J}$. To
obtain the marginal probability of the $j$-th sample labeled at $i$, we first
search for eligible rows in the $i$-th column of $\mathbf{S}$ that qualify
$\mathbf{S}[:,i]>0$, then we sum the corresponding elements in the $j$-th
column of $\mathbf{J}$, finally, we normalize the summation by dividing the
$j$-th element in $\mathbf{z}$.

In the training process, the observed label can be at any level of the
hierarchy, and we maximize the marginal likelihood of the observed ground
truth label. Given $m$ training samples
$\mathcal{D}=\\{\mathbf{x}^{(l)},\mathbf{y}^{(l)},g^{(l)}\\}$, $l=1,\ldots,m$,
where $\mathbf{y}^{(l)}$ is the complete ground truth label vector and
$g^{(l)}\in\\{1,\ldots,n\\}$ is the index of the observed label, the
probabilistic classification loss is defined as:

$\mathcal{L}_{Hier}(\mathcal{D})=-\frac{1}{m}\sum_{l}^{m}\ln(Pr(y^{(l)}_{g^{(l)}}=1|\mathbf{x}^{(l)}))$
(5)

####  3.2.3 Combinatorial Loss

The multi-class cross-entropy loss is commonly used in FGVC to separate fine-
grained categories. We add $\mathcal{L}_{CE}$ to our model to further increase
the discriminative power for fine-grained leaf classes. $\mathcal{L}_{CE}$
employs softmax outputs from $O_{CE}$, in which each node corresponds to a
fine-grained leaf label in the tree hierarchy. Softmax outputs imply mutually
exclusive relations among fine-grained classes, which is consistent with the
hierarchy constraint defined in the tree hierarchy. We combine
$\mathcal{L}_{CE}$ with $\mathcal{L}_{Hier}$ as follows:

$\mathcal{L}_{com}(\mathbf{x}^{(l)},y^{(l)}_{g^{(l)}})=\begin{cases}\mathcal{L}_{CE}+\mathcal{L}_{Hier},&\text{if
}g^{(l)}\text{ is in}\\\ &\text{leaf nodes}\\\
\mathcal{L}_{Hier},&\text{otherwise}\end{cases}$ (6)

Depending on whether $\mathbf{x}^{(l)}$ is labeled at fine-grained leaf
categories, the combined loss decides whether it needs to incorporate
$\mathcal{L}_{CE}$ or not. Finally, the total loss on $\mathcal{D}$ is:

$\mathcal{L}_{total}(\mathcal{D})=\sum_{l}\mathcal{L}_{com}(\mathbf{x}^{(l)},y^{(l)}_{g^{(l)}})$
(7)

##  4 Experiments

###  4.1 Implementaion Details

In all our experiments, we resize input images to $448\times 448$ and train
every single experiment for 200 epochs. Random horizontal flipping and random
cropping (random cropping for training and center cropping for testing) are
applied for data augmentation. We adopt ResNet-50 pre-trained on ImageNet as
our trunk network and use stochastic gradient descent (SGD) with a momentum of
0.9, weight decay of 0.0005 to optimize our model. The batch size is set to 8.
Meanwhile, the learning rates of the convolution layers and the FC layers
newly added for hierarchical interaction are initialized as 0.002 and adjusted
by the cosine annealing strategy [18]. The learning rates of the trunk layers
are maintained as 1/10 of the newly added layers. The code will be made
publicly accessible.

###  4.2 Datasets and Experimental Designs

We evaluate our proposed method on three widely used FGVC datasets, _i.e_.,
CUB-200-2011 [32], Aircraft [20], and Stanford Cars [16]. However,
CUB-200-2011 and Stanford Cars only provide one fine-grained label for each
image. To construct a taxonomy of label hierarchy for these two datasets, we
learn from the work of Chang _et al_. [4], in which they trace parent nodes in
Wikipedia pages. CUB-200-2011 is the most widely used benchmark for FGVC. It
covers 11,788 bird images that are re-organized into a three-level label
hierarchy with 13 orders, 38 families, and 200 species. Aircraft contains
10,000 images consisting of a three-level label hierarchy with 30 makers, 70
families, and 100 models. Stanford Cars contains 16,185 images of 196 classes
of cars. It is re-organized into a two-level label hierarchy with 9 car types
and 196 specific models. We do not use any bounding box/part annotations in
all our experiments and adopt the official train and test splits for
evaluation.

Besides assigning hierarchical multi-granularity labels for each image,
experimental designs simulate the aforementioned situation where samples are
observed at different levels of the hierarchy. To imitate the lack of domain
knowledge, we select 0%, 30%, 50%, 70%, and 90% samples from each fine-grained
class and relabel their last-level fine-grained classes to immediate parent
classes in the training set, respectively. Considering the impact of image
quality, we conduct another experiment by reducing the image resolution of
selected samples using the nearest-neighbor interpolation with a factor of 4
after relabeling. The extreme case 0% represents the conventional setting of
HMC or fine-grained classification that exploits the label hierarchy. Other
cases indicate that part of samples are observed at internal levels of the
tree hierarchy, and the rest owns the complete label hierarchy from the
highest level to the lowest fine-grained level. All images in the test set are
tested with the complete label hierarchy.

###  4.3 Evaluation Metrics

To reasonably evaluate the performance of HMC on FGVC datasets, we employ two
evaluation metrics. The first metric follows the convention of FGVC and uses
the overall accuracy (OA). The output of HMC models is a probability vector
for each class. Concerning the hierarchical label structure, we take the
maximum value of the output probability vector corresponding to each
hierarchical level as the predicted label and compute OA on the test set. The
second criterion commonly used in HMC literature [31, 33, 10] measures the
area under the average precision and recall curve $AU(\overline{PRC})$.
Instead of calculating the precision and recall curve (PRC) for each class,
$AU(\overline{PRC})$ computes an average PRC to evaluate the output
probability vector of all classes in the hierarchy. Specifically, for a given
threshold value, one point $(\overline{Prec},\overline{Rec})$ in the average
PRC is computed as:

$\displaystyle\overline{Prec}$
$\displaystyle=\frac{\sum^{n}_{i=1}TP_{i}}{\sum^{n}_{i=1}TP_{i}+\sum^{n}_{i=1}FP_{i}}$
(8) $\displaystyle\overline{Rec}$
$\displaystyle=\frac{\sum^{n}_{i=1}TP_{i}}{\sum^{n}_{i=1}TP_{i}+\sum^{n}_{i=1}FN_{i}}$

where $i$ ranges over all classes, and $TP_{i}$, $FP_{i}$, and $FN_{i}$ are
the numbers of true positives, false positives, and false negatives for class
label $i$, respectively. By varying the threshold, an average PRC is obtained
and $AU(\overline{PRC})$ denotes the area under this curve.
$AU(\overline{PRC})$ also has the advantage of being independent of the
threshold used to predict when a sample belongs to a particular class (which
is often heavily application-dependent).

###  4.4 Ablation Study

In this section, we first conduct ablation studies to investigate two key
designs of the proposed method on CUB-200-2011: HRN and combinatorial loss.
Then, we conduct visualization experiments by illustrating the attention
regions for each hierarchical level using Grad-Cam [27] to demonstrate the
role of HRN and combinatorial loss.

####  4.4.1 Significance of HRN

As displayed in Fig. 2, we analyze three components of HRN: granularity-
specific block (GSB), the linear combination of hierarchical features (LC),
and non-linear transformation of combined features (ReLU). We report OA on the
species level with the relabeling proportion of 0% in Tab. 1. The model that
only contains ResNet-50 and combinatorial loss obtained a result of 84.32. As
more components of HRN are integrated into the model, we gradually achieve
better results.

Table 1: OA on the species level with the relabeling proportion of 0% by gradually adding each component in HRN: granularity-specific block (GSB), the linear combination of hierarchical features (LC), and non-linear transformation of combined features (ReLU). Component | OA
---|---
Combinatorial Loss | 84.32
Combinatorial Loss + GSB | 85.77
Combinatorial Loss + GSB + LC | 86.17
Combinatorial Loss + GSB + LC + ReLU | 86.60

####  4.4.2 Contribution of Combinatorial Loss

In this subsection, we validate the effectiveness of combining the
probabilistic classification loss ($\mathcal{L}_{Hier}$) with the multi-class
cross-entropy loss ($\mathcal{L}_{CE}$). Tab. 2 records OA on the species
level with five relabeling proportions. In Tab. 2, it can be found that when
more training samples are relabeled to coarse-grained classes, the fine-
grained classification performance of $\mathcal{L}_{Hier}$ degenerates
drastically. In contrast, the combinatorial loss consistently outperforms
$\mathcal{L}_{Hier}$ by adding $\mathcal{L}_{CE}$ imposed on fine-grained leaf
classes.

Table 2: OA on the species level by analyzing the effectiveness of combining the probabilistic classification loss ($\mathcal{L}_{Hier}$) with the multi-class cross-entropy loss ($\mathcal{L}_{CE}$). Relabeling | 0% | 30% | 50% | 70% | 90%
---|---|---|---|---|---
$\mathcal{L}_{Hier}$ | 84.56 | 76.66 | 64.36 | 45.10 | 28.69
$\mathcal{L}_{Hier}$ \+ $\mathcal{L}_{CE}$ | 86.60 | 83.91 | 80.52 | 73.96 | 53.02

####  4.4.3 Visualization Experiments

We conduct visualization experiments to demonstrate that granularity-specific
blocks can capture different regions of interest while hierarchical knowledge
can be transferred across levels via hierarchical residual interaction and
supervision of combinatorial loss. To this end, we adopt Grad-Cam to visualize
different attention regions of each hierarchical level by propagating their
respective gradients back to feature maps generated from the trunk network.
Fig. 3 illustrates two species of birds: House Wren and Marsh Wren from the
same family (Troglodytidae) and order (Passeriformes). Domain knowledge
reveals that Marsh Wrens have a more distinctive eyebrow than House Wrens. In
Fig. 3, the order level focuses on the head and legs, and the family level
pays close attention to the head and feathers on the wing and tail. In
addition to feathers, the species level concentrates on the head. Three
hierarchical levels give similar attention to the head with a distinctive
eyebrow, while they have unique regions of interest.

(a) House Wren

(b) Marsh Wren

Figure 3: Visual attention maps of two species of birds from the same family (Troglodytidae) and order (Passeriformes). Table 3: OA(%)/$AU(\overline{PRC})$ results on CUB-200-2011 by comparing to state-of-the-art methods. Relabeling | Hierarchy |  HMC-LMLP [3] |  HMCN [33] |  C-HMCNN [10] |  Chang _et al_. [4] | Ours
---|---|---|---|---|---|---
0% | Order | 98.45 | 0.945 | 97.29 | 0.934 | 98.48 | 0.960 | 97.76 | 0.968 | 98.67 | 0.969
Family | 94.24 | 93.15 | 94.63 | 94.17 | 95.51
Specie | 79.60 | 79.75 | 81.58 | 85.56 | 86.60
30% | Order | 98.17 | 0.920 | 96.82 | 0.905 | 97.98 | 0.938 | 97.81 | 0.962 | 98.31 | 0.958
Family | 93.58 | 91.99 | 93.89 | 94.10 | 94.79
Specie | 71.30 | 71.68 | 74.91 | 82.53 | 83.91
50% | Order | 98.36 | 0.895 | 96.70 | 0.874 | 98.34 | 0.909 | 97.43 | 0.951 | 97.89 | 0.944
Family | 93.84 | 90.85 | 94.10 | 93.47 | 94.29
Specie | 64.34 | 64.29 | 67.52 | 79.30 | 80.52
70% | Order | 98.27 | 0.831 | 97.22 | 0.834 | 98.02 | 0.844 | 96.65 | 0.924 | 98.43 | 0.936
Family | 93.84 | 91.25 | 93.91 | 91.74 | 93.94
Specie | 47.98 | 52.90 | 50.05 | 70.03 | 73.96
90% | Order | 98.38 | 0.716 | 97.31 | 0.725 | 98.27 | 0.772 | 97.12 | 0.868 | 97.97 | 0.865
Family | 94.44 | 86.85 | 94.37 | 91.91 | 93.32
Specie | 22.89 | 30.69 | 26.16 | 49.36 | 53.02
Table 4: OA(%)/$AU(\overline{PRC})$ results on Aircraft by comparing to state-of-the-art methods. Relabeling | Hierarchy |  HMC-LMLP [3] |  HMCN [33] |  C-HMCNN [10] |  Chang _et al_. [4] | Ours
---|---|---|---|---|---|---
0% | Maker | 97.09 | 0.968 | 96.07 | 0.959 | 97.45 | 0.979 | 96.88 | 0.981 | 97.45 | 0.976
Family | 94.39 | 92.56 | 95.41 | 95.28 | 95.79
Model | 90.25 | 87.19 | 91.69 | 91.92 | 92.58
30% | Maker | 96.85 | 0.950 | 96.13 | 0.952 | 96.76 | 0.971 | 87.41 | 0.957 | 97.27 | 0.970
Family | 93.34 | 92.74 | 94.27 | 94.44 | 95.52
Model | 85.42 | 85.42 | 88.39 | 89.33 | 91.62
50% | Maker | 97.24 | 0.925 | 95.71 | 0.935 | 96.49 | 0.963 | 73.56 | 0.909 | 97.27 | 0.965
Family | 93.82 | 92.05 | 93.88 | 94.17 | 95.67
Model | 83.59 | 81.52 | 85.18 | 86.66 | 89.66
70% | Maker | 96.97 | 0.898 | 95.80 | 0.900 | 96.67 | 0.953 | 58.77 | 0.816 | 96.75 | 0.953
Family | 93.70 | 90.49 | 94.00 | 93.78 | 94.20
Model | 81.61 | 78.37 | 80.11 | 82.96 | 84.53
90% | Maker | 96.97 | 0.870 | 93.40 | 0.824 | 96.76 | 0.903 | 49.88 | 0.656 | 95.43 | 0.904
Family | 93.37 | 89.50 | 94.36 | 93.72 | 91.68
Model | 74.41 | 70.06 | 71.02 | 64.99 | 71.06

###  4.5 Comparison with State-of-the-art Methods

To fairly evaluate the proposed method, we compare it to state-of-the-art HMC
methods: HMC-LMLP [3], HMCN [33], and C-HMCNN [10], and the state-of-the-art
FGVC approach that exploits the label hierarchy: Chang _et al_. [4]. In our
hierarchical settings, we train all methods with different relabeling
proportions. Chang _et al_. [4] sum the multi-class cross-entropy loss from
each hierarchical level. When adapting their approach to hierarchical
settings, we neglect the last-level loss if a sample has been relabeled to its
parent class. We report OA of each hierarchical level and $AU(\overline{PRC})$
results on test sets of three FGVC datasets: CUB-200-2011, Aircraft, and
Stanford Cars, displayed in Tab. 3, Tab. 4, and Tab. 5, respectively.

From Tab. 3, Tab. 4, and Tab. 5, we can observe that the proposed method
achieves the best OA results of each hierarchical level and the best
$AU(\overline{PRC})$ results in most cases. In other cases, our results are
also comparable to the best results. Chang _et al_. [4] use level-specific
classification heads to disentangle coarse-level features with fine-grained
ones, but they only consider mutually exclusion in each hierarchical level
without examining subsumption relations between hierarchical levels in their
loss function. C-HMCNN only constrains subsumption relations. HMC-LMLP and
HMCN embed label hierarchy in their network architecture and train with the
binary cross-entropy loss that implies all classes are independent. By
contrast, in our framework, the tree hierarchy specifies the relation between
any two labels with mutually exclusion or subsumption, and the corresponding
probabilistic loss combined with the multi-class cross-entropy loss can
transfer hierarchical knowledge during training. In addition, the proposed HRN
disentangles hierarchical features by granularity-specific blocks, and these
features interact via residual connections to fuse attributes following the
hierarchy.

Table 5: OA(%)/$AU(\overline{PRC})$ results on Stanford Cars by comparing to state-of-the-art methods. Relabeling | Hierarchy |  HMC-LMLP [3] |  HMCN [33] |  C-HMCNN [10] |  Chang _et al_. [4] | Ours
---|---|---|---|---|---|---
0% | Type | 96.98 | 0.953 | 95.21 | 0.938 | 96.75 | 0.971 | 96.40 | 0.977 | 97.41 | 0.981
Maker | 87.65 | 88.71 | 90.64 | 93.65 | 94.03
30% | Type | 96.85 | 0.909 | 94.38 | 0.887 | 96.23 | 0.927 | 96.23 | 0.970 | 96.13 | 0.969
Maker | 79.16 | 81.59 | 81.92 | 91.61 | 90.55
50% | Type | 96.92 | 0.842 | 93.46 | 0.832 | 95.95 | 0.850 | 95.60 | 0.960 | 95.88 | 0.963
Maker | 66.45 | 73.03 | 70.22 | 88.10 | 88.72
70% | Type | 96.89 | 0.705 | 93.02 | 0.713 | 95.67 | 0.708 | 92.90 | 0.905 | 96.06 | 0.947
Maker | 41.52 | 52.66 | 43.17 | 76.13 | 83.72
90% | Type | 96.38 | 0.572 | 93.42 | 0.560 | 96.49 | 0.577 | 92.25 | 0.761 | 94.32 | 0.794
Maker | 13.51 | 19.89 | 13.54 | 45.79 | 49.30

###  4.6 Analyze State-of-the-art Methods by Reducing Image Resolution

Except for domain knowledge, samples captured at low-resolution can hardly be
identified with the last-level fine-grained categories, and thus they are more
likely to be inferred as upper-level coarse classes. Considering the practical
limitation of image quality, we reduce the image resolution of selected
samples corresponding to different relabeling proportions. Tab. 6 displays the
experimental results, and our method consistently outperforms compared methods
in most cases under two evaluation metrics.

Table 6: Compared OA(%)/$AU(\overline{PRC})$ results on CUB-200-2011 by reducing the image resolution after relabeling. Relabeling | Hierarchy |  HMC-LMLP [3] |  HMCN [33] |  C-HMCNN [10] |  Chang _et al_. [4] | Ours
---|---|---|---|---|---|---
0% | Order | 98.45 | 0.945 | 97.29 | 0.934 | 98.48 | 0.960 | 97.76 | 0.968 | 98.67 | 0.969
Family | 94.24 | 93.15 | 94.63 | 94.17 | 95.51
Specie | 79.60 | 79.75 | 81.58 | 85.56 | 86.60
30% | Order | 97.86 | 0.926 | 96.32 | 0.887 | 97.81 | 0.944 | 97.62 | 0.961 | 98.50 | 0.959
Family | 93.18 | 88.06 | 93.48 | 93.59 | 94.75
Specie | 74.32 | 70.78 | 76.04 | 82.33 | 84.13
50% | Order | 97.45 | 0.907 | 95.32 | 0.853 | 97.62 | 0.925 | 97.12 | 0.947 | 98.20 | 0.952
Family | 92.25 | 85.93 | 92.51 | 91.79 | 93.82
Specie | 68.10 | 62.70 | 70.37 | 78.30 | 81.18
70% | Order | 97.62 | 0.862 | 94.43 | 0.789 | 97.20 | 0.881 | 96.32 | 0.909 | 97.58 | 0.926
Family | 91.72 | 82.64 | 91.18 | 88.84 | 92.42
Specie | 53.11 | 45.75 | 55.14 | 68.06 | 73.98
90% | Order | 96.67 | 0.695 | 93.56 | 0.694 | 96.79 | 0.801 | 96.06 | 0.843 | 96.15 | 0.837
Family | 89.96 | 78.60 | 89.44 | 87.52 | 88.29
Specie | 20.78 | 22.52 | 28.32 | 46.58 | 50.10

###  4.7 Comparison with FGVC Methods Exploiting Hierarchical Knowledge

Considering hierarchical knowledge, FGVC approaches refine the feature
representation related to hierarchical levels in the feature space [35, 28, 6,
4], _e.g_., measuring the distance between classes in the hierarchy [35, 28],
learning finer-grained features with the prediction of higher level [6], or
disentangling coarse-level features with fine-grained ones [4]. Nevertheless,
they develop their loss functions based on the multi-class cross-entropy loss,
which implies mutual exclusion at the same hierarchical level. On the other
hand, encoding label relations like the parent-child correlation helps to
utilize samples observed at different levels. In contrast, the proposed method
specifies label relations with the tree hierarchy and computes the
combinatorial loss to effectively exploit samples labeled at different levels.

We record the best results reported in their works in which each sample has
complete hierarchical multi-granularity labels. In Tab. 7, Chang _et al_. [4]
achieve state-of-the-art performances in the traditional single-label FGVC
problem. Our approach reaches comparable results by simply replacing ResNet-50
with ResNeXt101-32$\times$4d [34]. Other techniques that enrich the feature
representation in the context of FGVC can be applied to boost the performance,
which is beyond our scope.

Table 7: OA on the species level by comparing to other FGVC approaches that exploit the label hierarchy. Method | CUB-200-2011 | Aircraft | Stanford Cars
---|---|---|---
Zhang _et al_. [35] |  |  | 88.4
Shi _et al_. [28] | 77.0 | 84.6 | 89.5
Chen _et al_. [6] | 88.1 |  |
Chang _et al_. [4] | 89.9 | 93.6 | 95.1
Ours | 88.6 | 94.1 | 95.1

##  5 Conclusion

We study the HMC problem in which different objects can be discerned at
various levels in the label hierarchy due to the differences in domain
knowledge or image quality. To address this problem, we propose combinatorial
loss and HRN. The combinatorial loss combines the probabilistic classification
loss defined on the tree hierarchy that encodes semantic relations between any
two hierarchical labels with the multi-class entropy loss imposed on the fine-
grained leaf categories. The probabilistic classification loss can transfer
hierarchical knowledge across levels, and the multi-class entropy loss
increases the discriminative power on the leaf classes. In addition, HRN
manages to perform hierarchical feature interaction via residual connections,
_i.e_., features from parent levels acting as skip connections are added to
features of children levels. Comprehensive experiments on three commonly used
datasets demonstrated the effectiveness of the proposed method compared to
state-of-the-art HMC and FGVC methods.

## References

  * [1] Z. Barutcuoglu, R. E. Schapire, and O. G. Troyanskaya.  Hierarchical multi-label prediction of gene function.  BMC Bioinform., 22(7):830–836, 2006.
  * [2] R. Cerri, R.C. Barros, and A. C. P. L. F. de Carvalho.  Hierarchical multi-label classification using local neural networks.  J. Comput. Syst. Sci., 80(1):39–56, 2014.
  * [3] R. Cerri, R. C. Barros, A. C. P. L. F. de Carvalho, and Y. Jin.  Reduction strategies for hierarchical multi-label classification in protein function prediction.  BMC Bioinformat., 17(1):373, 2016.
  * [4] Dongliang Chang, Kaiyue Pang, Yixiao Zheng, Zhanyu Ma, Yi-Zhe Song, and Jun Guo.  Your” flamingo” is my” bird”: Fine-grained, or not.  In CVPR, pages 11476–11485, 2021.
  * [5] B. Chen, X. Huang, L. Xiao, Z. Cai, and L. Jing.  Hyperbolic interaction model for hierarchical multi-label classification.  In AAAI, volume 34, pages 7496–7503, 2020.
  * [6] Tianshui Chen, Wenxi Wu, Yuefang Gao, Le Dong, Xiaonan Luo, and Liang Lin.  Fine-grained representation learning and recognition by exploiting hierarchical semantic embedding.  In ACM MM, pages 2023–2031, 2018.
  * [7] J. Deng, N. Ding, Y. Jia, A. Frome, K. Murphy, S. Bengio, Y. Li, H. Neven, and H. Adam.  Large-scale object classification using label relation graphs.  In ECCV, pages 48–64, 2014.
  * [8] I. Dimitrovski, D. Kocev, S. Loskovska, and S. Dzeroski.  Hierarchical annotation of medical images.  Pattern Recognit., 44(10-11):2436–2449, 2011.
  * [9] I. Dimitrovski, D. Kocev, S. Loskovska, and S. Džeroski.  Hierarchical classification of diatom images using ensembles of predictive clustering trees.  Ecological Informat., 7(1):19–29, 2012.
  * [10] E. Giunchiglia and T. Lukasiewicz.  Coherent hierarchical multi-label classification networks.  In NeurIPS, volume 33, 2020.
  * [11] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Identity mappings in deep residual networks.  In ECCV, pages 630–645. Springer, 2016.
  * [12] Andrew Howard, Mark Sandler, Grace Chu, Liang-Chieh Chen, Bo Chen, Mingxing Tan, Weijun Wang, Yukun Zhu, Ruoming Pang, Vijay Vasudevan, et al.  Searching for mobilenetv3.  In ICCV, pages 1314–1324, 2019.
  * [13] Andrew G Howard, Menglong Zhu, Bo Chen, Dmitry Kalenichenko, Weijun Wang, Tobias Weyand, Marco Andreetto, and Hartwig Adam.  Mobilenets: Efficient convolutional neural networks for mobile vision applications.  arXiv preprint arXiv:1704.04861, 2017.
  * [14] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  In CVPR, pages 4700–4708, 2017.
  * [15] C. N. Silla Jr and A. A. Freitas.  A survey of hierarchical classification across different application domains.  Data Min. Knowl. Discov., 22(1-2):31–72, 2011.
  * [16] Jonathan Krause, Michael Stark, Jia Deng, and Li Fei-Fei.  3d object representations for fine-grained categorization.  In CVPR workshops, pages 554–561, 2013.
  * [17] D. D. Lewis, Y. Yang, T. G. Rose, and F. Li.  Rcv1: A new benchmark collection for text categorization research.  J. Mach. Learn. Res., 5(Apr):361–397, 2004.
  * [18] Ilya Loshchilov and Frank Hutter.  Sgdr: Stochastic gradient descent with warm restarts.  arXiv preprint arXiv:1608.03983, 2016.
  * [19] Ningning Ma, Xiangyu Zhang, Hai-Tao Zheng, and Jian Sun.  Shufflenet v2: Practical guidelines for efficient cnn architecture design.  In ECCV, pages 116–131, 2018.
  * [20] Subhransu Maji, Esa Rahtu, Juho Kannala, Matthew Blaschko, and Andrea Vedaldi.  Fine-grained visual classification of aircraft.  arXiv preprint arXiv:1306.5151, 2013.
  * [21] Y. Mao, J. Tian, J. Han, and X. Ren.  Hierarchical text classification with reinforced label assignment.  In EMNLP-IJCNLP, pages 445–455, 2019.
  * [22] A. Mayne and R. Perry.  Hierarchically classifying documents with multiple labels.  In Proc. IEEE Symp. Comput. Intell. Data Mining, pages 133–139. IEEE, 2009.
  * [23] Y. Meng, J. Shen, C. Zhang, and J. Han.  Weakly-supervised hierarchical text classification.  In AAAI, volume 33, pages 6826–6833, 2019.
  * [24] H. Peng, J. Li, Y. He, Y. Liu, M. Bao, L. Wang, Y. Song, and Q. Yang.  Large-scale hierarchical text classification with recursively regularized deep graph-cnn.  In WWW, pages 1063–1072, 2018.
  * [25] J. Rousu, C. Saunders, S. Szedmak, and J. Shawe-Taylor.  Kernel-based learning of hierarchical multilabel classification models.  J. Mach. Learn. Res., 7:1601–1626, 2006.
  * [26] L. Schietgat, C. Vens, J. Struyf, H. Blockeel, D. Kocev, and S. Džeroski.  Predicting gene function using hierarchical multi-label decision tree ensembles.  BMC Bioinformat., 11(1):2, 2010.
  * [27] Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In ICCV, pages 618–626, 2017.
  * [28] Weiwei Shi, Yihong Gong, Xiaoyu Tao, De Cheng, and Nanning Zheng.  Fine-grained image classification using modified dcnns trained by cascaded softmax and generalized large-margin losses.  NeurIPS, 30(3):683–694, 2018.
  * [29] W. Huang et al. Hierarchical multi-label text classification: An attention-based recurrent network approach.  In Proc. 28th ACM Int. Conf. Inf. Knowl. Manag., pages 1051–1060, 2019.
  * [30] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In NeurIPS, pages 5998–6008, 2017.
  * [31] C. Vens, J. Struyf, L. Schietgat, S. Džeroski, and H. Blockeel.  Decision trees for hierarchical multi-label classification.  Mach. Learn., 73(2):185–214, 2008.
  * [32] Catherine Wah, Steve Branson, Peter Welinder, Pietro Perona, and Serge Belongie.  The caltech-ucsd birds-200-2011 dataset.  2011\.
  * [33] J. Wehrmann, R. Cerri, and R.C. Barros.  Hierarchical multi-label classification networks.  In ICML, pages 5075–5084, 2018.
  * [34] Saining Xie, Ross Girshick, Piotr Dollár, Zhuowen Tu, and Kaiming He.  Aggregated residual transformations for deep neural networks.  In CVPR, pages 1492–1500, 2017.
  * [35] Xiaofan Zhang, Feng Zhou, Yuanqing Lin, and Shaoting Zhang.  Embedding label structures for fine-grained feature representation.  In CVPR, pages 1114–1123, 2016.
  * [36] Feng Zhou and Yuanqing Lin.  Fine-grained image classification by exploring bipartite-graph labels.  In CVPR, pages 1124–1133, 2016.