# Random walks and moving boundaries: Estimating the penetration of diffusants
into dense rubbers

Surendra Nepala111corresponding author, email<EMAIL_ADDRESS>, Magnus
Ögrenb,c, Yosief Wondmagegnea, Adrian Munteana
a Department of Mathematics and Computer Science, Karlstad University,
Universitetsgatan 2, Karlstad, 65188, Sweden
b School of Science and Technology, Örebro University
SE-701 82, Örebro, Sweden
c HMU Research Center, Institute of Emerging Technologies,
GR-71004, Heraklion, Greece

###### Abstract

For certain materials science scenarios arising in rubber technology, one-
dimensional moving boundary problems (MBPs) with kinetic boundary conditions
are capable of unveiling the large-time behavior of the diffusants penetration
front, giving a direct estimate on the service life of the material. In this
paper, we propose a random walk algorithm able to lead to good numerical
approximations of both the concentration profile and the location of the sharp
front. Essentially, the proposed scheme decouples the target evolution system
in two steps: (i) the ordinary differential equation corresponding to the
evaluation of the speed of the moving boundary is solved via an explicit Euler
method, and (ii) the associated diffusion problem is solved by a random walk
method. To verify the correctness of our random walk algorithm we compare the
resulting approximations to results based on a finite element approach with a
controlled convergence rate. Our numerical experiments recover well
penetration depth measurements of an experimental setup targeting dense
rubbers.

Keywords: Moving boundary problem with a kinetic condition, explicit Euler
method, random walk approximation, finite element approximation
MSC 2020 Classification: 65M75, 65M60, 35R37

PACS: 81.05.Lg, 74.25.Ha, 05.40.Fb, 02.70.Dh

##  1 Introduction

Rubber is an intensively used material in nowadays technological solutions,
e.g., in the context of energy harvesting. Think for instance of offshore
windmills where rubber joints and ligaments are exposed to diffusants (e.g.
chlorides) able to penetrate the porous structure of the material. Even in a
dense form, rubbers are to some extent permeable and many factors affect their
permeability. For rubber, permeation is the rate at which small molecules of a
gas or liquid transfer through a rubber compound. Although these rates are
typically very low, they are important when designing a sealing process. As a
general rule, the higher the permeation, the more affected is the service life
of the material; see e.g. [3, 19, 26].

In this framework, we study a random walk method (RWM) to the solution of a
one-dimensional moving-boundary problem with a kinetic boundary condition that
describes the penetration of a population of diffusants into a dense rubber.
The target problem is well-posed (cf. [10, 11]), with the solution following
foreseeable paths for large times (cf. [1]). Furthermore, the approximation by
the finite element method (FEM) of its unique weak solution is well understood
(cf. [14, 15]). We did choose to explore the random walk approximation route
as an alternative to FEM driven by modeling reasons: It could well be that in
many practical situations even if the rubber sample is exposed to an infinite
reservoir of diffusants, these diffusants might find it difficult to penetrate
the material and hence only a tiny fraction of them succeed simply to ingress
rubber. Having in mind such situations, one may inquire why a continuum model
is needed to describe such ingress. Should such situations appear, then the
random walk approximation becomes the actual model, otherwise approaches like
[5, 25] are viable continuum-level alternative descriptions.

The RWM proposed in this framework is rather elementary, its non-standard part
is linked to the way we treat the kinetic boundary condition. We refer the
reader to [20] for the basic idea how the method works as well as, for
instance, the works [7], [24], [16], [17], [18], [23], [4], where concrete
problems have been successfully solved by RWM and related approaches. On the
computational side of things, it is worth noting that for large and
complicated physical domains in higher dimensions, the large number of mesh
points in a FEM formulation may result in equation systems that are
impractical for today’s computers. Instead, random walk methods need only a
small RAM memory and can deliver local approximation of solutions at points of
interest in the domain. For domains in practical applications connected to
processing digital images, such as photos or CT scans, filtering of the data
may result in topological errors that hinder the meshing procedure for FEM. On
the contrary, using a random walk method on a large domain with just a few
topological errors, will generally result in small statistical errors on the
calculated observables.

The paper has the following structure: In Section 2, we describe the governing
equations of our model. We present the dimensionless form of the proposed
model equations. The discussion of the setting of the model equations and the
identification of the characteristic scales of dimensionless variables is
based on [13]. In Section 3, we construct a numerical method that combines the
random walk method with the explicit Euler method to solve the dimensionless
model equations. Section 4 contains a couple of numerical experiments and
comparisons of random walk solution to the finite element approximation and
experimental lab data. Our conclusions on the constructed numerical scheme are
listed in Section 5 together with a couple of suggestions for further
investigations.

##  2 Setting of the model equations

For a one-phase Stefan problem with a kinetic condition, we consider the
following problem setting. For a fixed given observation time
$T_{f}\in(0,\infty)$, let the interval $[0,T_{f}]$ be the time span of the
involved physical processes. Let $t\in[0,T_{f}]$ denote the time variable,
$s(t)$ the position of the moving boundary at time $t$ and $x\in[0,s(t)]$ the
space variable, while $m(t,x)$ is the concentration of diffusant placed in
position $x$ at time $t$. The diffusants concentration $m(t,x)$ acts in the
non-cylindrical parabolic domain $Q_{s}(T_{f})$ defined by

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

The problem is: Find the profile of $m(t,x)$ simultaneously with the position
of the moving boundary $x=s(t)$, where the couple $(m(t,x),s(t))$ satisfies
the following system of equations:

$\displaystyle\displaystyle\frac{\partial m}{\partial
t}-D\frac{\partial^{2}m}{\partial x^{2}}=0\;\;\;\
\text{in}\;\;\;Q_{s}(T_{f}),$ (1) $\displaystyle-D\frac{\partial m}{\partial
x}(t,0)=\beta(b(t)-\text{H}m(t,0))\;\;\;\text{for}\;\;t\in(0,T_{f}),$ (2)
$\displaystyle-D\frac{\partial m}{\partial
x}(t,s(t))=s^{\prime}(t)m(t,s(t))\;\;\;\text{for}\;\;t\in(0,T_{f}),$ (3)
$\displaystyle
s^{\prime}(t)=a_{0}(m(t,s(t))-\sigma(s(t)))\;\;\;\;\text{for}\;\;\;t\in(0,T_{f}),$
(4) $\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)<\ell.$ (6)

