Application of finite element method to General Dynamic Equation of Aerosols - Comparison with classical numerical approximations

UEF//eRepository

DSpace https://erepo.uef.fi

Rinnakkaistallenteet Luonnontieteiden ja metsätieteiden tiedekunta

2022

Application of finite element method to General Dynamic Equation of Aerosols

- Comparison with classical numerical approximations

Salminen, Teemu

Elsevier BV

Tieteelliset aikakauslehtiartikkelit

© 2021 The Authors

CC BY http://creativecommons.org/licenses/by/4.0/

http://dx.doi.org/10.1016/j.jaerosci.2021.105902

https://erepo.uef.fi/handle/123456789/26742

Downloaded from University of Eastern Finland's eRepository

Journal of Aerosol Science 160 (2022) 105902

Available online 17 November 2021

0021-8502/© 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

Contents lists available atScienceDirect

Journal of Aerosol Science

journal homepage:www.elsevier.com/locate/jaerosci

Application of finite element method to General Dynamic Equation of Aerosols – Comparison with classical numerical approximations

Teemu Salminen

, Kari E.J. Lehtinen

, Jari P. Kaipio

, Vincent Russell

, Aku Seppänen

aDepartment of Applied Physics, University of Eastern Finland, P.O. Box 1627, FI-70211, Kuopio, Finland

bFinnish Meteorological Institute, P.O. Box 1627, FI-70211, Kuopio, Finland

cDepartment of Mathematics, University of Auckland, 38 Princes St, Auckland Central, Auckland, New Zealand

A R T I C L E I N F O

Editor: Dr. Chris Hogan Keywords:

Aerosol modeling Aerosol number distribution General dynamic equation Numerical approximation

Petrov–Galerkin finite element method

A B S T R A C T

This paper focuses on the numerical approximation of the general dynamic equation of aerosols (GDE), an integro-partial differential equation, which models the temporal evolution of aerosol number distributions. More specifically, we study the applicability of the finite element method (FEM) to numerically approximate the GDE at different conditions. In addition to applying the conventional Galerkin FEM to the GDE, we also introduce its extension, the Petrov–Galerkin FEM to improve the stability of the approximation. The FEM approximations are compared against the so-called sectional method, based on the finite difference approximation, which is commonly used for the numerical approximation of the GDE. The FEM schemes and the sectional method are validated with series of numerical simulations, where an accurate solution of GDE is available. In these simulations, we consider special cases of pure condensation and coagulation, as well as cases including both of these processes. The accurate solutions – to which numerical approximations are compared – consist of either the analytical solutions or an accurate solution to the discrete GDE, where the particle size range is discretized into intervals of the size of a monomer, and the condensation process is modeled as an interaction between aerosol particles and monomers. The results of this paper demonstrate that while in the solution of the coagulation problem, the conventional sectional method is more efficient than the FEM, in cases of the condensation equation and the condensation-dominated GDE, the FE approximations outperform the sectional method – yielding a desired level of accuracy of the solution using less discretization points and in computation time which can be even orders of magnitude lower than with the sectional method.

1. Introduction

Aerosols play a major role in many atmospheric processes (Seinfeld & Pandis,2006). Aerosol particles act as cloud condensing nuclei (CCN) indirectly affecting the albedo of the Earth (Riemer, West, Zaveri, & Easter,2009) and visibility (Byun & Schere,2006).

In addition, they scatter and absorb incoming solar radiation. Together these phenomena decrease the total radiative forcing of the Earth – a quantity describing Earth’s solar radiative energy budget (Riemer et al.,2009). Currently, however, the global climate models use only crude approximations of the aerosol dynamics; this is mainly because of the high computational demand of aerosol models when combined with global models. The inaccuracy of the aerosol modeling in the global scale leaves a significant gap in

∗ Corresponding author.

E-mail address: teemu.salminen@uef.fi(T. Salminen).

https://doi.org/10.1016/j.jaerosci.2021.105902

Received 26 August 2021; Received in revised form 16 October 2021; Accepted 18 October 2021

the knowledge on their effects: According to the climate report from the Intergovernmental Panel on Climate Change (IPCC), the magnitude of the cooling effect is highly uncertain due to the imperfect modeling of aerosol dynamics in the climate models (Boucher et al.,2013). To reduce this uncertainty, it is thus crucial to develop numerical methods that allow for modeling the behavior of aerosols – especially the processes related to their formation, growth and losses – in the atmosphere adequately without increasing the computational demand too much.

The dynamics of aerosol populations depend on several processes occurring simultaneously: condensation of gas molecules on the surface of particles, agglomeration of the aerosol particles into larger particles (coagulation), dry and wet deposition and formation of new aerosol particles from the gas phase (nucleation) (Seinfeld & Pandis,2006). The time evolution of the monodisperse aerosol number distribution can be described with an integro-partial differential equation, the General Dynamic Equation of aerosols (GDE) (Seinfeld & Pandis,2006)

𝜕𝑛(𝑣_p, 𝑡)

𝜕𝑡 = − 𝜕

𝜕𝑣_p[𝐼₀(𝑣_p)𝑛(𝑣_p, 𝑡)] +1 2∫

𝑣_p 0

𝛽(𝑣_p−𝑤_p, 𝑤_p)𝑛(𝑣_p−𝑤_p, 𝑡)𝑛(𝑤_p, 𝑡)d𝑤_p

−𝑛(𝑣_p, 𝑡)

∫

∞

0

𝛽(𝑣_p, 𝑤_p)𝑛(𝑤_p, 𝑡)d𝑤_p+𝑆(𝑣_p, 𝑡) −𝑅(𝑣_p, 𝑡)𝑛(𝑣_p, 𝑡),

(1)

where𝑡is time,𝑣_p is the particle volume, 𝑛(𝑣_p, 𝑡)is the aerosol number distribution,𝐼₀(𝑣_p)is the growth/evaporation rate and 𝛽(𝑣_p, 𝑤_p)is the coagulation kernel function. In this formulation, nucleation and possible emissions are included in the source term 𝑆(𝑣_p, 𝑡), while dry and wet deposition can be included in the removal rate𝑅(𝑣_p, 𝑡)(Seinfeld & Pandis,2006).