Here $a_{0}>0$ is a kinetic coefficient, $\beta$ is a positive constant
describing the resistance of the $x=0$ interface with respect to the ingress
of the diffusants, $D>0$ is an averaged diffusion constant, and $\text{H}>0$
is the Henry constant. Additionally, $\sigma$ is a function on $\mathbb{R}$,
$b$ is a given boundary diffusant concentration on $[0,T_{f}]$, $s_{0}$ is the
initial position of the moving boundary, while $m_{0}$ represents the
concentration of the diffusant at $t=0$. As an essential part of this type of
modeling, the function $\sigma$ incorporates in an ad hoc way eventual hyper-
elastic effects, i.e. without modeling explicitly the balance of linear
momentum for the target material.

By introducing the dimensionless quantities $z=x/x_{ref},\hskip
1.13809pt\tau=tD/x_{ref}^{2},\hskip 1.13809ptu=m/m_{ref},\hskip
1.13809pth=s/x_{ref}$, we can define the region for the dimensionless
diffusant concentration by

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

where $x_{ref}$ is a characteristic length scale, while $m_{ref}$ is a
reference value for the diffusant concentration. In dimensionless form, we can
write (1)–(6) as follows:

$\displaystyle\frac{\partial u}{\partial\tau}-\frac{\partial^{2}u}{\partial
z^{2}}=0\;\;\;\ \text{in}\;\;\;Q_{h}(T),$ (7) $\displaystyle-\frac{\partial
u}{\partial z}(\tau,0)=\text{Bi}\left(\frac{b(\tau)}{m_{ref}}-{\rm
H}u(\tau,0)\right)\;\;\;\text{for}\;\;\tau\in(0,T),$ (8)
$\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),$
(9) $\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),$
(10) $\displaystyle u(0,z)=u_{0}(z)\;\;\;\text{for}\;\;\;z\in[0,h(0)],$ (11)
$\displaystyle h(0)=h_{0}>0\;\;\text{with}\;\;0<h_{0}<h(\tau)<L,$ (12)

where $u_{0}(z):=m_{0}/m_{ref}$, $h_{0}:=s_{0}/x_{ref}$ and $L:=\ell/x_{ref}$.
It is worth mentioning that $\text{Bi}:=\beta x_{ref}/D$ is the standard mass
transfer Biot number while $A_{0}:=x_{ref}m_{ref}a_{0}/D$ is the Thiele
modulus (or the 2nd Damköhler number). We refer the reader to [13] for more
details on the derivation of dimensionless model equations and finite element
simulation results. The convergence analyses of the semi-discrete Galerkin
scheme in space and of the fully discrete Galerkin-Euler scheme for the model
equations are reported in [14, 15].

##  3 Random walk method

In this section, we study a symmetric random walk method (RWM) to solve
(7)–(12), with a suitable bias when treating (10)–(12). In this method, a
finite number of mass-carrying walkers are utilized to replace the continuous
diffusant concentration field. The random movement of the walkers ensures that
the overall behavior resembles that of macroscopic diffusion. By keeping track
of the detailed movements of the random walkers, where they start from, and
where they end up, one can understand the dynamics of the problem. We refer
the reader to [20] for more details on the theory of symmetric random walks.
Note though that our situation is somewhat special as the walkers positioned
at the location of the moving boundary are exposed to a biased random walk; we
will explain this particular feature later on in Section 3.3 and Section 3.4.

We discretize the space and time domains in the following way: Let
$N,M\in\mathbb{N}$ be given. We divide the interval $[0,L]$ into $N$
subintervals. We set $0=:z_{0}<z_{1}<\cdots<z_{N}:=L$ as discretization
points. We define $\Delta z:=z_{i+1}-z_{i}$ for $i\in\\{0,1,\cdots,N-1\\}$ to
be a uniform space mesh size. We consider a finite sequence of time nodes
$0=:\tau_{1}<\tau_{2}<\cdots<\tau_{M}:=T$ with the uniform time step
$\Delta\tau:=\tau_{j+1}-\tau_{j}$ for $j\in\\{0,1,\cdots,M-1\\}.$ We consider
walkers walking randomly along the $z$ axis starting at $z=z_{0}$. The
diffusant concentration $u$ is represented by a discrete number of walkers in
the domain. Each walker has a specified mass. For the sake of simplicity, and
without loss of generality, we assume here that all walkers have exactly the
same mass of the diffusants concentration of $1$ unit in the box of length
$\Delta z$. The rule is that each walker draws a uniformly distributed random
number and according to the drawn random number decides to go to the new
position on the lattice of discrete space nodes. When walkers reach a new
position, the new position gets an increase of $1$ unit of diffusants
concentration. More clearly, at each time step $\Delta\tau$, we consider a
unit mass of diffusants concentration that moves one step of $\Delta z$ unit
length randomly to the left or to the right equally likely, that is both with
a probability $P=1/2$.

Given the initial distribution of the walkers at $\tau=0$, we are interested
in finding the number of walkers starting at $z=z_{0}$ for $\tau>0$ and
reaching the point $z=z_{i}$ after $M$ moves. To be more precise, we denote
the number of walkers at $(\tau,z)=(\tau_{j},z_{i})$ by $N_{i}^{j}$. At the
next time step $\tau=\tau_{j+1}$, we then expect $PN_{i}^{j}$ walkers moving
to the left $z=z_{i-1}$, and the same number of walkers moving to the right
$z=z_{i+1}$. At the same time, $PN_{i-1}^{j}$ walkers move from the left to
$(\tau_{j+1},z_{i})$ and $PN_{i+1}^{j}$ walkers move from the right to
$(\tau_{j+1},z_{i})$. We can express this scenario in terms of the following
discrete balance equation:

$\displaystyle
N_{i}^{j+1}=N_{i}^{j}-PN_{i}^{j}-PN_{i}^{j}+PN_{i+1}^{j}+PN_{i-1}^{j}.$

We rearrange the terms to get

$\displaystyle\frac{N_{i}^{j+1}-N_{i}^{j}}{\Delta\tau}=d\frac{(N_{i-1}^{j}-2N_{i}^{j}+N_{i+1}^{j})}{(\Delta
z)^{2}},$ (13)

where $d:=(P/\Delta\tau)(\Delta z)^{2}$. We can see (13) as a discrete form of
the diffusion equation for $N(\tau,z)$ with the diffusion coefficient $d$.
Comparing diffusion coefficients in (13) and (7) gives the relation

$\displaystyle\Delta z:=\sqrt{2\Delta\tau}.$ (14)

We have described so far the random walk method for a setting where walkers
are free to move on an infinite domain. However, we are in fact interested in
confining their motion in between a fixed and a moving boundary as described
by our original problem posed at the continuum level. Consequently, our goal
is now to solve for $u$ and $h$ in the moving domain where the moving boundary
is driven by (10). To be able to solve (7)–(12), we first decouple the
equation for the diffusant concentration and for the position of the moving
boundary. We then solve them separately. The first step is to approximate the
solution of the ordinary differential equation (10) by using the explicit
Euler method. The second step is to solve the diffusion equation by the random
walk method. The details are presented in algorithm 1

Let $h_{j}=h(\tau_{j})$, i.e. the random walk approximation of $h$ at
$\tau_{j}$ and $u(\tau_{j},z_{i})$ be the random walk approximation of $u$ at
$\tau=\tau_{j}$ and $z=z_{i}$. Before continuing with the description of the
numerical algorithm, it is worth defining an index $k_{j}$ by

$\displaystyle k_{j}:=\left\lfloor\frac{h_{j}}{\Delta
z}\right\rfloor,\;\,\,j\in\\{0,1,\cdots,M-1\\},$ (15)

where $\left\lfloor x\right\rfloor$ rounds $x$ down towards the nearest
integer. By introducing (15), one can keep track of the position of a moving
boundary and identify the neighbour discrete space node within the domain,
which can be advantageous for implementing initial and boundary conditions. At
each time slice $\tau=\tau_{j}$, our system thus consists of a set of discrete
space nodes $z=z_{i}$ for $i\in\\{0,1,\cdots,k_{j}\\}$ in which walkers are
distributed. For the walker arrived at $z=z_{i}$ for
$i\in\\{1,2,\cdots,k_{j}-1\\}$ at $\tau=\tau_{j}$, a random step, based on the
value of $p\in\\{-1,1\\}$ drawn by each walker with equal probability, is
added to the position of all the walkers. We then update the number of walkers
at each discrete space node for the next time step by using the following
formula:

$\displaystyle u(\tau_{j+1},z_{i+p})=u(\tau_{j+1},z_{i+p})+1.$ (16)

Depending on the specified boundary condition, a set of rules for the walkers
needs to be defined when they reach the boundaries $z=z_{0}$ and
$z=z_{k_{j}}$.

###  3.1 Treatment of the initial condition

For the first time step, the initial distribution of walkers within the
boundaries is sampled from the initial condition. Otherwise, it is taken from
the distribution available from the previous time step. What concerns our
problem, the initial distribution of the walkers is determined by the initial
data $u(\tau_{0},z_{i})=u_{0}(z_{i})$ for any $i\in\\{0,1,\dots,k_{0}\\}$. To
get improved numerical results, we make use of the following scaling argument:
We multiply the number of walkers starting at all discrete points defined by
the initial value with a large integer $n$, and finally, we divide the
diffusants concentration profile by $n$. That means the total number of
walkers introduced at $\tau=0$ is given by $nu_{0}(z_{i})$ for
$i\in\\{0,1,\dots,k_{0}\\}$.

###  3.2 Treatment of Dirichlet and Robin boundary conditions at the fixed
boundary

To approximate numerically the non-homogeneous Dirichlet boundary condition on
the left boundary, we employ the following two-step procedure: The Step 1
involves removing or absorbing all walkers that reach the left boundary, which
is equivalent to a homogeneous Dirichlet condition. In Step 2, we adjust the
number of walkers at the left boundary during each time step to meet the non-
homogeneous boundary condition. The walkers that remain inside the boundary
after a time step are allowed to continue their random walk during the next
time step, while those who cross the left boundary are eliminated. The
described two-steps procedure is sometimes referred to as the Söderholm
procedure; this is a possible strategy for implementing the non-homogeneous
Dirichlet boundary condition in the random walk method. We refer the reader to
[21, 22] for more details.

To deal with a Robin boundary condition at the left boundary (an inflow
boundary condition), we introduce walkers starting at the left fixed boundary.
The number of walkers starting at the left boundary $u(\tau_{j},z_{0})$ is
computed based on (8). We use the following forward difference approximation
for approximating the space derivative of $u$,

$\displaystyle\frac{\partial u}{\partial
z}(\tau_{j},z_{0})\approx\frac{u(\tau_{j},z_{1})-u(\tau_{j},z_{0})}{\Delta
z}.$

From (8), we can write

$\displaystyle\frac{u(\tau_{j},z_{0})-u(\tau_{j},z_{1})}{\Delta
z}=\text{Bi}\left(\frac{b(\tau_{j})}{m_{0}}-\text{H}u(\tau_{j},z_{0})\right).$

We multiply the number of walkers defined by the boundary conditions with the
number $n$. We then update the number of walkers at the left boundary
according to the rule:

$\displaystyle u(\tau_{j},z_{0})=\left\lfloor\frac{n\Delta
z\text{Bi}\,b(\tau_{j})/m_{0}+u(\tau_{j},z_{1})}{1+\Delta
z\rm{Bi}\,\text{H}}\right\rceil\;\;\;\text{for}\;j\in\\{1,2,\cdots,M\\},$ (17)

where $\lfloor x\rceil$ rounds $x$ to the nearest integer.

###  3.3 Treatment of the ordinary differential equation