Previously, various methods for approximating time evolutions of aerosol number distributions numerically have been sug- gested (Barrett & Webb,1998;Cui, Fu, Liang, Cheng, & Wang,2014;Liang, Wang, & Cheng,2009). One approach is to discretize GDE into multiples of a single monomer. All phenomena can then be presented as interactions between existing aerosol particles and monomers of the gas phase (Lehtinen & Kulmala,2003). This approach is, of course, impractical if the particle size range extends to sizes that exceed the monomer size by orders of magnitude. In the sectional method, on the other hand, aerosol number interval is divided into bins representing size classes of different widths (Korhonen, Lehtinen, & Kulmala,2004;Lehtinen & Kulmala,2003). In addition to these methods, multiple numerical schemes, including the finite element method, have been applied for the GDE (Pilinis, 1990;Sandu,2006;Zhang & Wexler,2002). The FE approximation of the GDE was first proposed byPilinis(1990), and further developed in papers (Rigopoulos & Jones,2003;Sandu,2006;Sandu & Borden,2003). However, the FEM has not yet been widely used for the numerical approximation of the GDE, although it has proven useful in the approximation of a variety of other physics based partial differential equations (PDEs). Probably, the main reason for the lack of using the FEM for modeling the dynamics of aerosol particle populations is that its feasibility for this application has not yet been tested sufficiently. Especially, we are not aware of published works on comparing the FE approximation of the GDE to common numerical schemes. Another reason for the infrequent use of the FEM in modeling aerosol particle dynamics has probably been the absence of open source codes for the FEM approximation of the GDE.

In the present study, the performance of the FEM in the GDE approximation is compared with the widely used sectional method.

Particularly, we consider cases of the condensation and the coagulation equations, and cases where both processes are present.

We investigate how the accuracy and computational demand behave as functions of discretization point numbers in the different numerical schemes. Here, in addition to the conventional Galerkin FEM, an upwinding scheme, the Petrov–Galerkin finite element method (PGFEM) is introduced. Upwinding schemes are usually preferable choices for solving hyperbolic problems because they stabilize the solutions (Tezduyar & Ganjoo,1986;Westerink & Shea,1989). For this reason, we study how the PGFEM compares with the Galerkin FEM in cases of the condensation equation – which is a hyperbolic PDE – and the condensation dominated GDEs.

The Matlab code package, GDE-FEM 1.0, was developed for generating all the numerical simulations of this article and is publicly available for further development.

2. Application of FEM to GDE

In the FEM, the original integro-PDE is transformed into a variational form. Next, the unknown function of the equation is represented in some finite dimensional basis, which leads to a problem of solving a set of ordinary differential equations (ODEs) (Larson & Bengzon,2010). The ODE system is solved using some suitable numerical integration scheme. In order to apply the FEM, the GDE is written for a finite volume interval [𝑣_min, 𝑣_max]:

𝜕𝑛(𝑣_p, 𝑡)

𝜕𝑡 = − 𝜕

𝜕𝑣_p[𝐼₀(𝑣_p)𝑛(𝑣_p, 𝑡)] +1 2∫

𝑣_p−𝑣_min

𝑣_min

𝛽(𝑣_p−𝑤_p, 𝑤_p)𝑛(𝑣_p−𝑤_p, 𝑡)𝑛(𝑤_p, 𝑡)d𝑤_p

−𝑛(𝑣_p, 𝑡)

∫

𝑣_max 𝑣_min

𝛽(𝑣_p, 𝑤_p)𝑛(𝑤_p, 𝑡)d𝑤_p+𝑆(𝑣_p, 𝑡) −𝑅(𝑣_p, 𝑡)𝑛(𝑣_p, 𝑡), 𝑣_p∈ [𝑣_min, 𝑣_max], 𝑡∈ [0, 𝑡_max],

(2)

where the coagulation kernel𝛽= 0outside the size range [𝑣_min, 𝑣_max]. In addition, we write a Dirichlet boundary condition at the input boundary𝑣_minand specify the initial number distribution𝑛(𝑣_p,0):

𝑛(𝑣_min, 𝑡) =𝑛_min(𝑡),

𝑛(𝑣_p,0) =𝑛₀(𝑣_p). (3)

Although in this section, the GDE is written in terms of particle volume𝑣_p only, in Section4we also consider a test case, where the size variable is particle diameter𝑑_p. The change of variable from𝑣_p to𝑑_p in the formulation of FEM and other numerical approximations is straightforward,

2.1. Variational form and FE approximation

The variational form of the GDE (Eq.(2)) is (Sandu & Borden,2003)

∫

𝑣_max 𝑣_min

[𝜕𝑛(𝑣_p, 𝑡)

𝜕𝑡 + 𝜕

𝜕𝑣_p[𝐼₀(𝑣_p)𝑛(𝑣_p, 𝑡)] −𝑆(𝑣_p, 𝑡) +𝑅(𝑣_p, 𝑡)𝑛(𝑣_p, 𝑡)

− 1 2∫

𝑣_p−𝑣_min

𝑣_min

𝛽(𝑣_p−𝑤_p, 𝑤_p)𝑛(𝑣_p−𝑤_p, 𝑡)𝑛(𝑤_p, 𝑡)d𝑤_p

+𝑛(𝑣_p, 𝑡)

∫

𝑣_max 𝑣_min

𝛽(𝑣_p, 𝑤_p)𝑛(𝑤_p, 𝑡)d𝑤_p ]

𝜉(𝑣_p)d𝑣_p= 0,

(4)

where𝜉(𝑣_p)is a test function, which has the property that it vanishes at the integration boundaries (Larson & Bengzon,2010). The particle volume interval [𝑣_min, 𝑣_max] is divided into subintervals, called finite elements (FE). In order to find the approximative solution for the variational problem (Eq. (4)), the semi-discrete Galerkin approximation is applied for the aerosol number distribution, that is, the particle number density𝑛(𝑣_p, 𝑡)is approximated with a finite set of basis functions:

𝑛(𝑣_p, 𝑡) ≈𝑛^ℎ(𝑣_p, 𝑡) =

𝑁+1∑

𝑖=0

𝑛_𝑖(𝑡)𝜙_𝑖(𝑣_p), (5)

where𝜙_𝑖 is the𝑖th FE basis function,𝑁 is the number of complete basis functions𝜙_𝑖and𝑛_𝑖(𝑡)is the respective, time-dependent coefficient parameterizing the number distribution. Similarly, we also approximate the source term𝑆(𝑣_𝑝, 𝑡)using the basis functions, i.e.,𝑆(𝑣_𝑝, 𝑡) ≈ ∑𝑁+1