A key issue in solving moving boundary problems by the random walk method is
the handling of the moving boundary position and its velocity (or its speed
for the 1D case). Here, the movement of the boundary $h(\tau)$ needs to be
traced together with the approximation of the concentration of diffusants. We
discretize the ordinary differential equation defining the speed of the moving
boundary by using the explicit Euler method. We then obtain the random walk
approximation of $h$ at $\tau_{j+1}$ as

$\displaystyle h_{j+1}=h_{j}+\frac{\Delta
h_{j}}{n},\;\,\,j\in\\{0,1,2,\cdots,M-1\\},$ (18)

where the total increment of the position of the moving boundary for each time
step is defined by

$\displaystyle\frac{\Delta h_{j}}{n}:=N_{j}\left[\frac{\Delta\tau
A_{0}}{n}\left(u\left(\tau_{j},h_{j}\right)-\frac{\sigma(h_{j})}{m_{ref}}\right)\right].$
(19)

Here $N_{j}$ is the total number of walkers who arrived at the moving boundary
at the time $\tau=\tau_{j}$ and contribute to the increment of the moving
boundary. The term

$\displaystyle\frac{\Delta\tau
A_{0}}{n}\left(u\left(\tau_{j},h_{j}\right)-\frac{\sigma(h_{j})}{m_{ref}}\right),$
(20)

is the increment of the boundary for a walker at time $\tau=\tau_{j}$. In
(18), we divide by $n$ to adjust for the multiplication with the factor $n$ at
the starting points, as discussed in Section 3.1 and Section 3.2. It is worth
mentioning that $N_{j}$ is unknown at this moment. But, it will be identified
while discussing the Robin boundary condition in Section 3.4. Once we know the
value for $N_{0}$, we can compute $h_{1}$ as $\Delta h_{0}$ is computable from
given initial conditions. We then solve the diffusion equation in the fixed
domain $(0,T)\times(0,h_{1})$ by the random walk method. We compute $\Delta
h_{1}$ and $h_{2}$ and then solve for $u$. We continue the same process until
the final time $T$ is reached or until all the walkers arrived at the boundary
$L$.

The increment in the moving boundary $\Delta h_{j}$ varies over time and may
not be equal to the space mesh size $\Delta z$. Hence, at any time
$\tau_{j}=j\Delta\tau$ the current position of the moving boundary may not
coincide with a space mesh point. To get to know the distance between the
moving boundary position and its neighbouring discrete space point inside the
moving domain, we defined the index $k_{j}$ for $j\in\\{0,1,\cdots,M-1\\}$ in
(15). The index (15) defines the walker’s decision to move or not towards the
current location of the moving boundary. If the distance between the position
of boundary $h_{j}$ and the position of the walker is less than $\Delta z$,
then the walker does not have a place to move to the right. That is if the
inequality $0\leq h_{j}-z_{k_{j}}\leq\Delta z$ holds true, then the walker
moves to the left, or the boundary progresses according to (19), depending on
the boundary condition at the moving boundary. While implementing (19), we use
$u\left(\tau_{j},z_{k_{j}}\right)$ instead of $u\left(\tau_{j},h_{j}\right)$.

###  3.4 Treatment of the moving boundary condition

At the right boundary $z=h(\tau)$, we have a Robin-type boundary condition
(9). This type of boundary condition is suitable for a partial reflection
scenario; see e.g. [16]. In other words, a fraction of the walkers hitting the
boundary is adsorbed with some reaction probability $P_{b}$, while the rest of
the walkers are reflected to the left.

When the walkers succeed to arrive close to the moving boundary, i.e. at
$z=z_{k_{j}}$, then the following rule applies:

  1. Step 1:

The walker present at $z=z_{k_{j}}$ moves to the left if $p=-1$. If $p=1$,
then Step 2 and Step 3 below are followed.

  2. Step 2:

For all walkers arrived at node $z=z_{k_{j}}$ and $p=1$, we evaluate the
reaction probability $P_{b}(\tau_{j})$ at that node and compare it to a
uniformly generated random number $r$ between $(0,1)$. If $r$ is less than the
current evaluation of the reaction probability $P_{b}(\tau_{j})$, then the
walkers can potentially change the position of the moving boundary. The
advancement of the boundary continues unless the moving boundary crosses the
right-hand side neighboring the discrete space node $z=z_{k_{j}+1}$. Once the
moving boundary crosses this node, the new discrete node $z=z_{k_{j}+1}$
inside the domain is updated as a node close to the moving boundary.
Consequently, the remaining walkers at $z=z_{k_{j}}$ move either to the left
$z=z_{k_{j}-1}$ or to the right $z=z_{k_{j}+1}$ with equal probability. This
means that only a fraction of the walkers arrived at $z=z_{k_{j}}$, denoted by
$N_{j}$ in (19), contribute to the overall increase or decrease of the
boundary and the rest of them contribute to the diffusion process.

  3. Step 3:

Walkers that satisfy the condition $r\geq P_{b}(\tau_{j})$ are reflected to
the left and those that contribute to the increment of the moving boundary are
adsorbed. They stick to the boundary.

In the context of our problem, we make the following splitting: firstly, we
allow the walkers to diffuse on the left side from the point close to the
moving boundary, and secondly, the walkers are exposed to the interplay
between the kinetic condition and the diffusion process. This means that the
walkers are in principle able to ”push” the moving boundary unless
breaking/acceleration effects, due to the term $\sigma(\cdot)$, dominate the
diffusion process. When the cumulative effect of their individual
contributions leads to an overall increase and the walker crosses the
neighboring discrete space node on the right-hand side, the diffusant diffuses
further up to that neighboring discrete space node. In the third step, the
walkers decide their next movement based on the values of $r$ and
$P_{b}(\tau_{j})$.

The probability that the walkers stay at the same place rather than being
reflected depends on the reactivity of the boundary $h^{\prime}(\tau_{j})$.
Relying on the definition of reaction probability defined in [16], we define
the reaction probability $P_{b}$ for our problem by

$\displaystyle P_{b}(\tau_{j})=\sqrt{2\Delta\tau}h^{\prime}(\tau_{j}).$ (21)

Other choices for the definition of the reaction probability are possible, we
refer the reader to [6, 2] for more details on this matter. For each walker
who arrived at the boundary and obtains $p=1$, equation (19) leads to the
following formula for computing the probability