𝑖=0 𝑆_𝑖(𝑡)𝜙_𝑖(𝑣_p), where 𝑆_𝑖(𝑡) are time dependent nodewise parameters describing the source processes of the particles (Sandu & Borden,2003). For the removal term, we write an approximation𝑅(𝑣_𝑝, 𝑡)𝑛(𝑣_𝑝, 𝑡) ≈∑𝑁+1

𝑖=0 𝑅_𝑖(𝑡)𝑛_𝑖(𝑡)𝜙_𝑖(𝑣_𝑝), where 𝑅_𝑖(𝑡)denotes a nodal value of𝑅(𝑣_𝑝, 𝑡). In our study, we consider only piecewise linear basis functions, which are of the form

⎧⎪

⎪⎪

⎨⎪

⎪⎪

⎩

𝜙₀(𝑣_p) = ^{𝑣1−𝑣p}

𝑣₁−𝑣_min, 𝑣_p∈ [𝑣_min, 𝑣₁] 𝜙_𝑖(𝑣_p) =

⎧⎪

⎨⎪

⎩

𝑣_p−𝑣_𝑖−1

𝑣𝑖−𝑣_𝑖−1, 𝑣_p∈ [𝑣_𝑖−1, 𝑣_𝑖]

𝑣_𝑖+1−𝑣_p

𝑣_𝑖+1−𝑣_𝑖, 𝑣_p∈ [𝑣_𝑖, 𝑣_𝑖+1] 𝜙_𝑁+1(𝑣_p) = ^{𝑣p−𝑣𝑁}

𝑣_max−𝑣_𝑁, 𝑣_p∈ [𝑣_𝑁, 𝑣_max],

(6)

where𝑣_𝑖, 𝑖= 0,…, 𝑁+ 1are the nodes of the FE mesh (Hutton,2004). The approximation is linear inside each element when the piecewise linear basis functions are used. The node points can be distributed into the interval [𝑣_min, 𝑣_max] with any choice of spacing as the node point spacing only changes the widths of the linear basis functions𝜙_𝑖(𝑣_p). In the Galerkin FEM, the piecewise linear basis functions𝜙_𝑖are used also as test functions𝜉_𝑖. With these choices, the FE approximation of the variational form (Eq.(4)) is (Sandu & Borden,2003)

𝐴 ̄𝑛^′(𝑡) =𝐺 ̄𝑛(𝑡) +

⎡⎢

⎢⎢

⎢⎣

̄

𝑛^T(𝑡)(𝐵⁰−𝐶⁰)𝑛(𝑡)̄

̄

𝑛^T(𝑡)(𝐵¹−𝐶¹)𝑛(𝑡)̄

⋮

̄

𝑛^T(𝑡)(𝐵^𝑁+1−𝐶^𝑁+1)̄𝑛(𝑡)

⎤⎥

⎥⎥

⎥⎦

+𝐴𝑆(𝑡) −𝐴𝑅(𝑡)̄𝑛(𝑡), (7)

where𝑛̄= [𝑛₀(𝑡),…, 𝑛_𝑁+1]^Tand𝑛̄^′ is the time derivative of the vector𝑛. Matrices̄ 𝐴, 𝐺, 𝑅, 𝑆, 𝐵^𝑗, and𝐶^𝑗 are of the form (Sandu &

Borden,2003):

𝐴(𝑖, 𝑗) =

∫

𝑣_max 𝑣_min

𝜙_𝑖(𝑣_p)𝜉_𝑗(𝑣_p)d𝑣_p 𝐵^𝑗(𝑖, 𝑘) = 1

2∫

𝑣_max 𝑣_min

[

∫

𝑣_p−𝑣_min

𝑣_min

𝛽(𝑣_p−𝑤_p, 𝑤_p)𝜙_𝑖(𝑣_p−𝑤_p)𝜙_𝑘(𝑤_p)d𝑤_p ]

𝜉_𝑗(𝑣_p)d𝑣_p, 𝐶^𝑗(𝑖, 𝑘) =

∫

𝑣_max 𝑣_min

[

∫

𝑣_max 𝑣_min

𝛽(𝑣_p, 𝑤_p)𝜙_𝑘(𝑤_p)d𝑤_p ]

𝜙_𝑖(𝑣_p)𝜉_𝑗(𝑣_p)d𝑣_p, 𝐺(𝑖, 𝑗) =

∫

𝑣_max 𝑣_min

𝐼₀(𝑣_p)𝜙_𝑖(𝑣_p)𝜕𝜉_𝑗(𝑣_p)

𝜕𝑣_p d𝑣_p, 𝑅(𝑡) = diag[𝑅₀(𝑡), 𝑅₁(𝑡),…, 𝑅_𝑁(𝑡), 𝑅_𝑁+1(𝑡)], 𝑆(𝑡) = [𝑆₀(𝑡), 𝑆₁(𝑡),…, 𝑆_𝑁(𝑡), 𝑆_𝑁+1(𝑡)]^T.

(8)

The FE matrix𝐺corresponding to the condensation/evaporation term of Eq.(7)is obtained using integration by parts and applying the boundary conditions𝑛(𝑣_min, 𝑡) ≈ 0 and𝑛(𝑣_max, 𝑡) ≈ 0. The sign of 𝐼₀(𝑣_p) determines whether𝐺 describes condensation or evaporation. The coagulation term in Eq.(7)can be written with 3rd order tensors𝐵 and𝐶, and tensor product resulting in the form (Sandu & Borden,2003):

𝐴 ̄𝑛^′(𝑡) =𝐺 ̄𝑛(𝑡) + [(𝐵−𝐶) ×𝑛(𝑡)]̄ ⋅𝑛(𝑡) +̄ 𝐴𝑆(𝑡) −𝐴𝑅(𝑡)̄𝑛(𝑡). (9)

Fig. 1.Galerkin and Petrov–Galerkin test functions. The standard Galerkin basis/test function𝜙_𝑖is presented in the left. The Petrov–Galerkin test functions𝜉_𝑖 with three values of upwinding factor𝜖are presented in the right. It is notable that these basis and test functions get non-zero values in two elements and value of 1 in the node point between these elements (𝑣iin this figure).

2.2. Petrov–Galerkin FEM

In most of the numerical studies of this paper, we adopt the above described Galerkin FEM, and choose the set of test functions 𝜉_𝑗(𝑣_𝑝), 𝑗= 0,…, 𝑁+1to be equal to the set of basis functions𝜙_𝑖(𝑣_𝑝), 𝑖= 0,…, 𝑁+1. However, in cases of pure condensation, the GDE becomes a hyperbolic partial differential equation, which is often notoriously unstable (Tezduyar & Ganjoo,1986). Therefore, in the case of the condensation equation, we apply an upwinding scheme referred to as the Petrov–Galerkin FEM (PGFEM) to improve the stability of the numerical approximation. In the PGFEM, the test functions𝜉_𝑗(𝑣_𝑝), 𝑗= 0,…, 𝑁+ 1are of the form

𝜉_𝑗(𝑣) =𝜙_𝑗(𝑣) +3

2𝜖(𝜎_𝑗(𝑣) −𝜎_𝑗+1(𝑣)), (10)

where𝜖∈ [0,1]is the so-called upwinding factor and𝜎_𝑗(𝑣)is a quadratic polynomial of the form 𝜎_𝑗(𝑣) =

{4

ℎ²_𝑗(𝑣−𝑣_𝑗−1)(𝑣_𝑗−𝑣), when𝑣∈ [𝑣_𝑗−1, 𝑣_𝑗] 0, otherwise,

(11) whereℎ_𝑗=𝑣_𝑗−𝑣_𝑗−1. The Galerkin and Petrov–Galerkin test functions, corresponding to three choices of upwinding factor𝜖values, are illustrated inFig. 1.

2.3. Temporal discretization

To approximate the solution of the ODE system in Eq.(9)numerically, we apply the Crank–Nicolson (C–N) method (Crank &

Nicolson,1947). The time interval[0, 𝑡_m𝑎𝑥]is discretized into𝐾subintervals of the length𝛥𝑡, so that𝐾𝛥𝑡=𝑡_max. In the C–N method, the following approximation for the time derivative𝑛̄^′is written

̄ 𝑛_𝑘+1−𝑛̄_𝑘

𝛥𝑡 =1 2

(𝑛̄^′_𝑘+1+𝑛̄^′_𝑘)

, (12)

where𝑛̄_𝑘=𝑛(𝑘𝛥𝑡), and̄ 𝑛̄^′

𝑘= ^d_d𝑡^𝑛̄(𝑘𝛥𝑡). The derivatives𝑛̄^′

𝑘+1and𝑛̄^′

𝑘can be solved from Eq.(9). By inserting these into Eq.(12)and rearranging terms, we get

̄ 𝑛_𝑘+1−𝛥𝑡

2

(𝐴⁻¹𝐺+𝐴⁻¹[

(𝐵−𝐶) ×𝑛̄_𝑘+1]

−𝑅_𝑘+1)

̄ 𝑛_𝑘+1

=𝑛̄_𝑘+𝛥𝑡 2

(𝐴⁻¹𝐺+𝐴⁻¹[

(𝐵−𝐶) ×𝑛̄_𝑘]

−𝑅_𝑘)

̄ 𝑛_𝑘+𝛥𝑡

2

(𝑆_𝑘+𝑆_𝑘+1) .

(13)

By approximating𝑛̄_𝑘+1 ≈ 𝑛̄_𝑘 in the bracketed part of the second term in the left hand side of Eq. (13), we obtain a recursive time-evolution equation for𝑛̄_𝑘:

̄ 𝑛_𝑘+1=

[ 𝐼−𝛥𝑡

2

(𝐴⁻¹𝐺+𝐴⁻¹[

(𝐵−𝐶) ×𝑛̄_𝑘]

−𝑅_𝑘+1)]−1

⋅ [

̄ 𝑛_𝑘+𝛥𝑡

2

(𝐴⁻¹𝐺+𝐴⁻¹[

(𝐵−𝐶) ×𝑛̄_𝑘]

−𝑅_𝑘)

̄ 𝑛_𝑘+𝛥𝑡

2

(𝑆_𝑘+𝑆_𝑘+1)] .

(14)

Eq.(14)forms a recursion, which can be used for computing𝑛̄_𝑘, 𝑘= 1,…, 𝐾, after projecting the initial size density𝑛(𝑣_𝑝,0)to the selected FE basis so that𝑛(𝑣_𝑝,0) ≈∑𝑁+1

𝑖=0 𝑛_𝑖(0)𝜙(𝑣_𝑝), and denoting𝑛̄₀= [𝑛₀(0),…, 𝑛_𝑁+1(0)]^T.

The above described C–N method is a combination of the explicit and implicit Euler methods, and it has second-order convergence in time (Crank & Nicolson,1947). The implicit schemes (C–N and implicit Euler method) generally admit much larger time steps than the explicit methods. The downside is that matrix inversions are needed in the implicit schemes, which increases the computational cost, especially in the large scale problems.

3. Reference methods

In this section, we review two numerical solution methods for the GDE, which both are compared with the FE approximations in Section4. In the first reference method, referred to as thesectional method (Section 3.1), the particle number distribution is discretized by dividing particles into finite-sized particle classes, or bins, and particle dynamics are approximated as transfer of particles between bins. In this approximation, the bin widths typically differ from each other. This computational method is currently widely used in the numerical approximation of the particle number distribution dynamics.

The second reference method, referred to as thediscrete GDE(Section3.2) is also based on division into discrete particle classes.

The difference from the sectional method, however, is that here the size classes are even-sized, and their size is equal to the size of a single monomer (a gas molecule). In this approach, the phenomena occurring in the aerosol population are modeled either as interactions between particles (coagulation) or as interactions between aerosol particles and monomers (condensation) (Zhang et al.,2020). This approach where the particle dynamics are modeled in molecular level, results in accurate representation of the number distribution. The drawback of this approach is the large computational demand. In this paper, the discrete GDE is used for comparison only in one test case, where an analytical solution of the GDE is not available.

3.1. The sectional method

In the sectional model, the particle number distribution is represented as a set of size classes, bins. We calculate the changes of the aerosol number distribution caused by condensation with the 1st order difference method in each of the bins. For effects caused by coagulation, we apply the size splitting method (Lehtinen & Zachariah,2001).

3.1.1. First order difference method for the condensation equation

In this section, we describe the 1st order difference method for the condensation equation

𝜕𝑛(𝑣_p, 𝑡)

𝜕𝑡 = − 𝜕

𝜕𝑣_p[𝐼₀(𝑣_p)𝑛(𝑣_p, 𝑡)], (15)

which is obtained from the GDE (Eq.(2)) if other phenomena are assumed to be negligible. We write the 1st order difference approximation for the derivative with respect to𝑣_pin Eq.(15):

d𝑛_𝑖(𝑡)

d𝑡 ≈ 𝐼₀(𝑣_𝑖−1, 𝑡)𝑛_𝑖−1(𝑡) −𝐼₀(𝑣_𝑖, 𝑡)𝑛_𝑖(𝑡)

𝑣_𝑖−𝑣_𝑖−1 (16)

where𝑛_𝑖(𝑡)is the value of the number distribution for the bin𝑖. In this form of the sectional method, the number concentration of particles in each bin is assumed to be𝑁_𝑖=𝑛_𝑖𝛥𝑣_𝑖, where the𝛥𝑣_𝑖is the width of the bin𝑖. Eq.(16)can be expressed in terms of the number concentration of particles in each size bin𝑁_𝑖:

𝜕𝑁_𝑖(𝑡)

𝜕𝑡 = 𝛥𝑣_𝑖 𝛥𝑣_𝑖−1

1

𝑣_𝑖−𝑣_𝑖−1𝐼₀(𝑣_𝑖−1)𝑁_𝑖−1(𝑡) − 1

𝑣_𝑖−𝑣_𝑖−1𝐼₀(𝑣_𝑖)𝑁_𝑖(𝑡). (17)

3.1.2. Size Splitting method for the coagulation equation

The coagulation equation describes the rate of chance of the particle density at size𝑣_𝑝:

𝜕𝑛(𝑣_p, 𝑡)

𝜕𝑡 = 1 2∫

𝑣_p−𝑣_min

𝑣_min

𝛽(𝑣_p−𝑤_p, 𝑤_p)𝑛(𝑣_p−𝑤_p, 𝑡)𝑛(𝑤_p, 𝑡)d𝑤_p−𝑛(𝑣_p, 𝑡)

∫

𝑣_max 𝑣_min

𝛽(𝑣_p, 𝑤_p)𝑛(𝑤_p, 𝑡)d𝑤_p. (18) By dividing the particles into finite sized bins as in Section3.1.1, coagulation can be described with the discrete coagulation equation (Gelbard & Seinfeld,1979)

d𝑁_𝑖(𝑡) d𝑡 =1

2

∑𝑖−1

𝑗=1

𝛽_{𝑗,𝑖−𝑗}𝑁_𝑗(𝑡)𝑁_𝑖−𝑗(𝑡) −𝑁_𝑖(𝑡)

∑𝑁

𝑗=1

𝛽_𝑖,𝑗𝑁_𝑗(𝑡). (19)

In this formulation, the bins have to have an even spacing which is equal to the size of the smallest particle. However, this restriction can be avoided, if a slight modification is made to the coagulation formation term¹

2

∑𝑖−1

𝑗=1𝛽_{𝑗,𝑖−𝑗}𝑁_𝑗𝑁_𝑖−𝑗. For this purpose, we apply the following size-splitting operator

𝜒_𝑖𝑗𝑘=

⎧⎪

⎪⎨

⎪⎪

⎩

𝑣_𝑖+1−(𝑣_𝑘+𝑣_𝑗) 𝑣_𝑖+1−𝑣𝑖

, if𝑣_𝑖≤𝑣_𝑘+𝑣_𝑗< 𝑣_𝑖+1

(𝑣_𝑘+𝑣_𝑗)−𝑣_𝑖−1

𝑣𝑖−𝑣_𝑖−1 , if𝑣_𝑖−1≤𝑣_𝑘+𝑣_𝑗< 𝑣_𝑖 0, otherwise

(20)

for the coagulation formation term, and write the discrete coagulation equation as:

d𝑁_𝑖 d𝑡 =1

2

∑𝑁

𝑘,𝑗

𝜒_𝑖𝑗𝑘𝛽_𝑗,𝑘𝑁_𝑗(𝑡)𝑁_𝑘(𝑡) −𝑁_𝑖

∑𝑁

𝑗=1

𝛽_𝑖,𝑗𝑁_𝑗(𝑡). (21)

In this formulation, any spacing for the bins can be used as the size splitting operator enforces the conservation of the total volume of the number distribution by transferring particle volume, which would end up between two bins, into the adjacent bins (Lehtinen

& Zachariah,2001).

If we consider a case where both condensation and coagulation are affecting the number distribution, the 1st order difference (Eq.(17)) and the size splitting operator (Eq.(21)) can be applied to approximate the effect of condensation and coagulation respectively. With these approximations, the sectional method for the GDE becomes

d𝑁_𝑘(𝑡) d𝑡 =1

2

∑𝑁

𝑖,𝑗

𝜒_𝑖𝑗𝑘𝛽_𝑖,𝑗𝑁_𝑖(𝑡)𝑁_𝑗(𝑡) −𝑁_𝑘

∑∞ 𝑗=1

𝛽_𝑘,𝑗𝑁_𝑗(𝑡) + 𝛥𝑣_𝑘 𝛥𝑣_𝑘−1

1

𝑣_𝑘−𝑣_𝑘−1𝐼₀(𝑣_𝑘−1)𝑁_𝑘−1− 1

𝑣_𝑘−𝑣_𝑘−1𝐼₀(𝑣_𝑘)𝑁_𝑘. (22) 3.2. Discrete GDE

A rigorous approach to model the dynamics of the particle size distribution is to use the so-called discrete GDE, where particles are represented as integer multiples of a structural unit, typically a molecule. In this approach, condensation is considered as a part of the coagulation process — namely coagulation between particles and molecules. In this form, the GDE is a set of ordinary differential equations for the number densities of each possible particle size (Gelbard & Seinfeld,1979). The discrete GDE is usually presented in terms of the number concentration of particles of each size,𝑁_𝑘. This is obtained by multiplying the size distribution of the size class (𝑛_𝑘) with the volume of the monomer𝑣₁ (𝑁_𝑘=𝑛_𝑘𝑣₁). Hence, the discrete GDE can be written as

d𝑁_𝑘(𝑡) d𝑡 =1

2

𝑘−𝑔∑^∗

𝑗=𝑔^∗

𝛽_{𝑗,𝑘−𝑗}𝑁_𝑗(𝑡)𝑁_𝑘−𝑗(𝑡) −𝑁_𝑘(𝑡)

∑∞ 𝑗=𝑔^∗

𝛽_𝑘,𝑗𝑁_𝑗(𝑡) +𝑝_𝑘−1𝑁_𝑘−1(𝑡) − (𝑝_𝑘+𝛾_𝑘)𝑁_𝑘(𝑡) +𝛾_𝑘+1(1 +𝛿_1,𝑘)𝑁_𝑘+1(𝑡) +𝑆_𝑘−𝑅_𝑘,

(23)

where𝛽_𝑘,𝑗is a value of the discrete coagulation kernel function,𝑝_𝑘is the rate of condensation,𝛾_𝑘is the rate of evaporation,𝑆_𝑘is the source term for size class𝑘, and𝑅_𝑘is the removal rate for size class𝑘(Gelbard & Seinfeld,1979;Seinfeld & Pandis,2006).

In addition,𝑘≥𝑔^∗≥2as otherwise the collision between the monomers and the smallest particles would be counted both in the condensation term and in the coagulation formation term (Seinfeld & Pandis,2006). It is notable, that this formulation is numerically heavy when the particle size range spans over multiple orders of magnitude as all the particle size class volumes are multiples of the volume of the monomer (Zhang et al.,2020). The discrete GDE is an accurate way to present particle dynamics as every phenomenon changes the number concentration of the monomers in the each size class, and monomers are, in this case, considered to be atoms or molecules (Gelbard & Seinfeld,1979;Zhang et al.,2020).

4. Numerical simulations

In this section, numerical approximations of the GDE are evaluated with four numerical simulation studies, where the FEM/

PGFEM-based particle number distributions are compared with analytical solutions (Cases 1–3) or an accurate solution obtained by the computationally demanding discrete GDE (Case 4). For comparison, the numerical approximations for the evolution of the number distribution are also computed using a reference method, the sectional method. In Cases 1 and 2, pure condensation and coagulation processes are considered, respectively. Further, in Cases 3–4, numerical approximation of the GDE with non-negligible condensation and coagulation is investigated. The specifications of the models as well as choices of numerical approximations and parameters within them are listed inTable 1.

We calculate the error for the approximated aerosol number distribution at time instant𝑡_𝑖with the following equation 𝜖(𝑡_𝑖) = 100%⋅∫_𝑣^𝑣_min^max|𝑛_acc(𝑣_p, 𝑡_𝑖) −𝑛_num(𝑣_p, 𝑡_𝑖)|d𝑣_p

∫_𝑣^𝑣_min^max𝑛_acc(𝑣_p, 𝑡_𝑖)d𝑣_p , (24)

where𝑛_acc(𝑣_p, 𝑡_𝑖)is the accurate number distribution at the time𝑡_𝑖, and𝑛_num(𝑣_p, 𝑡_𝑖)its numerical approximation. The accurate number distribution is given by either the analytical solution of the GDE (Cases 1–3) or the molecular resolution discrete GDE (Case 4). In each case, the numerical approximations are computed with various levels of discretization, that is, using a large range of numbers of elements/bins in the FE/sectional method. For each discretization, the timestep-wise error is calculated with Eq.(24). The average of the timestep-wise errors is used to approximate the error for whole time evolution:

𝜖_evolution=

∑𝐾 𝑖=1𝜖(𝑡_𝑖)

𝐾 . (25)

Computation times for each time evolution are recorded and compared to assess the efficiency of each of the numerical approximations. Simulations were conducted with MATLAB R2018b software (MATLAB,2018). Computations are done with a laptop equipped with Intel(R) Core(TM) i7-8850H CPU @ 2.60 GHz, 2592 MHz, 6 Core(s), 12 Logical Processor(s) and 16 GB of RAM. Obviously, the absolute computation times depend on the used hardware. Further, the computation times considered here do not include the pre-computation times (i.e. FE and sectional method matrix formation) as these have to only be formed once if the coagulation kernel function and/or growth rate do not change during the time evolution. The pre-computation times for the condensation FE matrices and for the coagulation FE matrices are discussed further in Sections4.1and4.2respectively.

Table 1

Aerosol process rates, initial size distributions, discretizations, temporal evolution description and initializations for the test cases. Used values for the parameters in the equations of the table are given in text.

Property/Parameter Case 1 Case 2 Case 3 Case 4

Phenomenon/model Condensation eq. Coagulation eq. GDE GDE

Growth rate 𝐼_𝑑(𝑑p) =_𝑑^𝐴

p – 𝐼(𝑣p) =𝐴₁𝑣_p 𝐼(𝑣p) =𝑁̃₀𝑣₀(3

4𝜋

)¹₆(6𝑘𝑏𝑇 𝜌𝑝

)¹₂

⋅ (1

𝑣p+_𝑣¹

0

)¹₂( 𝑣

1 3 p+𝑣

1 3 0

)2

Coagulation kernel – 𝛽(𝑣_p, 𝑤_p) =𝛽₀ 𝛽(𝑣_p, 𝑤_p) =𝛽₀ 𝛽(𝑣_p, 𝑤_p) =(3

4𝜋

)¹₆(6𝑘𝑏𝑇 𝜌𝑝

)¹₂

⋅ (₁

𝑣p

+ ¹

𝑤p

)¹

2( 𝑣

1 3 p +𝑤

1 3 p

)2

Initial distribution √ ^𝑁0 2𝜋𝑑pln𝜎_𝑔exp

(

−^ln

2(𝑑p∕𝑑̄pg) 2 ln²𝜎𝑔

) _𝑁

0𝑣p 𝑣²₀ exp

(

−^𝑣_𝑣^p

0

) _𝑁

0𝑣p 𝑣0 exp

(

−^𝑣_𝑣^p

0

) 𝑁̃₁exp(

−(𝑣p−𝜇1 𝜎1

)2) + 𝑁̃₂exp(

−(𝑣p−𝜇2 𝜎2

)2)

Evolution time 96 h 96 h 96 h 0.5 s

Particle size variable Diameter Volume Volume Volume

Size range [10⁻⁷,…,10⁻⁵]m [10⁻⁸,…,10⁻²] μm³ [10⁻⁸,…,10⁻²] μm³ [10⁻⁹,…,10⁻⁵] μm³

Accurate solution Analytical Analytical Analytical Discrete GDE

Finite element approximations

FEM, PGFEM FEM FEM FEM, PGFEM

Time integration Implicit Euler, Crank–Nicolson

Crank–Nicolson Crank–Nicolson Crank–Nicolson

Reference method Sectional Sectional Sectional Sectional

Time step lengths 0.1 h and 1 h 0.1 h and 1 h 0.1 h and 1 h 0.1 s

Size class spacing Logarithmic Logarithmic Logarithmic Logarithmic

Number of size bins/elements

26 choices, from 25 to 5000

16 choices, from 25 to 1000

16 choices, from 25 to 1000

20 choices, from 25 to 1000

4.1. Case 1: The condensation equation

In the first test case, we only consider condensation (Eq.(15)). An analytical solution for the condensation equation is presented inSeinfeld and Pandis(2006), when the independent variable is the particle diameter𝑑_pinstead of the particle volume𝑣_p, and the growth rate is of the form𝐼_𝑑(𝑑_p) = ^𝐴

𝑑_p, where𝐴is constant. With these assumptions, Eq.(15)becomes

𝜕𝑛(𝑑_p, 𝑡)

𝜕𝑡 = − 𝜕

𝜕𝑑_p[𝐼_𝑑(𝑑_p)𝑛(𝑑_p, 𝑡)]. (26)

An analytical solution is obtained for Eq.(26)with the method of the characteristics (Seinfeld & Pandis,2006):

𝑛_cond,a(𝑑_𝑝, 𝑡) = 𝑑_𝑝

√ 𝑑²_𝑝− 2𝐴𝑡

𝑛₀(√

𝑑_𝑝²− 2𝐴𝑡)

, (27)

where𝑛₀is a function defining the initial distribution. In this test, we choose the initial distribution to be of the lognormal form:

𝑛₀(𝑑_𝑝) = 𝑁₀

√2𝜋𝑑_𝑝ln𝜎_𝑔 exp

(

−ln²(𝑑_𝑝∕𝑑̄_𝑝𝑔) 2 ln²𝜎_𝑔

)

, (28)

where𝑑̄_𝑝𝑔and𝜎_𝑔, respectively, are the mean and standard deviation of the particle diameter at the initial state. In this case, the analytical solution is of the form (Seinfeld & Pandis,2006)

𝑛_cond,a(𝑑_𝑝, 𝑡) = 𝑑_𝑝𝑁₀

√2𝜋ln𝜎_𝑔(𝑑_𝑝²− 2𝐴𝑡) exp

⎛⎜

⎜⎜

⎝

− ln(√

𝑑²_𝑝− 2𝐴𝑡∕𝑑̄_𝑝𝑔) 2 ln²𝜎_𝑔

⎞⎟

⎟⎟

⎠

. (29)

For the initial parameters, we choose𝐴= 1⋅10⁻¹⁵m²h⁻¹,𝑑̄_𝑝𝑔= 0.3⋅10⁻⁶m,𝜎_𝑔= 1.2, and𝑁₀= 180. The chosen interval for the number distribution is[10⁻⁷,10⁻⁵] mand the total evolution time is 96 h.

The time evolution for the particle number distribution is calculated using the FEM and PGFEM. For the PGFEM, the upwinding factor of𝜖=0.5 is used. For both FEM schemes, two choices for time steps are tested in the C–N method: 0.1 h and 1 h. In order to test the effect of the temporal discretization scheme, we also calculate solutions using the implicit Euler method (Butcher,2008). In

Fig. 2.Case 1: The condensation equation. (a) The analytical solution and the numerical approximations based on the PGFEM and sectional method with 200 bins/elements and the C–N method with time step 0.1 h. The leftmost peak drawn with dashed line represents the initial particle number distribution, and the other two peaks correspond to instants 48 h and 96 h. (b) The relative errors of PGFEM and sectional method as functions of the evolution time, corresponding to choices 200 bins/elements and the C–N method with time step 0.1 h. (c) Temporally averaged relative errors for the time evolutions: the implicit Euler method and the C–N method for the standard Galerkin FEM and PGFEM. Here, the time step is 1 h for both temporal discretization methods. (d) The temporally averaged relative errors for PGFEM and the sectional method vs. numbers of bins/elements, corresponding to time step lengths𝛥𝑡=1 h and 0.1 h.

all FE meshes, the node points are spaced logarithmically in the chosen size interval. Pre-computation times for the condensation FE matrix formations are less than 1 s with every discretization.

For comparison, the particle number distribution evolutions are computed using the sectional method (see Section3.1). Same tests are made as for the FEM solution (seeTable 1). However, only the C–N method is used for the temporal discretization. In the sectional method, the middle points of the size bins are spaced logarithmically. The pre-computation times for the sectional method are almost equal compared to the pre-computation times of the condensation FE matrix formations.

InFig. 2a, the solutions obtained with the PGFEM and sectional method are compared with the analytical solution. The figure shows the initial number distribution, and analytical solutions as well as the PGFEM and sectional approximations corresponding to times 48 h and 96 h of evolution. InFig. 2a, both of the numerical solutions correspond to a discretization with 200 elements/bins and the time step of 0.1 h. In this example case, the PGFEM approximates the analytical solution rather well, although at time 96 h, the peak value of the size density is slightly underestimated. Also a minor negative peak is present in the 96 h PGFEM approximation.

The sectional method, on the other hand, suffers from the significant numerical diffusion and underestimates the size density peaks heavily both at 48 h and 96 h. The relative errors of the above two numerical solutions are plotted as functions of time inFig. 2b, confirming the above observation on the difference of the accuracies of PGFEM and sectional method when the above mentioned size- and time-discretizations are used.

Fig. 2c shows the temporally averaged relative errors for the PGFEM and the Galerkin FEM as functions of number of elements, when the time step of 1 h is chosen. Here, the results corresponding to the C–N time-integration are compared with those obtained with the implicit Euler method. Based on this figure, the C–N method yields better accuracy in the both FEM schemes, especially when the number of elements is above 100. Here, differences between the accuracies of the PGFEM and the Galerkin FEM are relatively small, but in the C–N based approximation, the PGFEM leads overall to slightly more accurate solutions than the Galerkin FEM, especially, when under 200 elements is used.

Fig. 2d shows the comparison between the PGFEM and the sectional method in cases of both time steps: 0.1 h and 1 h. Again, the temporally averaged relative errors are plotted as functions of numbers of elements/bins. All these solutions correspond to the C–N method. This figure reveals that when the PGFEM approximation is used, decreasing the time step from 1 h to 0.1 h improves the approximation remarkably — indeed, the relative error drops from 2% to about 0.2%, when a time step of 0.1 h is used. In the sectional method-based evolutions, the time step has only a very small effect on the accuracy. Further, when the sectional method is used, the relative error decreases rather slowly as function of the bin number and is above 3% even when the number of bins is 5000.

In conclusion, the FE-based approximations outperform the sectional method in this example case. However, it is worth to recall that this simulation represented a special case of the GDE — pure condensation. In the following sections, the numerical solution schemes are further compared in cases of coagulation equation and the GDE where the effects of condensation and coagulation are combined.

We also tested how well the numerical approximations of the condensation equation conserve the number concentration of particles. When the element/bin number is 200, the relative error of the number concentration is 1.6% in both PGFEM and sectional approximations. With 400 elements/bins, the relative error is only 0.4%. Hence, although none of the above schemes enforce the conservation of the number concentration, they conserve it well.

In order to understand the difference between PGFEM and FEM methods, we also plot the respective approximations of the size distribution at time 96 h corresponding to four FE mesh densities, seeFig. 3. The different FE meshes consist 75, 100, 125 and 150 elements, and the time step in the Crank–Nicolson method is 1 h. Clearly, the standard Galerkin FEM approximation contains oscillations that are not present in the PGFEM approximation. When the number of elements increases, the magnitude of oscillations decreases but with 150 elements they are still significant.

Finally, we test the numerical approximations in the case of the reverse process of condensation, i.e., evaporation. For this case, the particle volume𝑣_pis used as a spatial variable instead of the particle diameter𝑑_p, and the evaporation rate is𝐼₀(𝑣_p) =𝐴_Evap𝑣_p, where𝐴_Evap= −4.3⋅10⁻⁶h⁻¹. In the sectional approximation, we choose the backward 1st order difference approximation instead of the forward 1st order difference approximation used in the case of condensation, because of the direction of the peak convection is backwards in evaporation. The aerosol number distribution evolution obtained with sectional method with 5000 bins is used as a reference for the numerical methods. The results of the evaporation case are shown inFig. 4. As expected, the outcome of this test is similar to that of Case 1 (condensation): Again, the sectional approximation leads to heavy numerical diffusion, while the PGFEM retains a high accuracy with only 200 elements.

4.2. Case 2: The coagulation equation

An analytical solution for the coagulation equation (Eq. (18)) is obtained when the coagulation kernel function is constant 𝛽(𝑣_p, 𝑤_p) =𝛽₀, and the initial number distribution is of the form

𝑛_0,coag=𝑁₀𝑣_p 𝑣²

0

exp (

−𝑣_p 𝑣₀ )

, (30)

where𝑣₀ and𝑁₀are the mean particle volume and the total number of particles at time instant𝑡= 0. The analytical solution for the coagulation equation is then (Ramabhadran, Peterson, & Seinfeld,1976)

𝑛_coag,a(𝑣_p, 𝑡) = 𝑁₀(1 −𝑇)² 𝑣₀√

𝑇 exp

(

− 𝑣_p 𝑣₀ )

sinh

⎛⎜

⎜⎝ 𝑣_p√

𝑇 𝑣₀

⎞⎟

⎟⎠

, where𝑇= 1 −𝑀₀

𝑁₀, 𝑀₀= 2𝑁₀

2 +𝛽₀𝑁₀𝑡, (31)

where𝑀₀is the total number of particles at time𝑡. As the initial parameters we choose𝛽₀ = 8⋅10⁻⁶cm³h⁻¹,𝑁₀ = 1⋅10⁴, and 𝑣₀= 2⋅10⁻⁵μm³. In addition, size interval [10⁻⁸,10⁻²]μm³is chosen. For other choices of model and tests parameters, seeTable 1.

InFig. 5a, the initial size distribution, the analytical solution, the FEM and the sectional method approximations corresponding to time instants of 48 h and 96 h are presented. The numerical solutions correspond to discretizations with 100 elements/bins, and to a time step of 0.1 h in the C–N method. In this example, the FEM slightly overestimates the effect of the coagulation as the obtained numerical approximation is slightly ahead of the analytical size distribution. The sectional method is more accurate than FEM in this case as the size splitting operator conserves the particle volume (see Section3.1.2), and thus, approximates the coagulation equation with high precision. The relative errors of the approximation methods are plotted as a function of the evolution time in Fig. 5b. The relative error is higher for the FEM than for the sectional method.

Fig. 5c shows the comparison between the FEM and sectional method approximations. Again, the accuracy of the FE solutions improves when the number of elements is increased. However, in this case, with any number of discretization points, the sectional method gives more accurate approximation of the size distribution than the FEM. In addition,Fig. 5c shows that the more accurate approximations are obtained with the smaller time step.Fig. 5d shows the average relative errors as functions of the computation

Fig. 3. Case 1, The condensation equation. The analytical solution (dashed line) compared with PGFEM (blue line) and FEM (red line) approximations in cases of four FE mesh densities. The rightmost peaks correspond to time 96 h. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

time for the aerosol number distribution time evolution. This figure reveals that the computation times for the time evolutions are larger for the FEM than for the sectional method. In any case, the FE approximations are also feasible, when enough elements are used. For example, with 500 elements, the relative error of about 1.5% is obtained. With this choice, the computation time is about 10 s.

4.3. Case 3: GDE with condensation and constant coagulation kernel

An analytical solution for the GDE, where both condensation and coagulation are included, is obtained when the growth rate is of the form𝐼(𝑣_p) =𝐴₁𝑣_p, where𝐴₁is constant, and the coagulation kernel function is constant (𝛽(𝑣_p, 𝑤_p) =𝛽₀). In the case of initial number distribution

𝑛₀(𝑣_p) =𝑁₀𝑣_p 𝑣₀ exp

(

−𝑣_p 𝑣₀ )

, (32)

the analytical solution is 𝑛_GDE,a(𝑣_p, 𝑡) = 𝑁(𝑡)²

𝑀(𝑡)√

1 −𝑁(𝑡)∕𝑁₀exp (

− 𝑁₀𝑣_p 𝑀(𝑡) )

sinh(√

1 −𝑁(𝑡)∕𝑁₀ 𝑁₀𝑣_p 𝑀(𝑡) )

, (33)