$\displaystyle
P_{b}(\tau_{j})=\frac{\sqrt{2\Delta\tau}A_{0}}{n}\left(u\left(\tau_{j},z_{k_{j}}\right)-\frac{\sigma(h_{j})}{m_{ref}}\right).$
(22)

We remind the reader that $h_{j}=h(\tau_{j})$. If $r<P_{b}(\tau_{j})$, then we
update the position of the moving boundary by adding the increment of the
moving boundary for a walker (20), defined in (19).

The number of walkers for the next time step is updated by

$\displaystyle u(\tau_{j+1},z_{k_{j}})=u(\tau_{j+1},z_{k_{j}})+1.$ (23)

If $r\geq P_{b}(\tau_{j})$, then the walker is reflected and moves to the
left. For a small diffusion coefficient $D$, $P_{b}(\tau_{j})$ can have a
large value, see the definition of $A_{0}$ under (12). Hence, a suitable time
step size $\Delta\tau$ must be chosen in order to obtain the value of
$P_{b}(\tau_{j})$ between zero and one.
Claim: assume the time step size $\Delta\tau\in(0,1)$ satisfies (14) and the
following condition

$\displaystyle\sqrt{\Delta\tau}\leq\frac{n}{\sqrt{2}A_{0}\max_{j}\left(u(\tau_{j},z_{k_{j}})-\frac{\sigma(h_{j})}{m_{ref}}\right)}.$
(24)

Then it holds that $P_{b}(\tau_{j})$ in (22) is positive and bounded above by
$1$. Moreover, the increment of the moving boundary for each walker defined in
(19) satisfies the following inequality

$\displaystyle\frac{\Delta\tau
A_{0}}{n}\max_{j}\left(u\left(\tau_{j},z_{k_{j}}\right)-\frac{\sigma(h_{j})}{m_{ref}}\right)<\Delta
z.$ (25)

The positivity of $P_{b}(\cdot)$ depends on the choice of the function
$\sigma(\cdot)$. To preserve the positivity, we have to choose $\sigma(\cdot)$
in such a way that the following inequality holds:

$\displaystyle\min_{j}\left(u(\tau_{j},z_{k_{j}})-\frac{\sigma(h_{j})}{m_{ref}}\right)>0.$

In this case, all walkers who satisfy $r<P_{b}(\tau_{j})$ push the boundary.
In addition, it also holds that $h_{j+1}>h_{j}$ for all
$j\in\\{0,1,\cdots,M-1\\}$. The upper bound for $P_{b}(\cdot)$ (i.e. by $1$)
follows from (22).
Using (24), we estimate the maximum increment of the moving boundary for each
walker by

$\displaystyle\frac{\Delta\tau
A_{0}}{n}\max_{j}\left(u\left(\tau_{j},z_{k_{j}}\right)-\frac{\sigma(h_{j})}{m_{ref}}\right)\leq\Delta\tau\frac{1}{\sqrt{2\Delta\tau}}=\frac{\Delta
z}{2}<\Delta z.$

This completes the proof of the claim. It is worth mentioning that if the
following holds

$\frac{\Delta\tau
A_{0}}{n}\left(u\left(\tau_{j},z_{k_{j}}\right)-\frac{\sigma(h_{j})}{m_{ref}}\right)\gg\Delta
z,$

for some $j$, then the boundary will move several $\Delta z$ steps from the
contribution of a walker. This leads to a poor approximation of the position
of the moving boundary. To solve this issue, we have to choose the number $n$
in such a way that (25) holds true.

$z$$\tau$$z_{1}$$z_{2}$$z_{3}$$z_{4}$$z_{5}$$z_{6}$$\tau_{1}$$\tau_{2}$$\tau_{3}$$\tau_{4}$$\tau_{5}$$\tau_{6}$$\tau_{7}$
Figure 1: Sketch of the dynamics of the walkers inside the moving domain – a
vizualization of Step 1– Step 3 in the algorithm of Section 3.4.

We present a sketch in Figure 1 to illustrate the dynamics of the walkers
inside the domain. The solid line connecting two black-black or black-white
nodes represents the movement of the walkers with an equal probability, i.e.
$P=1/2$. Walkers who arrived at the white nodes, i.e. $z=z_{k_{j}}$, follow
the moving boundary conditions described in Step 1–Step 3. The solid line
connecting two white-black nodes represents the movement of the walkers with a
probability greater than or equal to $1/2$, depending on the conditions Step 1
and Step 3. The dashed line connecting two white-white nodes represents the
movement of the walkers after the moving boundary crosses the neighbourhood
node. The dashed-dotted line connecting two white-white nodes indicates the
movement of the walkers who push the moving boundary.

Algorithm 1 Procedure to compute the random walk solution to (7)–(12).

1:Choose data $a_{0},m_{0},s_{0},D,b,\sigma$.

2:Select a uniform spatial step-size $\Delta z$ generating the nodes
$z_{i}=i\Delta z,0\leq i\leq N$, in $[0,L]$ such that $N\Delta z=L$. $N$ is
the total number of discrete space nodes in the domain $[0,L]$.

3:Choose a time step $\Delta\tau$ satisfying the condition (14).

4:Consider a partition of the temporal interval $[0,T]$ generating the nodes
$\tau_{j}=j\Delta\tau,0\leq j\leq M$, in $[0,T]$ such that $M\Delta\tau=T$.
$M$ is the total number of discrete space nodes in the domain $[0,T]$.

5:Initialize a matrix $u$ for diffusants and a vector $h$ for the moving
boundary.

6:Initialize the index for space $i$, time $j$ and moving boundary $k_{j}$.

7:while $j<M-1$ and $k_{j}<N-1$ do

8: if $j==0$ then

9: Set the given initial value for $u$, i.e., $u[0,0]=u_{0}n$.

10: else

11: Compute the number of walker at the left boundary $z=z_{0}$ by using (17).

12: end if

13: for all $i\in\\{0,1,\cdots,N\\}$ do

14: for all $k\in range(u[j,i])$ do

15: Generate a number randomly $p=-1$ or $p=+1$

16: if $i+p>0$ and $i+p<=k_{j}$ and $i<=N-1$ then

17: Update $u$ i.e., $u[j+1,i+p]=u[j+1,i+p]+1$

18: else if $i+p==k_{j}+1$ then

19: Compute $\displaystyle\frac{\Delta h_{j}}{n}$ using (19)

20: Compute $P_{b}$ using (22).

21: Generate a random number $r$ uniformly distributed in the interval
$(0,1)$.

22: if $r<P_{b}$ then

23: Update $h_{j+1}$ using (18)

24: Walkers stay at the boundary, i.e., $u[j+1,k_{j}]=u[j+1,k_{j}]+1$.

25: else

26: Walkers move one step to the left, i.e., $u[j+1,k_{j}-1]=u[j+1,k_{j}-1]+1$

27: end if

28: Update $k_{j}$, i.e., $k_{j}=\text{floor}(h_{j}/\Delta z)$.

29: end if

30: end for

31: end for

32: Update time index, i.e., $j=j+1$

33:end while

34:Transfer dimensionless form of $u$ and $h$ to dimensional form $m$ and $s$.

##  4 Numerical results

Our goal in this section is to present the simulation results obtained by our
random walk algorithm. Firstly, we solve the dimensionless model equations
with a non-homogeneous Dirichlet boundary condition at $z=0$ using the random
walk method. As a next step, we solve the problem presented in (7)–(12) and
exhibit the results of the simulation. Lastly, we solve the dimensional form
of the equation described in (1)–(6) to represent the laboratory experimental
data. An analytical solution to the problem is not available. Thus, we test
the performance of the random walk algorithm in each case by solving the same
problem using the finite element method and comparing the resulting random
walk solution with the finite element solution. To approximate with finite
elements, we first transform the moving boundary problem to a problem posed in
a fixed domain using the transformation $y=z/h(\tau)\in(0,1)$. We then solve
the transformed equations using the finite element methods. For details on the
finite element approximation and simulation of (1)–(6), we refer the reader to
our previous work [14].

###  4.1 Simulation results for non-homogeneous Dirichlet boundary condition

In this section, we present the random walk solution to (7)–(12) except for
the boundary condition at $z=0$. Instead of a Robin-type boundary condition
defined in (8), we work with a non-homogeneous Dirichlet boundary condition
defined by

$u(\tau,0)=u_{D}\;\;\;\text{for}\;\;\tau\in(0,T),$

where $u_{D}$ is a constant. The aim of treating this case first is to make
the problem simpler and understand the qualitative behaviour of the random
walk solution. The observation time is set to $T=0.0001$. We take the value
$10$ mm and $0.5$ gram/mm3 for the characteristic length scale $x_{ref}$ and
for reference diffusant concentration $m_{ref}$ respectively. The function
$\sigma(h(\tau))$ is chosen as $\frac{h(\tau)}{20}$. This linear choice for
$\sigma(h(\tau))$ is taken from our previous work [13]. With choosing
$h_{0}=0.001$, we compute the index $k_{0}:=\left\lfloor\frac{h_{0}}{\Delta
z}\right\rfloor$. Initially, we place $n$ walkers at each discrete node
$z=z_{i}$ for $i\in\\{0,1,\cdots,k_{0}\\}$. We choose $u_{D}=10$ so that at
each time step, we introduce $10n$ walkers at the left boundary $z=0$. With
our choice of parameters, the dimensionless number $A_{0}$ defined in Section
2 becomes $2500$.

In Figure 2, we present the concentration profile at time $\tau=0.00005$ for
different values of $n$ and compare them to the finite element solution. The
plots in Figure 2 show that the random walk solution approaches the finite
element solution as $n$ increases. Additionally, increasing the value of $n$
leads to a better approximation for the position of the right (moving)
boundary, see Figure 3.

Figure 2:  Concentration profile at $\tau=0.00005$ by FEM and RWM for
different values of $n$, for $n=100$ (left), and $n=1000$ (right), with
$\Delta\tau=5\times 10^{-8}$.

In Figure 3, we present the corresponding simulation results for the position
of the moving boundary for different values for $n$. As a clear trend, we
observe that the quality of the approximation improves with increasing $n$ and
decreasing $\Delta\tau$.

Figure 3: Comparison of the moving front by RWM for different values of $n$
and FEM with $\Delta\tau=5\times 10^{-8}$ (left) and $\Delta\tau=2.5\times
10^{-8}$ (right).

We calculate the total mass of the concentration $M(\tau)$ in the domain
$(0,h(\tau))$ by

$\displaystyle M(\tau)=\int_{0}^{h(\tau)}u(\tau,z)dz.$ (26)

Figure 4 shows the evolution in time of the total mass for $n=2000$ and
$\Delta\tau=5\times 10^{-8}$. Comparing the random walk solution to the finite
element solution, it can be seen that the total mass of concentration obtained
by the finite element method is reproduced by the random walk solution.

Figure 4:  Comparison of the total mass over time by RWM and FEM, with
$n=2000$, $\Delta\tau=5\times 10^{-8}$.

###  4.2 Simulation results for Robin boundary condition

In this section, we present the random walk solution to (7)–(12) and compared
it with the finite element solution. We again take the observation time
$T=0.0001$. All parameters are taken to be the same as in Section 4.1 except
for the parameters arising at the left boundary. We now take $5000$ for Bi,
$2.5$ for H and $10$ for $b(\tau)$. We choose $u_{0}=1$ so that initially, we
place $n$ walkers at each discrete node at $z=z_{i}$ for
$i\in\\{0,1,\cdots,k_{0}\\}$. For the next time step, we introduce the walkers
at the left boundary $z=0$, given by (17). It is worth noting that in Section
4.1 we prescribed an equal number of walkers at $z=0$ for each time step.
However, the number of walkers at $z=0$ now depends on time. Therefore, to get
a better approximation, it is crucial to introduce correctly the number of the
walkers at $z=0$. In Figure 5, we compare the concentration profile at the
left boundary $z=0$ obtained by RWM and FEM. The plots in Figure 5 illustrate
that with increasing $n$ the approximations by the two methods are in good
agreement.

Figure 5: Numerical approximations of the concentration profile at the left
boundary $z=0$, for $n=100$ (left), and $n=1000$ (right), with
$\Delta\tau=5\times 10^{-8}$.

In Figure 6, we show the concentration profile at $\tau=0.00005$ obtained by
RWM for different values of $n$ and compare them to the solution obtained by
FEM. Comparing the plots in Figure 6 and Figure 2, we observe that the
concentration profiles, for the non-homogeneous Dirichlet and Robin-boundary
conditions, have a similar shape.

Figure 6: Numerical approximations of the concentration profile at
$\tau=0.00005$, for $n=100$ (left), and $n=1000$ (right), with
$\Delta\tau=2.5\times 10^{-8}$.

In Figure 7, we present the moving front obtained by RWM for different values
of $n$ and compare them with the FEM solution.

Figure 7: Comparison of the moving front by RWM for different values of $n$
and FEM, with $\Delta\tau=5\times 10^{-8}$ (left), and $\Delta\tau=2.5\times
10^{-8}$ (right).

We compute the finite element solution on a time mesh size $\Delta\tau=5\times
10^{-9}$ by taking $100$ finite elements. During the random walk simulation,
we observed that as long as the time step $\Delta\tau$ satisfies the
requirement (25), we can ensure that $P_{b}(\tau_{j})<1$. However, if the
walkers at the moving boundary do not satisfy $0<P_{b}(\tau_{j})<1$, the
walkers contribute neither to the increment of the moving boundary nor to the
increment of the concentration profile at the neighborhood discrete node. If
the total number of walkers who do not satisfy $0<P_{b}(\tau_{j})<1$ is
negligibly small, then the contribution of those walkers does not affect the
results of the solution. The simulations were done on a Mac with 8 cores and
16 GB RAM, using a code written in Python.

$\Delta\tau$ $n$ |  $500$ |  $1000$ |  $1500$ |  $2000$
---|---|---|---|---
$1.0\times 10^{-7}$ |  1.1139 |  2.2317 |  3.2732 |  4.4231
$5.0\times 10^{-8}$ |  3.1012 |  6.1230 |  9.4097 |  12.2056
$2.5\times 10^{-8}$ |  8.7836 |  17.3675 |  28.3242 |  37.1905
Table 1: Computational time in minutes for different values of $n$ and
$\Delta\tau$.

The computational cost of the simulation mainly depends on the value of $n$
and on the choice of the time-step size $\Delta\tau$. The effect of time
discretization and $n$ is tested by running the simulation for different
values of $\Delta\tau$ and $n$. We list the computational time for different
values of $n$ and $\Delta\tau$ in Table 1. To visualize the dependency of
computational time over the chosen $n$, we plot the corresponding
computational time over $n$ in Figure 8. We observe that the computational
time increases linearly with increasing $n$.

Figure 8: Computational time of RWM for different values of $n$ and
$\Delta\tau$.

###  4.3 Simulation results for capturing laboratory-controlled penetration
depths

In this section, we compare the experimentally measured positions of the
moving boundary with the results of our simulations based on the RWM and FEM.

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

The typical values of our reference parameters and the dimension of the
parameters are listed in Table 2. If not stated otherwise, the numerical
results included in this section are produced using the values of the
reference parameters. The reference set of parameters for the current
simulations is not precisely the same as in our previous work [13]. Some of
them are chosen based on the experimental setup, while some of them are chosen
to capture the laboratory penetration front. We perform our simulations to
recover 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). We take $T_{f}=31$ minutes for the final time. We take
the value $10$ mm and $0.5$ gram/mm3 for the characteristic length scale
$x_{ref}$ and the reference value for the concentration profile $m_{ref}$,
respectively. With this choice of parameters, the dimensional numbers Bi and
$A_{0}$ defined in Section 2 are of order $10^{4}$ and $10^{2}$, respectively.
Initially, the diffusant is uniformly distributed within the rubber up to
$0.01$ mm. Within a short time, the diffusant enters the rubber from the left
boundary $x=0$ and diffuses further into the material. In Figure 9, we present
the approximated positions of the moving boundary computed by both random walk
and finite element methods. We compare these computed profiles of penetration
depths versus time against the experimental data reported in [13].

Figure 9: Comparison of the experimental data versus the computed position of
the moving boundary by both FEM and RWM.

In Figure 10, we show the concentration profile at $t=3$ min obtained by RWM
for different values of $n$ and compare them to the solution obtained by FEM.

Figure 10: Numerical approximations of the concentration profile at $t=3$ min,
for $n=100$ (left), and $n=500$ (right), with $\Delta t=0.0005$ min.

##  5 Conclusion and outlook

We reported a random walk method (RWM) capable to approximate, with
controllable accuracy, the FEM approximation of the weak solution to our
moving boundary problem with a kinetic condition describing the penetration of
diffusants into rubber. Since in our previous work [15], we did prove
convergence rates (in suitable norms) of the FEM approximation to the original
solution, the quality of the FEM approximation transfers to the approximation
by RWM.

Trusting [12, 8, 9], it seems to be possible to analyse rigorously the
convergence property of the proposed RWM to approximate the targeted moving
boundary problem. In [12, 8], the authors solve a reaction-diffusion equation
by a fractional step method, which acts like an operator splitting between
computing a deterministic ordinary differential equation and running a random
walk method for the diffusion equation. They prove that the expected value of
the computed solution tends to the finite difference approximation of the
(strong) solution to the reaction-diffusion equation with a suitable control
on the variance. Additionally, we refer the reader also to the slightly more
sophisticated approach mentioned in [9], which fits very well to our setting
when $\sigma(r)=0$ for all $r\geq 0$. The main advantage of [9], compared to
e.g. [12], is that in the first work one points out that, at least in 1D,
moving boundary problems with kinetic conditions may admit probabilistic
interpretations, which is in principle not a priori obvious due to the non-
dissipative feature of the kinetic condition. As further work, careful
mathematical analysis needs to be performed to shed light on an eventual
probabilistic interpretation of our problem as described in Section 2.

Last but not least, the proposed RWM is able to compute the large time
behavior of the moving boundary so that the experimental results reported in
[13] are qualitatively recovered. The attention here was focused exclusively
on the dense rubber case.

## Acknowledgments

The activity of S.N. and A.M. is financed partially by the Swedish Research
Council’s project ”Homogenization and dimension reduction of thin
heterogeneous layers”, grant nr. VR 2018-03648. A.M. also thanks the Knowledge
Foundation for the grant KK 2019-0213, which led to the formulation of this
problem setting. The authors thank Dr. N. Suciu (Cluj, Romania) for inspiring
discussions on closely related matters.

## References

  * [1] T. Aiki, K. Kumazaki, and A. Muntean.  A free boundary problem describing migration into rubbers – Quest for the large time behavior.  Z. Angew. Math. Mech. (ZAMM), 102:e202100134, 2022.
  * [2] G. Boccardo, I.M. Sokolov, and A. Paster.  An improved scheme for a Robin boundary condition in discrete-time random walk algorithms.  Journal of Computational Physics, 374:1152–1165, 2018.
  * [3] C.S.P. Borges, A. Akhavan-Safar, E.A.S. Marques, R.J.C. Carbas, C. Ueffing, P. Weißgraeber, and L.F.M. da Silva.  Effect of water ingress on the mechanical and chemical properties of polybutylene terephthalate reinforced with glass fibers.  Materials, 14(5):1261, 2021.
  * [4] M.-C. Casabán, R. Company, and L. Jódar.  Numerical difference solution of moving boundary random Stefan problems.  Mathematics and Computers in Simulation, 2022.
  * [5] S. A. Chester, C. V. Di Leo, and L. Anand.  A finite element implementation of a coupled diffusion-deformation theory for elastomeric gels.  International Journal of Solids and Structures, 52:1–18, 2015.
  * [6] R. Erban and S. J. Chapman.  Reactive boundary conditions for stochastic simulations of reaction–diffusion processes.  Physical Biology, 4(1):16 – 28, 2007.
  * [7] A. Haji-Sheikh and E. M. Sparrow.  The solution of heat conduction problems by probability methods.  Journal of Heat Transfer, 89(2):121–130, 1967.
  * [8] O. H. Hald.  Convergence of random methods for a reaction-diffusion equation.  SIAM Journal on Scientific and Statistical Computing, 2(1):85–94, 1981.
  * [9] V. Kaushansky, C. Reisinger, M. Shkolnikov, and Z. Q. Song.  Convergence of a time-stepping scheme to the free boundary in the supercooled Stefan problem.  The Annals of Applied Probability, 33(1):274–298, 2023.
  * [10] 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.
  * [11] 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.
  * [12] W. Lu.  Convergence of a random walk method for a partial differential equation.  Mathematics of Computation, 67(222):593–602, 1998.
  * [13] S. Nepal, R. Meyer, N. H. Kröger, T. Aiki, A. Muntean, Y. Wondmagegne, and U. Giese.  A moving boundary approach of capturing diffusants penetration into rubber: FEM approximation and comparison with laboratory measurements.  Kautschuk Gummi Kunststoffe, 5:61–69, 2021.
  * [14] S. Nepal, Y. Wondmagegne, and A. Muntean.  Error estimates for semi-discrete finite element approximations for a moving boundary problem capturing the penetration of diffusants into rubber.  International Journal of Numerical Analysis & Modeling, 19(1):101–125, 2022.
  * [15] S. Nepal, Y. Wondmagegne, and A. Muntean.  Analysis of a fully discrete approximation to a moving-boundary problem describing rubber exposed to diffusants.  Applied Mathematics and Computation, 442:127733, 2023.  https://doi.org/10.1016/j.amc.2022.127733.
  * [16] M. Ögren.  Local boundary conditions for NMR-relaxation in digitized porous media.  The European Physical Journal B, 87(11):1–6, 2014.
  * [17] M. Ögren.  Stochastic solutions of Stefan problems with general time-dependent boundary conditions.  https://doi.org/10.1007/978-3-031-17820-7_29 In: Malyarenko, A., Ni, Y., Rančić, M., Silvestrov, S. (eds) Stochastic Processes, Statistical Methods, and Engineering Mathematics. SPAS 2019. Springer Proceedings in Mathematics & Statistics, vol 408. Springer, Cham. (arXiv:2006.04939), 2022.
  * [18] M. Ögren, D. Jha, S. Dobberschütz, D. Müter, M. Carlsson, M. Gulliksson, S.L.S. Stipp, and H.O. Sørensen.  Numerical simulations of NMR relaxation in chalk using local Robin boundary conditions.  Journal of Magnetic Resonance, 308:106597, 2019.
  * [19] E. Rostami-Tapeh-Esmaeil, A. Vahidifar, E. Esmizadeh, and D. Rodrigue.  Chemistry, processing, properties, and applications of rubber foams.  Polymers, 13(10):1565, 2021.
  * [20] S. Salsa.  Partial Differential Equations in Action: From Modelling to Theory, volume 99.  Springer, 2016.
  * [21] M. Schwind.  Some remarks on boundary conditions for random walk – the Söderholm condition.  Scripta Materialia, 48(4):461–465, 2003.
  * [22] M. Schwind and J. Ågren.  A random walk approach to Ostwald ripening.  Acta Materialia, 49(18):3821–3828, 2001.
  * [23] N. Suciu, D. Illiano, A. Prechtel, and F. A. Radu.  Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media.  Advances in Water Resources, 152:103935, 2021.
  * [24] S. Talebi, K. Gharehbash, and H. R. Jalali.  Study on random walk and its application to solution of heat conduction equation by Monte Carlo method.  Progress in Nuclear Energy, 96:18–35, 2017.
  * [25] J. Wilmers and S. Bargmann.  A continuum mechanical model for the description of solvent induced swelling in polymeric glasses: Thermomechanics coupled with diffusion.  European Journal of Mechanics-A/Solids, 53:10–18, 2015.
  * [26] N. Yasser, A. Abdelrahman, M. Kohail, and A. Moustafa.  Experimental investigation of durability properties of rubberized concrete.  Ain Shams Engineering Journal, page 102111, 2023.  https://doi.org/10.1016/j.asej.2022.102111.