University of Southern California Dissertations and Theses

Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties

(USC Thesis Other)

Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

Modeling of Bicontinuous Nanoporous Materials and Investigation of their Mechanical and Hydraulic
Properties

by

Chang Liu

A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MATERIALS SCIENCE)

August 2022

Copyright 2022 Chang Liu
ii

Acknowledgments
First and foremost, I am extremely grateful to my supervisor, Dr. Paulo S. Branicio, for his
invaluable advice, continuous support, and patience during my Ph.D. studies over the past several years.
His immense knowledge and plentiful experience have encouraged me all the time in my academic research.
He is not only my academic mentor but also a man of virtue who gave me precious advice that I can benefit
from in my whole life.
I would like to sincerely thank my qualifying and defense committee, Drs. Priya Vashista, Aiichiro
Nakano, Andrea Hodge, and Birendra Jha, for their advice, support, and contribution to this dissertation.
I would like to thank all the instructors of all the courses that I have taken during these years. The
knowledge and skills that I have learned will enable me to pursue a successful career and life. I truly want
to thank USC, the Mork Family Department, the Graduate School, and Dr. Branicio for funding me
throughout my whole Ph.D. program, which allowed me to receive the best education at USC.
I would like to thank all the USC staff that has been guiding and assisting me for all these years,
including advisors from the Mork Family Departement, Computer Science Department, VASE, OIS, and
Graduate School.
At last, I would like to thank my wife, with whom I have a cat and a lovely daughter. I also want
to thank my parents and parents-in-law. It would be impossible for me to complete my studies without the
strong support of my family.

iii

Table of Contents
Acknowledgments ......................................................................................................................................... ii
List of Tables ................................................................................................................................................ v
List of Figures .............................................................................................................................................. vi
Abbreviations .............................................................................................................................................. vii
Abstract ...................................................................................................................................................... viii
Chapter 1: Introduction ................................................................................................................................. 1
1.1 Nanoporous Materials ......................................................................................................................... 1
1.2 Metallic Glasses .................................................................................................................................. 2
1.3 Fluid Flow in Nanoporous Media ....................................................................................................... 3
Chapter 2: Methodology ............................................................................................................................... 4
2.1 Molecular Dynamics ........................................................................................................................... 4
2.1.1 Overview of the Molecular Dynamics Method ............................................................................ 4
2.1.2 Interatomic Potential .................................................................................................................... 5
2.1.3 Velocity-Verlet Algorithm ........................................................................................................... 6
2.1.4 Periodic Boundary Conditions ..................................................................................................... 7
2.1.5 Shifted Potential ........................................................................................................................... 9
2.2 Dissipative Particle Dynamics ............................................................................................................ 9
2.2.1 Formulation of the Dissipative Particle Dynamics Method ......................................................... 9
2.2.2 Formulation of the Many-Body Dissipative Particle Dynamics Method .................................. 11
2.2.3 Time Integration ......................................................................................................................... 12
Chapter 3: Generation of Non-Cubic Stochastic Periodic Bicontinuous Nanoporous Structures .............. 14
3.1 Background ....................................................................................................................................... 14
3.2 Generation of Nanoporous Structures in Periodic Cells ................................................................... 16
3.2.1 Cubic Periodic Cells ................................................................................................................... 16
3.2.2 Arbitrary Parallelepiped Periodic Cells ..................................................................................... 17
3.3 Structural Analysis ............................................................................................................................ 22
3.3.1 Ligament Size Distribution ........................................................................................................ 22
3.3.2 Mean and Gaussian Curvatures .................................................................................................. 23
3.3.3 Nanoporous Structure Genus ..................................................................................................... 25
3.4 Extension of the Method: Gradient Nanoporous Structure ............................................................... 27
3.5 Conclusions ....................................................................................................................................... 29
Chapter 4: Bicontinuous Nanoporous Design Induced Homogenization of Strain Localization in Metallic
Glasses ........................................................................................................................................................ 31
4.1 Background ....................................................................................................................................... 31
4.2 Method .............................................................................................................................................. 32
4.2.1 Model Generation ...................................................................................................................... 32
4.2.2 Measurement of Surface Curvature ........................................................................................... 33
4.2.3 Simulation Details ...................................................................................................................... 33
4.3 Results ............................................................................................................................................... 35
4.3.1 Analysis of Bicontinuous Topologies ........................................................................................ 35
4.3.2 Mechanical Response to Tensile Loading ................................................................................. 35
4.3.3 Realignment of Ligaments Towards Loading Direction ............................................................ 39
4.3.4 Analysis of the Atomic Local Strain during Tensile Loading ................................................... 40
4.4 Discussion ......................................................................................................................................... 42
iv

4.4.1 Local Heating ............................................................................................................................. 42
4.4.2 Comparison to Nanoporous Metallic Glasses with Simple-Shaped Pores ................................ 42
4.4.3 Effect of Metallic Glasses Ligaments ........................................................................................ 43
4.4.4 Comparison to Bicontinuous Nanoporous Crystalline Metals ................................................... 44
4.4.5 Bicontinuous Nanoporous Metallic Glass vs. Collection of Metallic Glass Nanopillars .......... 45
4.4.6 Effect of High Strain Rate Used in Molecular Dynamics Simulations ...................................... 46
4.5 Conclusions ....................................................................................................................................... 47
Chapter 5: Structural and Compositional Effects on the Mechanical Properties and Scaling Laws of
Bicontinuous Nanoporous Metallic Glasses ............................................................................................... 48
5.1 Background ....................................................................................................................................... 48
5.2 Method .............................................................................................................................................. 49
5.2.1 Model Generation and Morphology Characterization. .............................................................. 49
5.2.2 Simulation Details ...................................................................................................................... 51
5.2.3 Genetic Programming ................................................................................................................ 51
5.3 Results ............................................................................................................................................... 54
5.3.1 Geometrical and morphological features of bicontinuous NPMG ............................................. 54
5.3.2 Mechanical Properties of bicontinuous CuZr NPMG ................................................................ 56
5.3.3 Relative Density Effect on Deformation and Failure Modes ..................................................... 59
5.3.4 Scaling Laws Uncovered by Genetic Programming Assisted Symbolic Regression ................ 61
5.4 Discussion ......................................................................................................................................... 63
5.4.1 System Size Effect on Mechanical Properties ........................................................................... 63
5.4.2 Composition Effect on Mechanical Properties .......................................................................... 64
5.4.3 Existing Scaling Laws with Relative Density ............................................................................ 65
5.4.4 Comparison of GPSR Uncovered and Existing Scaling Laws .................................................. 67
5.5 Conclusions ....................................................................................................................................... 68
Chapter 6: Pore-Size Dependence of Permeability in Bicontinuous Nanoporous Media ........................... 70
6.1. Background ...................................................................................................................................... 70
6.2 Method .............................................................................................................................................. 73
6.2.1 Bicontinuous Nanoporous Model and its Tortuosity ................................................................. 73
6.2.2 Configurations of Solid and Liquid Phases ................................................................................ 74
6.3 Results ............................................................................................................................................... 78
6.3.1 Self-Closed Void During Relaxation ......................................................................................... 78
6.3.2 Pressure Distribution and Variation in the System .................................................................... 79
6.3.3 Motion of Pistons ....................................................................................................................... 80
6.3.4 Local Velocities ......................................................................................................................... 82
6.4 Discussion ......................................................................................................................................... 87
6.4.1 Pressure Evolution ..................................................................................................................... 87
6.4.2 Boundary Condition at the Solid-Fluid Interface ....................................................................... 89
6.4.3 Pore Size Effect on Permeability Derived from Hagen-Poiseuille Law .................................... 89
6.5 Conclusions ....................................................................................................................................... 94
Chapter 7: Conclusions ............................................................................................................................... 95
References ................................................................................................................................................... 98

v

List of Tables
Table 1 Pseudocode of an MD simulation. ................................................................................................... 5
Table 2 Pseudocode of the Velocity-Verlet algorithm. ................................................................................. 7
Table 3 Scaling laws for Young’s modulus obtained by GP. ..................................................................... 62
Table 4 Scaling law for ultimate tensile strength obtained by GP. ............................................................. 62
Table 5 Conversion between DPD values and real values. ......................................................................... 77

vi

List of Figures
Figure 1. Schematic of an atom trajectory simulated with the MD method. ................................................ 4
Figure 2. Lennard Jones Potential. ................................................................................................................ 6
Figure 3. Schematic of PBC. ......................................................................................................................... 7
Figure 4. Examples of random BNP microstructures in different periodic cells… .................................... 18
Figure 5. A sphere in reciprocal space showing a Fibonacci grid with 2000 points… .............................. 20
Figure 6. Ligament size distribution for different periodic cells… ............................................................ 23
Figure 7. Numerical measurements of Mean and Gaussian curvature plots… ........................................... 26
Figure 8. Scaled genus density versus relative density for different periodic cells… ................................ 27
Figure 9. GNP Structure. ............................................................................................................................. 29
Figure 10. Double GNP. ............................................................................................................................. 29
Figure 11. Morphology of the BNPMG sample investigated… ................................................................. 34
Figure 12. Tensile loading global and local stresses evolution… ............................................................... 36
Figure 13. Snapshot of the sample at 0 tensile strain… .............................................................................. 37
Figure 14. Probability distribution of the components of a unit normal vector… ...................................... 38
Figure 15. Strain localization and failure of NMG… ................................................................................. 40
Figure 16. The local temperature in a deformed BNPMG at 𝜀 =0.12. ....................................................... 42
Figure 17. Stress-strain curves for the BNPMG sample and MG ligaments of different diameters… ...... 44
Figure 18. Snapshots of failure by necking process of the MG ligaments. ................................................ 44
Figure 19. Initial structures for relative density 𝜌 =0.3, 0.4, 0.5, 0.6, and 0.7 samples… ......................... 50
Figure 20. Ligament size distributions for the initial structures… ............................................................. 51
Figure 21. Binary tree representations of symbolic expressions… ............................................................ 52
Figure 22. Tensile loading engineering stress-strain relationships… ......................................................... 57
Figure 23. Young’s modulus and ultimate tensile strength for samples… ................................................. 58
Figure 24. Young’s modulus and ultimate tensile strength for Cu 64Zr 36 samples… .................................. 59
Figure 25. Deformation and failure of Cu 64Zr 36 BNPMG samples… ......................................................... 61
Figure 26. Comparison between existing and GP-derived scaling laws accuracy… .................................. 63
Figure 27. DPD simulation models. (a) Unit nanoporous media morphology used… ............................... 74
Figure 28. Illustrations of the S40 model during the initial relaxation process… ...................................... 78
Figure 29. Pressure evolution in the source and sink regions… ................................................................. 80
Figure 30. Local pressure distribution in the systems during the fluid flow process… ............................. 81
Figure 31. Position of the pistons during the simulation… ........................................................................ 82
Figure 32. Fluid velocity distribution in the different models… ................................................................ 83
Figure 33. Mean square displacement (MSD) of fluid particles… ............................................................. 85
Figure 34. Particle flux per unit area under different pressure differences for different samples. ............. 87
Figure 35. Permeability dependence on pore size compared with predictions… ....................................... 89
Figure 36. Fluid flow in cylindrical channels of different diameters… ...................................................... 91

vii

Abbreviations
Abbreviations Explanation
BNP Bicontinuous Nanoporous
GNP Gradient Nanoporous
FCC Face Centered Cubic
MG Metallic glass
BMG Bulk Metallic glass
PMG Porous Metallic Glass
NPMG Nanoporous Metallic Glass
BNPMG Bicontinuous Nanoporous Metallic Glass
MD Molecular Dynamics
PBC Periodic Boundary Conditions
DPD Dissipative Particle Dynamics
mDPD Many-Body Dissipative Particle Dynamics
MSD Mean Square Displacement
AI Artificial Intelligence
ML Machine Learning
GP Genetic Programming
GPSR Genetic Programming Assisted Symbolic Regression

viii

Abstract
Nanoporous materials have been drawing increasing attention in the last decades due to their
lightweight and large surface area to volume ratio. The materials have proved their usefulness in many
applications, such as catalysts, sensors, and lightweight structural design. This research is focused on the
mechanical and hydraulic properties of nanoporous materials, which are investigated using a combination
of multi-scale simulation methods. First, we propose a direct method to generate three-dimensional (3D)
stochastic bicontinuous nanoporous (BNP) microstructures on general periodic parallelepiped cells,
providing flexibility to the design of nanoporous simulation models in different length scales. The method
combines the ideas of Cahn with an adjustable Fibonacci grid to ensure the generation of a stochastic
structure, while strictly enforcing periodicity in an arbitrary parallelepiped cell. The topological and
morphological features, including ligament size, mean and Gaussian curvatures, and genus, are investigated
to validate the method. The proposed method offers an efficient framework for the generation of 3D
stochastic periodic bicontinuous microstructures. The method is particularly convenient in the design of
adjustable nanoporous models for atomistic simulations. With the support of the novel method, we
investigated the mechanical properties of BNP metallic glasses (BNPMGs), which were recently
experimentally synthesized.
We conduct molecular dynamics (MD) simulations of tensile loading deformation and failure of
Cu 64Zr 36 BNPMG with 55% porosity and 4.4 nm ligament size. Results indicate an anomalous mechanical
behavior featuring delocalized plastic deformation preceding ductile failure. The deformation follows two
mechanisms: i) Necking of ligaments aligned with the loading direction and ii) progressive alignment of
randomly oriented ligaments. Failure occurs at 0.16 strain, followed by a massive rupture of ligaments. This
work indicates that a BNP design is able to effectively delocalize strain localization in a metallic glass (MG)
due to a combination of size effects on the ductility of MGs, resulting in nano ligaments necking and
progressive asynchronous alignment of ligaments. We further investigate the effect of porosity levels and
compositions on the mechanical behavior of such materials. MD simulations are employed to study the
ix

mechanical properties of nanoporous CuxZr 1-x MGs with five different compositions, 𝑥 = 0.28, 0.36, 0.5,
0.64, and 0.72, and porosity in the range 0.1 < 𝜙 < 0.7. Results from tensile loading simulations indicate a
strong dependence of the Young’s modulus (E) and ultimate tensile strength (UTS) on porosity and
composition. By increasing the porosity from 𝜙 = 0.1 to 𝜙 = 0.7, the topology of the nanoporous MG shifts
from closed cell to open-cell bicontinuous. The change in nanoporous topology enables a brittle-to-ductile
transition in deformation and failure mechanisms from a single critical shear band to necking and rupture
of ligaments. Genetic Programming (GP) is employed to find scaling laws for E and UTS as a function of
porosity and composition. A comparison of the GP-derived scaling laws against existing relationships
shows that the GP method is able to uncover expressions that can predict accurately both the values of E
and UTS in the whole range of porosity and compositions considered.
Other than mechanical properties, a good understanding of the hydraulic behavior of nanoporous media
is necessary for its wide application field, including electrokinetic energy conversion devices, water
desalination using nanosized membranes, and extraction of hydrocarbons from ultra-tight shales. To
address this, we employed Many-body Dissipative Particle Dynamics simulations to shed light on the fluid
flow process through BNP media, which are representative models for a wide class of nanoporous materials.
BNP models are constructed considering a defined morphology and porosity level and varying pore sizes.
Fluid flow is generated through the nanoporous media by the action of two pistons added to the ends of the
system. Simulation results, obtained by imposing different pressure differences along with the nanoporous
model, indicate a linear pressure drop within the nanoporous structure, regardless of pore size. Regardless
of the complexities and different scales of the porous media considered, the steady-state fluid flow through
the nanoporous models is proportional to the pressure gradient applied, in agreement with Darcy’s law. The
calculated pore-size dependence of permeability is well described by the Hagen-Poiseuille law, derived for
single tubes, considering a single shape correction factor that accounts for the flow resistance due to the
complex nanoporous morphology.
1

Chapter 1
Introduction
1.1 Nanoporous Materials
In the past decades, nanoporous materials have been applied in diverse technological areas, including
sensors, catalysis, actuation, fuel cells, microfluidics, and biomedical devices.
1–6
The wide applicability is
driven by their outstanding set of properties, such as high thermal and electrical conductivity, high surface-
to-volume ratio, and remarkably ten times the theoretical strength predicted from macroporous foam scaling
laws.
7–10
The ease of manufacturing nanoporous stochastic structures further makes them attractive.
11,12
An
usual approach to fabricate nanoporous materials is to chemically dealloy a bulk binary compound.
12–16
The
structural properties, e.g., porosity and ligament size, can be tailored through parameters such as
temperature, dealloying time, and initial composition.
17–27

Many experimental and theoretical efforts have been devoted to exploring the properties of
nanoporous materials. The assessment of the mechanical properties of nanoporous materials in experiments
is often done by nanoindentation
7,28
instead of a direct evaluation of tensile yield strength by conducting a
tensile loading experiment, which is challenging due to the small size of experimental samples. As a result,
numerical studies, especially atomistic simulations, are conducted to study the mechanical properties and
failure mechanisms of nanoporous materials. To achieve this, a realistic model representing nanoporous
material is necessary. Phase-field simulations of the spinodal decomposition process,
10,29,30
kinetic Monte
Carlo
29,31,
and Molecular Dynamics (MD) simulations
32
have been used to produce bicontinuous
nanoporous (BNP) structures that are similar to experimental samples. However, the methods mentioned
all require a considerable amount of computational resources.
In this work, a method combining the ideas of Cahn
33
with an adjusted Fibonacci grid, to ensure
the generation of a stochastic structure while strictly enforcing periodicity in an arbitrary parallelepiped
2

cell, is proposed. Mechanical and hydraulic properties of nanoporous materials models generated with this
method are then studied with multi-scale numerical approaches.
1.2 Metallic Glasses
Metallic glasses (MGs) are an intriguing class of metallic materials that are typically obtained by
rapid quenching of a molten metal alloy at a very high cooling rate. The first MG ever produced is an AuSi
alloy obtained by Klement et al.
34
at Caltech using a 10
6
K/s quenching rate. At such a high quenching rate,
atoms do not have sufficient time to arrange themselves into crystalline structures and form an amorphous
structure instead. The absence of long-range order in such a structure eliminates common crystalline defects
such as dislocations and grain boundaries. Due to the absence of crystal plasticity and crystalline defects,
MGs often possess higher strength, hardness, elasticity, and elastic energy storage capacity compared to
their crystalline counterparts. For some compositions, MGs may also have high corrosion resistance.
35–38

As a result, MGs are considered ideal candidates for structural applications. However, the negligible
plasticity displayed by most MG alloys prevents their widespread applications. Efforts have been placed on
improving the ductility of MGs. Introducing pores is one of the effective, proven ways to enhance their
plasticity.
39
In particular, bicontinuous porous MGs synergize the outstanding properties of MGs and open-
cell porous materials and have the potential to enhance the applicability of MGs in catalysis, sensors, and
lightweight structural designs. Advanced functional use of porous MGs may even demand a larger specific
surface area than currently available. This demands a reduction in the size of the pores into nanoscales.
Recently, Jiao et al.
40
successfully fabricated a 3D BNP MG (BNPMG) by selective corrosion of spinodally
decomposed MG precursors and passivation. The stochastic nature of the spinodal process and the nano-
length scale of the porous features resulted in a high specific surface area, making such material an excellent
candidate for gas absorption.
40
In this work, an investigation of the mechanical properties of BNPMGs is
conducted by employing MD simulations. Results indicate that a BNP design is able to effectively
delocalize strain localization in an MG due to a combination of size effects on the ductility of MGs, resulting
in nano ligaments necking and progressive asynchronous alignment of ligaments.
3

1.3 Fluid Flow in Nanoporous Media
Human society demands an ever-increasing supply of energy to support the rapid development of
the world economies. Studies estimate a 50% increase in energy demand by 2050.
41,42
To satisfy this large
amount of energy demand and concurrently reduce or limit the carbon footprint, new technologies to
efficiently and sustainably produce energy are required. To list a few major applications that need
innovative technologies, electrokinetic energy conversion devices,
43
water desalination using nanosized
membranes,
44
storage of hydrogen,
45
and extraction of hydrocarbons from ultra-tight shales
46
are attracting
more and more attention. All these technologies can be improved by enhancing our understanding of fluid
flow through nanoporous media. Indeed, many experimental and theoretical efforts have been placed on
the investigation of fluid flow in nanoporous media, especially those with extremely low permeability.
47–49

Scientific modeling and simulations offer flexible and versatile tools to understand, define, quantify,
and visualize the physical behavior of a fluid flow. Continuum scale computational fluid dynamics and
nanoscale MD are widely used for modeling fluid flow in porous media. However, continuum scale
simulation might not be able to capture the physical and chemical behaviors of fluid particles in nanoporous
media,
50,51
and atomistic simulations are often limited by their small length and time scales. Mesoscale
modeling approaches, such as dissipative particle dynamics (DPD), can bridge the gap between continuum
and nanoscale methods. In this work, the DPD method is employed to study fluid properties in BNP media
and the effect of the geometrical features of the BNP media on fluid flow.
4

Chapter 2
Methodology
2.1 Molecular Dynamics
2.1.1 Overview of the Molecular Dynamics Method
MD is a versatile simulation method to study physical and chemical processes of materials at atomic
and nano length scales. In MD simulations, the trajectory of each atom is calculated at each time step based
on Newton’s equations of motion. The movement of atoms can reveal information on various properties of
materials, such as modulus of elasticity, strength, and ductility.

Figure 1. Schematic of an atom trajectory simulated with the MD method.
In an MD simulation, a system with 𝑁 particles is described by the corresponding positions and
velocities, i.e., '(𝑟
!
(𝑡)
----------⃗
,𝑣
!
(𝑡)
-----------⃗
12𝑖 =0,1,2,3,…,𝑁−1:. The initial positions and velocities at 𝑡 = 0 must
be provided to initialize the simulation. Then, at each time step, the positions and velocities are updated
based on the force applied to each atom 𝐹
"
--⃗
(𝑡). A typical MD process is illustrated in Table 1.
To calculate the force on each atom, an interatomic potential describing the chemical and physical
interaction between atoms in the system is required. This is a key ingredient in MD simulations and it must
be accurate enough for the MD simulation to be able to reproduce realistic physical and chemical
phenomena. The force is calculated as the derivative of the potential with respect to spatial coordinates.
Step three is the time integration step, and common algorithms for step 3 are Leap-Frog and Velocity-Verlet.
Introductions of interatomic potentials and the Velocity-Verlet algorithm are given in the following sections.
5

Input: '(𝑟
!
(0)
-----------⃗
,𝑣
!
(0)
------------⃗
12𝑖 =0,1,2,3,…,𝑁−1:
Step 𝑡 =1,𝑠𝑡𝑒𝑝_𝑚𝑎𝑥:
1. Compute 𝐹
"
--⃗
(𝑡) for all atoms
2. Compute 𝑎
"
---⃗(𝑡) for all atoms based on 𝑚𝑎 ⃗ =𝐹
⃗

3. Update position and velocity, (𝑟
!
(𝑡)
----------⃗
,𝑣
!
(𝑡)
-----------⃗
1for each atom, using 𝑎
"
---⃗(𝑡) and
(𝑟
!
(𝑡−1)
-------------------⃗
,𝑣
!
(𝑡−1)
--------------------⃗
1
Table 1 Pseudocode of an MD simulation.
2.1.2 Interatomic Potential
In MD simulation, the potential energy, 𝑉(𝑟
!
----⃗), is the starting point and it depends only on the
atomic positions. The force experienced by the 𝑘
th
atom at a given positional configuration can be calculated
as follows:
𝐹
#
----⃗
=−
𝜕
𝜕𝑟
#
---⃗
𝑉(𝑟
!
)
------⃗
=−E
𝜕𝑉(𝑟
!
)
------⃗
𝜕𝑥
#
,
𝜕𝑉(𝑟
!
)
------⃗
𝜕𝑦
#
,
𝜕𝑉(𝑟
!
)
------⃗
𝜕𝑧
#
H 1
where the potential energy 𝑉(𝑟
!
----⃗) can be expressed, in one of the most simple and fundamental ways, as
𝑉(𝑟
!
)
------⃗
=I𝑢(𝑟
$%
1
$&%

2
with 𝑟
$%
= 2𝑟
"'
---⃗2 the distance of atom 𝑖 and atom 𝑗. In this case, 𝑢(𝑟
$%
1 is a simple pairwise interatomic
potential. To give an example, a widely used pairwise interatomic potential is the Lennard Jones potential,
which can describe well materials such as liquid argon and inert gases. It takes the form
𝑢(𝑟
$%
1=4𝜖EN
𝜎
𝑟
$%
P
()
−N
𝜎
𝑟
$%
P
*
H 3
It has two tunable parameters, 𝜎 and 𝜖, which can be obtained by fitting the cohesive energy and lattice
constant of different materials.
It can be seen from the curve of Lennard Jones that it consists of a short-range repulsion and long-
range attraction. The short-range repulsion becomes infinitely large when the distance is close to zero. This
6

is consistent with the Pauli exclusion principle, such that no two particles can be in the same state, i.e.,
overlap their electron cloud. The attractive long-range term is used to describe van der Waals interactions
between atoms.
However, interatomic potentials can take much more complex forms and describe interactions
beyond simple pairwise functionals. For example, some interatomic potentials consider explicitly the angle
made by chemical bonds 𝜃
$%#
.
52
Potentials such as the embedded atom model (EAM) also takes into account
the electron density to correctly describe the behavior of metals and metal alloys.
53

Figure 2. Lennard Jones Potential.
2.1.3 Velocity-Verlet Algorithm
In order to describe the trajectories of atoms in an MD simulation, the trajectories are discretized,
and become a sequence of points in position-velocity space: '(𝑟
!
(0)
-----------⃗
,𝑣
!
(0)
------------⃗
1, (𝑟
!
(Δ)
-----------⃗
,𝑣
!
(Δ)
------------⃗
1,
(𝑟
!
(2Δ)
--------------⃗
,𝑣
!
(2Δ)
---------------⃗
1,…: where Δ stands for the time step. The popular and robust Velocity-Verlet algorithm
can be used to integrate the equations of motion. The algorithm can be derived using a Taylor expansion
and it is summarized in Table 2.
The advantage of the Velocity-Verlet algorithm is that it rigorously conserves the position-velocity
volume, regardless of the time step. This significantly enhances the stability of the simulation. However, it
is still important to choose an appropriate time step since a time step too large will introduce inaccuracy,
7

while a time step too small will consume more time and computational resources. The time step of MD
simulations should be well below the period of the highest vibrational frequency of the atoms in the system.
Input: (𝑟
"
--⃗(0),𝑣
"
---⃗(0)) for all 𝑖
Compute 𝑎
"
---⃗ based on {𝑟
"
--⃗(0)}
Step 𝑡 =1,𝑠𝑡𝑒𝑝_𝑚𝑎𝑥:
1. 𝑣
"
---⃗ ← 𝑣
"
---⃗+𝑎
"
---⃗Δ/2 for all 𝑖
2. 𝑟
"
--⃗← 𝑟
"
--⃗+𝑣
"
---⃗Δ for all 𝑖
3. Compute 𝑎
"
---⃗ based on {𝑟
"
--⃗}
4. 𝑣
"
---⃗ ← 𝑣
"
---⃗+𝑎
"
---⃗Δ/2 for all 𝑖
Table 2 Pseudocode of the Velocity-Verlet algorithm.
2.1.4 Periodic Boundary Conditions
The number of atoms is always finite in MD simulations. To study bulk materials properties, a
system with no surfaces is required. To achieve this, periodic boundary conditions (PBC) are commonly
applied to the simulation box.

Figure 3. Schematic of PBC.
Under PBC, an atom is always brought back to the other side of the simulation box, whenever it
crosses the boundary of the box. While we simulate only the atoms in the simulation box, the outcome
properties are representative of a bulk material. In practice, some properties may be affected by the size of
8

the simulation cell, such as phonon properties. Nevertheless, the use of PBC prevents any effects coming
from the presence of free surfaces. The following equation is how a PBC in MD is implemented:
⎩
⎪
⎨
⎪
⎧𝑥
$
=𝑥
$
−SignRb
𝐿
+
2
,𝑥
$
d−SignR(
𝐿
+
2
,𝑥
$
−𝐿
+
)
𝑦
$
=𝑦
$
−𝑆𝑖𝑔𝑛𝑅i
𝑦
2
,𝑦
$
j−SignR(
𝐿
,
2
,𝑦
$
−𝐿
,
)
𝑧
$
=𝑧
$
−SignRb
𝐿
-
2
,𝑧
$
d−SignR(
𝐿
-
2
,𝑧
$
−𝐿
-
)
4
where
SignR(𝑎,𝑏)=l
𝑎 𝑖𝑓 𝑏 >0
−𝑎 𝑖𝑓 𝑏 ≤0
5
As a result of PBC, one atom will not only interact with all the atoms in the cell but also interacts
with images of all atoms in other simulation boxes, which are duplicates of the original box and stack along
all three directions. However, the interaction between atoms that are far away from each other is ignored
since the magnitude of such interaction is negligible. In practice, for a simulation box (𝐿
+
×𝐿
,
×𝐿
-
), we
can ignore the contribution from atoms that are further than 𝑟 >
./0 12
!
,2
"
,2
#
4
)
. This means, for each pair of
atoms 𝑖 and 𝑗 there is at most one interaction to be considered among all the possible periodic images.
Whether it is the interaction between these two atoms in the original cell, or it is the interaction between
the original 𝑖 and an image of 𝑗 depends on the difference in their position in the original cell, 𝑟 ⃗
$%
=
(𝑥
$%
,𝑦
$%
,𝑧
$%
). The final positional difference between 𝑖 and 𝑗 can be expressed as follows:
⎩
⎪
⎨
⎪
⎧𝑥
$%
=𝑥
$%
−SignRb
𝐿
+
2
,𝑥
$%
+
𝐿
+
2
d−SignRb
𝐿
+
2
,𝑥
$%
−
𝐿
+
2
d
𝑦
$%
=𝑦
$%
−SignRb
𝐿
,
2
,𝑦
$%
+
𝐿
,
2
d−SignRb
𝐿
,
2
,𝑦
$%
−
𝐿
,
2
d
𝑧
$%
=𝑧
$%
−SignRb
𝐿
-
2
,𝑧
$%
+
𝐿
-
2
d−SignRb
𝐿
-
2
,𝑧
$%
−
𝐿
-
2
d
6
This is called the minimum image convention.
9

2.1.5 Shifted Potential
The interatomic potential always decays and can be neglected when the distance between two atoms
is large. It is physically intuitive since two particles that are far away from each other should not affect each
other. In practice, during MD simulations, there is always an interatomic cutoff distance (𝑟
5
) to be set. This
is required to define when to neglect the interaction of atoms that are further from each other. In such a way,
computation time and resources can be saved. Commonly, interatomic potentials can be expressed as
follows:
𝑢(𝑟) = l
𝑢(𝑟) 𝑟 ≤ 𝑟
5
0 𝑟 >𝑟
5
7
However, this will introduce a problem. Suppose that 𝑢(𝑟 =𝑟
5
) ≠0. Then the potential function is not
continuous and thus not differentiable with respect to 𝑟. Then the force at 𝑟
5
is not well defined. To avoid
this problem, one must modify the potential to 𝑢’(𝑟) such that 𝑢’(𝑟
5
) = 0,
67’(:;:
$
)
6:
= 0 and 𝑢’(𝑟) is
differentiable over for all 𝑟. The modified potential is the so-called shifted potential, and it takes the
following form:
𝑢’(𝑟) =v
𝑢(𝑟)−𝑢(𝑟
5
)−(𝑟−𝑟
5
)
d𝑢
d𝑟
x
:;:
$
𝑟 ≤𝑟
5
0 𝑟 >𝑟
5
8

2.2 Dissipative Particle Dynamics
2.2.1 Formulation of the Dissipative Particle Dynamics Method
DPD is a coarse-grained particle dynamics model. Instead of single atoms, the particles represent
clusters of molecules that interact via conservative, dissipative, and fluctuating forces. The DPD method
can overcome the small spatial and time scale of MD. The formulation of the DPD method can be shown
in the following expression:
10

𝑭
$%
=𝑭
$%
=
+𝑭
$%
>
+𝑭
$%
?
9
where 𝑭
$%
=
represents a random stochastic term, 𝑭
$%
>
a dissipative term, and 𝑭
$%
?
a conservative term.
Considering 𝒓
$
, the position of the i
th
particle, and 𝒗
$
, the velocity of the i
th
particle, the conservative force
is described as follows:
𝑭
$%
?
=𝐴𝑤(𝑟
$%
)𝒓
$%

10
where 𝐴 is the magnitude of the force, 𝑟
$%
denotes the distance between the particle 𝑖 and particle 𝑗 and 𝒓
$%

is the unit vector pointing from particle 𝑖 to particle 𝑗. The weighting factor, 𝑤(𝑟
$%
1, varying from 0 to 1 is
expressed as follows:
𝑤(𝑟
$%
1=1−𝑟
$%
/𝑟
5
11
where 𝑟
5
is the cutoff for the DPD potential beyond which the conservative force is zero. The dissipative
part is described as:
𝑭
$%
>
=−𝛾𝑤
>
(𝑟
$%
1(𝒓
$%
∙𝒗
$%
)𝒓
$%
12
where 𝛾 is the magnitude of the dissipative force. Finally, the random term is given by:
𝑭
$%
=
=𝜎𝑤
=
(𝑟
$%
1𝛼(∆𝑡)
(/)
𝒓
$%
13
where ∆𝑡 is the simulation time step, 𝛼 is Gaussian white noise that is equal for 𝑭
$%
=
and 𝑭
%$
=
to ensure
momentum conservation. It satisfies also the following stochastic properties:

〈𝛼
$%
(𝑡)〉=0
〈𝛼
$%
(𝑡)𝛼
$A%A
(𝑡)〉= (𝛿
$$
%𝛿
%%
% +𝛿
$%
%𝛿
%$
%1𝛿(𝑡−𝑡′)
14
Based on Newton’s second law, we can derive a Langevin equation from Eq. (9) – Eq. (13):
d𝒑
$
= I𝑭
$%
5
(𝒓
$%
1
$B%
+I−𝛾𝑤
>
(𝒓
$%
1( 𝒓
$%
∙𝒗
$%
1 𝒓
$%
$B%
d𝑡+I𝜎𝑤
=
(𝑟
$%
1𝒓
$%
d𝑊
$%
$B%
15
where d𝑊
$%
=𝛼
$%
d𝑡. It satisfies the following relationship:
d𝑊
$%
d𝑊
$A%A
=(𝛿
$$
%𝛿
%%
% +𝛿
$%
%𝛿
%$
%)d𝑡 16
11

Suppose we have particle distribution function 𝜌(𝒓,𝒑,𝑡). We can derive the expression of 𝜌(𝒓,𝒑,𝑡+Δ𝑡).
Knowing the stochastic features of d𝑊
$%
and ignoring the term that is infinitesimal of more than order 2,
we can obtain the Fokker-Plank
54,55
equation:
𝜕𝜌(𝒓,𝒑,𝑡)
𝜕𝑡
=𝐿
?
𝜌(𝒓,𝒑,𝑡)+𝐿
>
𝜌(𝒓,𝒑,𝑡)
17
where the operators are defined as follows:
⎩
⎪
⎨
⎪
⎧
𝐿
?
≡−I
𝒑
$
𝑚
𝜕
𝜕𝒓
$
$
+ I 𝑭
$%
?
𝜕
𝜕𝒑
$
$,%B$

𝐿
?
≡ I 𝒓
$%
𝜕
𝜕𝒑
$
𝛾𝑤
>
(𝒓
$%
1( 𝒓
$%
∙𝒗
$%
1+
𝜎
)
2
𝑤
=
(𝑟
$%
1 𝒓
$%
N
𝜕
𝜕𝒑
$
−
𝜕
𝜕𝒑
%
P
$,%B$
18
When a steady-state is reached, the derivative of 𝜌 should be zero, and the distribution should be the
equilibrium distribution 𝜌
CD
. If the state of the system is described by the Gibbs canonical ensemble,
𝜌
CD
=
1
𝑍
exp−EI
𝑝
$
)
2𝑚
$
+𝑉(𝑟)
$
H/𝑘
E
𝑇 19
where 𝑉 is the potential function that gives rise to the conservative force 𝑭
?
, 𝑘
E
the Boltzmann’s constant,
𝑇 the equilibrium temperature, and 𝑍 the normalizing partition function. The Liouville operator 𝐿
?
on the
canonical distribution will result in 𝐿
?
𝜌
CD
(𝒓,𝒑,𝑡) = 0. We can further suppose the following relation to
make 𝐿
>
𝜌
CD
(𝒓,𝒑,𝑡)= 0, which eventually leads to a steady-state:
v
𝑤
=
(𝑟) =𝑤
>
(𝑟)
(
)
𝜎 = (2𝑘
E
𝑇𝛾)
(
)

20
This is known as the fluctuation-dissipation theorem for the DPD method.
2.2.2 Formulation of the Many-Body Dissipative Particle Dynamics Method
The original DPD model interactions are not sufficient to model multi-phase systems, such as
systems containing liquid-liquid interfaces, liquid-vapor interfaces, and liquid-vacuum interfaces. To
address this, an attractive term can be added to the conservative force. In the formulation of Warren,
34
the
12

repulsive term is also enhanced by density-dependent contributions. The modified conservative term is
given by:
𝑭
$%
?
=𝐴
$%
𝑤(𝑟
$%
1𝒓
$%
+𝐵
$%
(𝜌̅$
+𝜌̅%
1𝑤
F
(𝑟
$%
1𝒓
$%
21
where 𝐴
$%
<0 is the magnitude of the attractive term and 𝐵
$%
>0 is the magnitude of the repulsive term.
Both terms are dependent on the types of particle 𝑖 and particle 𝑗. 𝑤
F
(𝑟
$%
) is a weighting term with a shorter
cutoff 𝑟
F
<𝑟
5

𝑤
F
(𝑟
$%
1=1−𝑟
$%
/𝑟
F
22
In our simulations, 𝑟
F
is chosen to be 0.75 𝑟
5
. The local density at the position of particle i is computed as
follows:
𝜌̅$
= I𝑤
G
(𝑟
$%
1
$B%

23
with
𝑤
G
(𝑟) =v
15
2𝜋𝑟
F
H
b1−
𝑟
𝑟
F
d 𝑟 <𝑟
F
0 𝑟 >𝑟
F
24
The 𝑤
G
(𝑟
$%
) is defined in such a way that satisfies the normalization condition:
4𝜋𝑟
)
𝑤
G
(𝑟)
I
K
d𝑟 =1 25
To simulate a multiphase environment, one can adjust the attractive and repulsive terms for different
interactions between different types of particles. This model is easy to parameterize and can correctly
represent the surface tension of a fluid.
2.2.3 Time Integration
The time integration algorithm is crucial for DPD simulations. A poor time integration strategy will
result in inaccuracy. Early implementations of the DPD method made use of the Euler scheme:
13

v
𝒓
$
(𝑡+Δ𝑡) =𝒓
$
(𝑡)+Δ𝑡𝒗
$
(𝑡)
𝒗
$
(𝑡+Δ𝑡) =𝒗
$
(𝑡)+Δ𝑡𝒇
$
(𝑡)
𝒇
$
(𝑡+Δ𝑡) =𝒇
$
(𝒓
$
(𝑡+Δ𝑡),𝒗
$
(𝑡+Δ𝑡))
26
where 𝒓
$
is the position of atom 𝑖, 𝒗
$
is the velocity of atom 𝑖 and 𝒇
$
is the force experienced by atom 𝑖.
All the above quantities are dimensionless. The problem with the Euler scheme is the time irreversibility,
and this can lead to energy drift during the DPD simulations. To overcome this drawback, a modified
version of the Velocity-Verlet algorithm is proposed:
56

⎩
⎪
⎨
⎪
⎧ 𝒓
$
(𝑡+Δ𝑡)=𝒓
$
(𝑡)+Δ𝑡𝒗
$
(𝑡)+
1
2
(Δ𝑡)
)
𝒇
$
(𝑡)
𝒗
$
(𝑡+Δ𝑡) =𝒗
$
(𝑡)+𝜆Δ𝑡𝒇
$
(𝑡)
𝒇
$
(𝑡+Δ𝑡) =𝒇
$
(𝒓
$
(𝑡+Δ𝑡),𝒗
$
(𝑡+Δ𝑡))
𝒗
$
(𝑡+Δ𝑡)=𝒗
$
(𝑡)+
1
2
Δ𝑡(𝒇
$
(𝑡)+𝒇
$
(𝑡+Δ𝑡))
27
where 𝒗
$
(𝑡+Δ𝑡) is the prediction of the velocity at time 𝑡+Δ𝑡 and 𝜆 is an empirical parameter, which
represents the effects of stochastic interactions. Following this integration method, the predicted velocity is
used to obtain the force. Then the velocity is corrected in the last step. For a system where force depends
only on positions, the standard Velocity-Verlet algorithm can be recovered at 𝜆 =1/2. Optimal 𝜆 resides
in the range from 0 to 1 for different DPD systems.
14

Chapter 3
Generation of Non-Cubic Stochastic Periodic Bicontinuous
Nanoporous Structures
3.1 Background
Open-cell nanoporous materials have drawn large attention because of their outstanding set of
properties, including high thermal and electrical conductivity, high surface-to-volume ratio, and remarkably
ten times the theoretical strength predicted from macroporous foam scaling laws.
7–9,57–61
This set of
properties offers the possibility of promising applications in actuators,
62–64
catalysis,
65–67
and surface
chemistry-based sensing.
19,27,68,69
In addition, further applications include those where macroporous metals
have been conventionally employed e.g., heat sinks,
70,71
light structural sandwich panels,
72
and energy
absorption devices.
73,74
A common way to fabricate nanoporous structures is through the chemical
dealloying of a bulk binary compound.
12–16
The structural properties, e.g., porosity and ligament size, can
be controlled through parameters such as temperature, dealloying time, and initial composition.
17–27

The evaluation of the mechanical properties of nanoporous metals in experiments is often done by
nanoindentation.
28,75
Atomistic investigations, e.g., using MD simulations, have been increasingly used to
complement experiments and unveil the underlying deformation mechanisms.
11,76–81
Commonly, the
nanoporous topology required to perform atomistic simulations is generated in phase-field simulations of
the spinodal decomposition process.
10,30,37
Alternatively, kinetic Monte Carlo
29,31
and even direct MD
simulations
32
have been used for that purpose. While those methods are able to produce bicontinuous
structures similar to those produced experimentally, they require considerable effort and computational
resources, especially for large systems. This also creates unnecessary constraints on atomistic investigations
of different aspects related to nanoporous metals in particular those requiring the study of many different
configurations, the effect of ligament size, porous size, relative density, etc.
15

It is desirable to have an available method to generate BNP structures directly based on a few
chosen features such as relative density and ligament size. In its seminal work, Cahn
33
suggested that the
approximate solution of the now-known as Cahn-Hillard equation could be written as a set of superimposing
sinusoidal waves:
𝑓(𝒓) =
¡
2
𝑁
Icos(𝒒
𝒊
∙𝒓+𝜑
$
)
!
$;(

28
where 𝒓 is the position vector, §
)
!
the normalization factor, 𝑁 the number of waves, 𝒒
𝒊
spherically
symmetric random wave vectors, and 𝜑
$
a set of random phases evenly distributed in the range (0,2𝜋).
The wavenumber 𝒒
𝒊
should be chosen such that |𝒒
𝒊
| =𝑞
K
, where 𝑞
K
is a constant that controls the ligament
size. As a result of that, 𝑓(𝒓) is a Gaussian random field.
82,83
This approximate solution could be used to
generate the bicontinuous topology of a nanoporous structure. Using a given level cut 𝜉, the points where
𝑓(𝒓) <𝜉 can be defined as belonging to the solid phase. 𝑓(𝒓) =𝜉 indicates the interface (surface of the
solid phase) and 𝑓(𝒓)>𝜉 indicates the porous region. Soyarslan et al.
84
modified this model to generate
periodic bicontinuous structures in a cubic cell by selecting a finite number of waves of which 𝒒
𝒊
are in the
form:
𝒒
𝒊
=
2𝜋
𝑎
(ℎ
$
,𝑘
$
,𝑙
$
) 29
where 𝑎 is the cubic cell edge length and ℎ
$
, 𝑘
$
, 𝑙
$
are integers, with the requirement that
𝐻 =
§
ℎ
$
)
+𝑘
$
)
+𝑙
$
)

30
is constant. The method proposed by Soyarslan et al.
84
offers a convenient and direct way to generate models
for atomistic MD simulations, which commonly are based on periodic systems. PBC is strictly enforced
by the use of Eq. (29) and Eq. (30), while Eq. (28), and the level cut 𝜉 are used to define the solid region of
the nanoporous structure. The method allows for the efficient generation of larger periodic BNP systems
16

and the convenience of generating as many independent configurations as desired, since it involves a minor
computational effort.
While the method of Soyarslan et al.
84
offers an efficient method for the generation of bicontinuous
microstructures in cubic cells, very often, such microstructures are desired or required in non-cubic cells.
For example, high aspect ratio samples are commonly used in mechanical tests.
11,75,76
Non-cubic cells are
also desirable in studies of membranes,
85
composites,
86
and functionally graded samples.
87,88
A non-cubic
cell is also required in certain cases where one chooses to build a periodic bicontinuous structure model
composed of a single crystal with non-cubic symmetry, such as Bi, As, Sb, Te, etc. Therefore, a method for
the generation of random bicontinuous structures in general non-cubic cells is highly desirable to provide
flexibility to studies of different materials, composites, and designs.
In this chapter, based on the original ideas of Cahn
33
and Soyarslan et al.,
84
we propose a direct and
efficient computational method to generate stochastic BNP structures with periodicity in an arbitrary
parallelepiped cell, i.e., with no constraints on dimensions and shape. We first discuss the proposed method,
then we validate the method by comparing topologies of the structures in cubic, rectangular prisms, and a
general parallelepiped to that of other methods and experiments. Finally, we present our conclusions.
3.2 Generation of Nanoporous Structures in Periodic Cells
3.2.1 Cubic Periodic Cells
Before presenting our method, it is instructive to further discuss the method proposed by Soyarslan
et al.
84
for cubic periodic cells. As a reference, we use their method and reproduce a sample representative
BNP structure in a cubic cell with PBC. For that, we choose to fill a cubic cell of edge 𝑎 = 100 nm with a
gold Face Centered Cubic (FCC) crystal structure, with lattice parameter 4 Å. The value of the 𝐻 constant,
defined by Eq. (30), is set to be √161, which results in a constant wave number:
|𝒒|=𝑞
K
=
2𝜋𝐻
𝑎
=
𝜋√161
50
nm
M(

31
17

The value 𝐻 = √161 is chosen following the work of Soyarslan et al.
84
to maximize the number of
possible wave vectors while keeping the expected number of elements about 10 per edge length. With that
𝐻 choice, there are four possible sets of wave vectors satisfying Eq. (31): 〈12,4,1〉, 〈11,6,2〉, 〈10,6,5〉, and
〈9,8,4〉. Each set represents 24 independent vectors, e.g., 〈12,4,1〉 includes vectors (12,4,1), (12,1,4),
(4,12,1), (−12,4,1), (12,−4,1), etc. The vectors (−12,−4,-1), (12,−4,−1), etc., are not included since they
can be generated from the vectors of the set by a constant multiplication, i.e., −1. The combination of the
possible sets results in 96 independent wave vectors. For each wave, we assign a random phase 𝜑
$
, as
defined in Eq. (28). With the wave vectors and random phases assigned, the function 𝑓(𝒓), as defined in
Eq. (28), is calculated at the position of each atom and its value is compared to the level cut 𝜉, which
controls the relative density. To obtain a nanoporous structure with 50% relative density, we use 𝜉 =0.
The resulting cubic BNP structure is demonstrated in Figure 4(a). The topology of the nanostructure is
highlighted with the use of a smooth triangular surface mesh as implemented in OVITO.
89

3.2.2 Arbitrary Parallelepiped Periodic Cells
With a reference cubic system generated, we describe the general method that we propose to
generate BNP structures in arbitrary periodic cells, including cubic cells as just described in section 3.2.1.
For a general parallelepiped system with edge vectors 𝒂, 𝒃 and 𝒄, 𝑎 ≠𝑏 ≠𝑐,𝛼 ≠𝛽 ≠𝛾 where 𝑎,𝑏,𝑐 are
the length of 𝒂, 𝒃, 𝒄, and 𝛼,𝛽,𝛾 are the angles between (𝒃,𝒄), (𝒂,𝒄) and (𝒂,𝒃) respectively, the PBC is
satisfied if
𝑓(𝒓) =𝑓(𝒓+𝑚𝒂+𝑛𝒃+𝑝𝒄) 32
where 𝑚, 𝑛, and 𝑝 are integers. The above condition is fulfilled if we have:
v
cos(𝒒∙𝒓)=cos(𝒒∙𝒓+𝒒∙𝒂)
cos(𝒒∙𝒓) =cos(𝒒∙𝒓+𝒒∙𝒃)
cos(𝒒∙𝒓) =cos(𝒒∙𝒓+𝒒∙𝒄)
33
which implies
18

µ
𝒒∙𝒂=𝑖
(
×2𝜋
𝒒∙𝒃=𝑖
)
×2𝜋
𝒒∙𝒄=𝑖
H
×2𝜋
34
where 𝑖
(
, 𝑖
)
, and 𝑖
H
are integers.

Figure 4. Examples of random BNP microstructures in different periodic cells. (a) cubic, (b) rectangular
prism, and (c) general parallelepiped cell.
The cartesian coordinates of 𝒂, 𝒃, and 𝒄 can be written as follows:
(𝒂,𝒃,𝒄) =𝑴
M(
𝑰 35
where 𝑰 is the identity matrix and 𝑴
M(
is the orthogonalization matrix:
19

𝑴
M(
=
⎝
⎜
⎜
⎛
𝑎 𝑏cos𝛾 𝑐cos𝛽
0 𝑏sin𝛾
𝑐(cos𝛼−cos𝛽cos𝛾)
sin𝛾
0 0
𝑐(1−cos
)
𝛼−cos
)
𝛽−cos
)
𝛽+2cos𝛼cos𝛽cos𝛾)
(/)
sin𝛾 ⎠
⎟
⎟
⎞

= E
𝑀
((
𝑀
()
𝑀
(H
0 𝑀
))
𝑀
)H
0 0 𝑀
HH
H
36
with that, q=(𝑞
(
,𝑞
)
,𝑞
H
) must satisfy
𝒒∙𝑴
M(
2𝜋
=(𝑖
(
,𝑖
)
,𝑖
H
)
37
Similar to Eq. (31), Eq. (37) implies that a small value of the |𝒒| constant will lead to a small number of
(𝑖
(
,𝑖
)
,𝑖
H
). In order to generate a BNP structure composed of randomly oriented ligaments, we must
consider a large set of (𝑖
(
,𝑖
)
,𝑖
H
), therefore a large |𝒒| value. However, larger |𝒒| values result in a larger
number of elements per edge length, corresponding to smaller ligament sizes for a given simulation box.
That may create a practical problem since setting a large |𝒒| to ensure randomness for a relevant ligament
size may require extremely large cells. In the proposed method, in order to address that, we consider a range
of possible values around a desired |𝒒| value, instead of setting a single value for |𝒒|.
In the proposed method, to maximize the randomness associated with the stochastic BNP structure,
one should consider a large number of independent waves. For the validation of the method, we start by
considering 2000 waves 𝑖 with a desired |𝒒
$
A
| value
𝒒
$
A
=(𝑞
$(
A
,𝑞
$)
A
,𝑞
$H
A
) 38
uniformly distributed over the solid angle 4𝜋 using a spherical Fibonacci grid.
90,91

One should note that there are many approaches available to address the problem of distributing
evenly a large set of points over the surface of a sphere.
92
However, the spiral set of points as provided by
a Fibonacci grid offers an elegant and simple way to solve the problem of evenly distributing grid points
without imposing unnecessary artificial symmetries. The Fibonacci grid is constructed using the angle
based on the golden ratio, 𝜓 =(1+√51 2 ⁄ , as described below. Fibonacci grids offer a computationally
20

efficient way to distribute the load in parallel computing of climate models.
91,93
They are also useful in
diverse fields such as quasicrystals,
94,95
photonic waveguides,
94
and evacuation strategies.
96
Fibonacci grids
can be conceived as mathematically generalizations of patterns spontaneously occurring in nature, e.g., in
leaves, seeds, and spirals of different plants. Well-known examples are the arrangements of seeds in
sunflower heads and pineapples. These arrangements are uniform and highly isotropic, features useful in
the generation of the independent wave vectors required in this work.

Figure 5. A sphere in reciprocal space showing a Fibonacci grid with 2000 points, indicating the location
of the wave vectors before and after modification. (a) Original Fibonacci grid over a spherical surface of
radius |𝒒
𝒊
| =𝒒
𝟎
. Lines connect the nearest points. (b) The modified Fibonacci grid is shown in red with
the original grid in black for reference.
To generate the Fibonacci grid, we start by using polar coordinates. From the center (0,0), the 2D
Fibonacci grid is constructed as a series of points where the polar coordinate of the 𝑖th point is (𝑟
$
,𝜃
$
) =
(𝑐
:
√𝑖,𝑖𝜔1, where 𝑐
:
is a constant and 𝜔 the golden angle, ~137.508°, defined as 2𝜋(1−𝜓
M(
), where 𝜓 is
21

the golden ratio aforementioned. Slight modifications are applied to obtain a Fibonacci grid on the surface
of a 3D sphere. We start by creating a grid on a unit sphere. The cylindrical coordinate of the starting point
is (0,0,−1) and the coordinate of the 𝑖th point is (𝑟
$
,𝜃
$
,𝑧
$
) = N§1−𝑧
$
)
,𝑖𝜔,−1+2(𝑖−1)/(𝑛−1)P,
where 𝑛 = 2000 is the total number of grid points we use. After converting the cylindrical coordinates to
cartesian coordinates, we multiply the resulting vectors by |𝒒
$
A
| to find the actual wave vectors. The quasi-
equidistant grid obtained in such a way ensures the ideal isotropy of the waves required to generate the
random bicontinuous structure.
To enforce periodicity in the final structure we take each generated wave, find the closest integer
𝑖
(
to 𝑞
$(
A
𝑀
((
/2𝜋, and obtain 𝑞
$(
as follows:
𝑞
$(
=𝑖
(
×2𝜋/𝑀
((
39
Subsequently, we find the closest integer 𝑖
)
to (𝑞
$(
𝑀
)(
+𝑞
$)
A
𝑀
))
)/2𝜋, then obtain 𝑞
$)
as follows:
𝑞
$)
=(𝑖
)
×2𝜋−𝑞
$(
𝑀
)(
)/𝑀
))
40
Similarly, we find the closest integer 𝑖
H
to (𝑀
H(
+𝑞
$)
𝑀
H)
+𝑞
$H
A
𝑀
HH
)/2𝜋 and obtain 𝑞
$H
as follows:
𝑞
$H
=(𝑖
H
×2𝜋−𝑞
$(
𝑀
H(
−𝑞
$)
𝑀
H)
)/𝑀
HH
41
We apply the method above to generate periodic BNP structures in a rectangular prism cell (𝑎 =
100 nm,𝑏 = 80 nm,𝑐 = 200 nm), and a general parallelepiped system with 𝑎 = 100 nm,𝑏 = 80 nm,𝑐 =
120 nm,𝛼 = 80°,𝛽 = 85°,𝛾 = 75°. For consistency with the cubic cell model generated in section 3.2.1,
we choose |𝒒
$
A
| =𝑞
K
=
O√(*(
QK
nm
M(
. In the case of the general parallelepiped system, the application of Eq.
(39) – Eq. (41) to the initial 2000 waves result in a set of 1711 unique 𝒒
$
=(𝑞
$(
,𝑞
$)
,𝑞
$H
) wave vectors with
|𝒒
$
| in a range of values around |𝒒
$
A
|.
To ensure the BNP structure has a well-defined ligament size distribution, we constrain the values
of |𝒒
$
| to a limited range (± 3%) around the desired 𝑞
K
value and only consider the waves with
22

𝜋√151
50
nm
M(
<|𝒒
$
| <
𝜋√171
50
nm
M(

42
The calculated average values of |𝒒
$
| are
O√(*K.SK
QK
nm
M(
and
O√(*( .TQ
QK
nm
M(
for the rectangular
prism and the parallelepiped cells, respectively. Please see Figure 5 for the Fibonacci grid on the surface of
a reciprocal sphere of radius |𝒒
$
A
| and the final distribution of points corresponding to the final set of waves.
To generate a bicontinuous BNP structure with relative density 𝜙
E
= 0.5, we use the level cut 𝜉 = 0. The
resulting nanoporous structures are illustrated in Figure 4(b) and 4(c).
3.3 Structural Analysis
It is instructive to compare the topology of the BNP structures generated using the proposed method
for the rectangular prism and parallelepiped cells to that of the cubic cell generated using the method of
Soyarslan et al.
84
For this purpose, we calculate and compare the ligament size distribution, mean and
Gaussian curvatures, and the genus density of the different topologies.
3.3.1 Ligament Size Distribution
The ligament size distribution is an important characteristic of a given BNP structure. BNP
structures generated by different methods should have comparable distributions of ligament size when
adjusted for different average ligament sizes. We employ the AQUAMI package developed by Stuckner et
al.
97
to evaluate the average ligament size and ligament size distribution. The required binary images for
the AQUAMI processing are generated using the surface mesh tool in OVITO.
98

We calculated the ligament size distribution for the BNP cubic cell produced with the method of
Soyarslan et al.
84
and compare it to the distributions calculated for the rectangular prism and parallelepiped
cells produced with the proposed method. In the calculation, we consider three 2 nm-thick slices of atoms
taken at random positions from each sample. These atoms are used to generate binary cross-section images
using the surface mesh tool in OVITO,
89
setting the mesh color black and background color white. To
further improve the statistics, we generate three independent samples for each type of cell. In total, we
23

process nine binary images with AQUAMI. The algorithm implemented in the software first search and
identify all ligaments and their medial axis. The shortest distance between positions along the ligament
medial axis and ligament surface atoms is noted. Double of the distance is equal to the ligament diameter
or ligament size. The diameter measurement is made on several points in every ligament generating
thousands of measurements from each image. A histogram on ligament size is obtained with 0.22 nm bin-
width, considering the results from all nine images. The average ligament size is the average of all
measurements.
The resulting distribution of ligament size for the three cells is shown in Figure 6. The data show
very good agreement on the average ligament size and ligament size distribution. We find that the data for
all cells follow a Gaussian distribution, with a mean of 4.60 nm and a standard deviation of 1.45 nm. There
is a mild peak at ~8.5 nm. This seems to be a result of thicker sections on branches connecting ligaments.

Figure 6. Ligament size distribution for different periodic cells. A Gaussian distribution is shown to guide
the eye.
3.3.2 Mean and Gaussian Curvatures
Another important characteristic of a BNP structure topology is the average curvature of its surface.
Different methods can be used to estimate the curvature of a BNP topology. We consider the mean and
Gaussian curvatures. As the average curvature of a BNP structure is expected to depend on the relative
densities, we consider the experimental range of relative densities: 0.25 - 0.75. In the numerical calculation,
24

we start by generating a triangular surface mesh for the particular BNP topology. For each surface triangular
element, we calculate the curvatures at its center using the following equations:
99

𝜅
U
=
𝑓
+
)
(𝑓
,,
+𝑓
--
1+𝑓
,
)
(𝑓
++
+𝑓
--
)+𝑓
-
)
(𝑓
++
+𝑓
,,
1
2∙(𝑓
+
)
+𝑓
,
)
+𝑓
-
)
)
H/)
−
𝑓
+
𝑓
,
𝑓
+,
+𝑓
+
𝑓
-
𝑓
+-
+𝑓
,
𝑓
-
𝑓
,-
(𝑓
+
)
+𝑓
,
)
+𝑓
-
)
)
H/)
43

𝜅
V
=2
𝑓
+
𝑓
,
(𝑓
+-
𝑓
,-
−𝑓
+,
𝑓
--
1+𝑓
+
𝑓
-
(𝑓
+,
𝑓
,-
−𝑓
+-
𝑓
,,
1+𝑓
,
𝑓
-
(𝑓
+,
𝑓
+-
−𝑓
,-
𝑓
--
1
(𝑓
+
)
+𝑓
,
)
+𝑓
-
)
)
)
+
𝑓
+
)
(𝑓
,,
𝑓
--
−𝑓
,-
)
1+𝑓
,
)
(𝑓
++
𝑓
--
−𝑓
+-
)
)+𝑓
-
)
(𝑓
++
𝑓
,,
−𝑓
+,
)
1
(𝑓
+
)
+𝑓
,
)
+𝑓
-
)
)
)

44
where 𝜅
U
and 𝜅
V
denote the mean and Gaussian curvatures respectively, 𝑓
$
the first-order derivative of the
implicity surface function 𝑓 defined in Eq. (32) along with the 𝑖 component, and 𝑓
$%
the second-order
derivative. The average mean and Gaussian curvatures of the topology, 〈𝜅
U
〉 and 〈𝜅
V
〉, are then calculated
by the area-weighted average of the 𝜅
U
and 𝜅
V
of all the triangular elements. The calculated data is
displayed in Figure 7. As a reference, we also plot the theoretical curvatures for a random field structure
where the wave vectors are perfectly spherically symmetric. Values for different relative densities are given
by
100

〈𝜅
U
〉 =−√2erf
M(
(𝜙
E
)×
𝑞
K
2
§
𝜋
6
45

〈𝜅
V
〉 =
𝑞
K
)
6
iÈ√2erf
M(
(𝜙
E
)É
)
−1j
46

where "erf" denotes the error function. For comparison, we use dimensionless quantities by normalizing
the mean and the Gaussian curvatures by 𝑞
K
and 𝑞
K
)
, respectively. The data in Figure 7 displays excellent
agreement on the curvature-relative density dependence within the three BNP topologies. An excellent
25

agreement is also indicated between the models and the theoretical prediction. The zeroing of the mean
curvature at a relative density 𝜙
E
=0.5 agrees with the experimental data.
101

3.3.3 Nanoporous Structure Genus
The continuous connectivity of regions is an essential feature of a BNP structure. In order to
characterize the connectivity, we calculate the genus densities with respect to relative density for all the
three different cell types. The genus of a connected surface, 𝑔, can be expressed as
32

𝑔 =1−i
𝜒
2
j 47
where the Euler Characteristic, 𝜒, an integer associated with the topology of the given shape, is defined as
𝜒 =𝑣−𝑒+𝑓 48
where 𝑣, 𝑒, and 𝑓 denote the number of vertices, edges, and faces of a polyhedral or meshed surface. To
facilitate the comparison of the genus of BNP structures of different relative densities, we use the scaled
genus density, 𝑔
W
ÌÌÌÌ, defined as
𝑔
W
ÌÌÌÌ =
𝑔
𝑉
𝑆
W
MH
49
where 𝑉 is the total volume of the simulation box and 𝑆
W
the specific area of the foam solid (solid area/solid
volume). We calculate and plot in Figure 8 the 𝑔
W
ÌÌÌÌ for BNP structures with relative densities in the range of
0.25<𝜙
E
<0.55.
26

Figure 7. Numerical measurements of mean and Gaussian curvature plots as a function of relative density
for different periodic cells. Continuous lines correspond to theoretical predictions. The displayed values for
mean and Gaussian curvature are scaled by dividing 𝑞
K
and 𝑞
K
)
respectively for both numerical
measurements and theoretical predictions.
The numerical results indicate consistent values for the cubic, general parallelepiped, and
rectangular prism cells. For comparison, we also include reported values of 𝑔
W
ÌÌÌÌ of experimental samples,
102–
104
and foams generated using phase-field
102
and MD
32
simulations. Overall, the experimental and numerical
values agree well and follow the theoretical 𝑔
W
ÌÌÌÌ
X
curve. The values of 𝑔
W
ÌÌÌÌ
X
are calculated assuming an
infinite number of waves, homogenously distributed over the surface of a spherical volume of radius 𝑞
K
.
𝑔
W
ÌÌÌÌ
X
is calculated from the theoretical specific area, 𝑆
W
X
, and genus density, 𝑔
W
X
, which have been derived
previously
84
and are given by
𝑆
W
X
=−
1
𝜙
E
2𝑞
K
𝜋√3
𝑒
Y
&
)
50
𝑔
W
X
=−
1
4𝜋
〈𝜅
V
〉𝑆
W
X
51
Combining Eq. (46) and Eq. (51), 𝑔
W
ÌÌÌÌ
X
is given by
𝑔
W
ÌÌÌÌ
X
=𝜙
E
H
𝜋
32
(1−𝜉
)
)𝑒
Y
&
52
where 𝜉 is a function of 𝜙
E

𝜉 =√2erf
M(
(2𝜙
E
−1) 53
27

One should note that the numerical values 𝑔
W
ÌÌÌÌ for different cells are calculated using a finite mesh of the
BNP structure. The deviation of 𝑔
W
ÌÌÌÌ from 𝑔
W
ÌÌÌÌ
X
at large relative densities can be minimized by using a finer
mesh of the surface.
3.4 Extension of the Method: Gradient Nanoporous Structure
A promising new design combines BNP materials of different length scales to generate seamless
gradient nanoporous (GNP) materials.
87,88
These designs may be used to explore the combination of
properties offered by BNP materials with different pore sizes. GNP designs have the potential to induce
unpredictable properties originating from the non-trivial interaction of distinct structural elements displayed
by the gradient structure. Recently, preliminary GNP metal samples were synthesized using two techniques:
1) free corrosion segmental dealloying at different temperatures
88
and 2) free corrosion dealloying of
samples with gradient composition.
87
The results have indicated that GNP metals can be used effectively in
surface-enhanced Raman scattering surpassing the performance of nanoporous metals suggesting many
applications, e.g., hierarchical “nature-like” bone graft coating for orthopedic implants, gas diffusion layers
for fuel cells, and catalysis.

Figure 8. Scaled genus density versus relative density for different periodic cells. The line represents the
theoretical prediction given by Eq. (52). Results from MD,
32
Phase Field,
102
and experimental
reports
102,103,105
are shown for comparison.
A completely unexplored area is the use of GNP metal designs to generate materials with novel
mechanical properties. Our idea here is to harness the possible benefits of combining gradient and NP
28

structures. For example, gradient grain-sized designs have been shown to generate materials with
unprecedented mechanical properties, displaying an outstanding combination of strength, ductility, and
hardness.
106–111
Furthermore, BNP materials have demonstrated a strength/density ratio that was much
higher than the predictions from macroscopic theories. Therefore, GNP metals could leverage the high
strength/density ratio of BNP metals with the non-trivial effects originating from the gradient structural
design to potentially: 1) generate materials that could exceed the reported strength-to-density ratio value of
nanoporous metals; 2) generate materials with better stability and toughness than displayed by brittle
unimodal BNP metals.
The method proposed here can generate a GNP structure for atomistic simulations with slight
modifications. To generate a gradient of ligament size in the 𝑥-direction, two different sets of 𝒒 are used,
𝒒
Z[\]]
and 𝒒
]\:^C
. 𝒒
Z[\]]
indicates that the corresponding wavenumber |𝒒
Z[\]]
| is small while that of
2𝒒
]\:^C
2 is large. The fields generated by such two sets of waves are:
𝑓
]\:^C
(𝒓)=¡
2
𝑁
]\:^C
I cos(𝒒
𝒍𝒂𝒓𝒈𝒆
∙𝒓+𝜑
]\:^C
1
!
'()*+
$;(
small ligament
54
𝑓
Z[\]]
(𝒓) =¡
2
𝑁
Z[\]]
I cos(𝒒
𝒔𝒎𝒂𝒍𝒍
∙𝒓+𝜑
Z[\]]
)
!
,-(''
$;(
large ligament 55

The gradient is created by weighted summation, in which the weight of two fields is related to the 𝑥
coordinate:
𝐹(𝒓) =i1−
𝑥
𝑎
j×𝑓
]\:^C
(𝒓)+
𝑥
𝑎
×𝑓
Z[\]]
(𝒓)
56
where 𝑎 is the length of the simulation box. The structure generated is shown in Figure 9.
29

Figure 9. GNP Structure.
One can note that the structure shown in Figure 9 preserves the periodicity in the 𝑦 and 𝑧-directions, while
along the 𝑥- direction the periodicity is broken to allow for the gradient in ligament size to be defined. To
preserve PBC in the 𝑥-direction, a symmetric gradient can be designed, which results in a small ligament
size at x = 0, growing and reaching the largest value in the middle, and returning to the starting level at x =
a. The corresponding field is
𝐹(𝒓)=Ð
b1−
2𝑥
𝑎
d×𝑓
]\:^C
(𝒓)+
2𝑥
𝑎
×𝑓
Z[\]]
(𝒓) 𝑥 <𝑎/2
b
2𝑥
𝑎
−1d×𝑓
]\:^C
(𝒓)+b2−
2𝑥
𝑎
d×𝑓
Z[\]]
(𝒓) 𝑥 >𝑎/2

57
Eq. (57) creates the double gradient structure shown in Fig. 10.

Figure 10. Double GNP.
3.5 Conclusions
In this work, we proposed an efficient method to generate periodic BNP structures in general
parallelepiped cells of different length sizes and shapes. The relative density of the structures is defined by
choosing a level set value. The calculated mean and Gaussian curvatures in periodic cubic, rectangular
30

prism, and general parallelepiped cells are in excellent agreement with theoretical predictions. In addition,
the calculated dependence of the scaled genus density with relative density agrees well with reported
experimental and numerical data using different methods. Compared to the commonly used phased field
simulation of the spinodal decomposition process, this method offers an attractive alternative for the
generation of BNP models for atomistic simulations.

31

Chapter 4
Bicontinuous Nanoporous Design Induced Homogenization of Strain
Localization in Metallic Glasses
4.1 Background
MG has unique and contrasting mechanical, chemical, and physical properties when compared to
their crystalline counterparts.
112–116
The combination of such properties with the high surface-to-volume
ratio, i.e., specific area, produced with the introduction of a high porosity level in the material makes porous
MGs (PMG) promising candidates for lightweight structural applications, membranes, porous electrodes,
and catalysts in a large range of fields including aerospace, medicine, electrochemistry, and environmental
protection.
15,117–123

While research on PMGs has suggested many applications, its functional use, e.g., in catalysis and
hydrogen storage, demands a much larger specific surface area than currently available. The morphology
of BNP structures, with their inherent high specific surface area, offers an elegant path to enhance the
functional capabilities of PMGs. Recently, Jiao et al.
40
successfully fabricated a 3-dimensional BNPMG by
selective corrosion of spinodally decomposed MG precursors and passivation. The stochastic nature of the
spinodal process and the nano-length scale of the porous features result in a high specific surface area
making such material an excellent candidate for gas absorption.
40

Besides providing a high specific area, porous designs have also been used to alleviate the strength-
ductility trade-off that significantly limits the application of MGs.
118,124
Sopu et al.
39
demonstrated that
Cu 64Zr 36 NPMG with 3% porosity displays enhanced plasticity with an optimized heterostructure. It is
worth noting that the newly developed BNP structures are morphologically different from traditional closed
and open-cell PMGs previously investigated. It would be interesting to study the mechanical properties,
which closely depend on the structural and morphological features.
32

In this work, we generate Cu 64Zr 36 BNPMG using the level-set method presented in Chapter 3. We
then perform tensile loading MD simulations to describe the deformation and failure processes and their
underlying mechanisms.
4.2 Method
4.2.1 Model Generation
To prepare the samples, initially, a CuZr MG is generated following a similar procedure as reported
previously.
125
A cubic simulation cell with a 5 nm side is filled with 9,826 atoms (~64% Cu and ~36% Zr)
in an FCC crystal structure. PBC is applied along the three cartesian directions. The system is then heated
up and thermalized at 2,000 K and 0 GPa for 2 ns. Nose-Hoover thermostat and barostat are applied to
control temperature and pressure. The thermalized liquid is then quenched to 50 K using a -0.01 K/ps
cooling rate. The small resulting amorphous cube is then replicated 10 times along the 𝑥, 𝑦, and 𝑧-directions.
The large system containing 1,000 small cubes is then relaxed at 800 K for 0.250 ns to minimize periodic
patterns in the sample. The final MG model is obtained after subsequent quenching to 50 K using a -15
K/ps cooling rate.
To generate a BNP structure based on the MG sample we use the level set method described in
Chapter 3 and reported in the literature.
126
We use a Fibonacci grid to generate 3,000 wave vectors with
wavenumber
!𝑞
0
! = 2𝜋
√40
𝑎

58
where 𝑎 = 537.941 nm is one side of the simulation box. After applying PBCs and subsequently eliminating
repeated wave vectors, 644 waves remain. Random phases are then attached to each wave. For each atom,
the value is calculated using Eq. (28). The level set is chosen to be 𝜉 =0.138, i.e., atoms with 𝑓 <𝜉 are
deleted. Taking three 2-nm slices at random positions in each direction, we employ the AQUAMI package
97

to evaluate the average ligament size and ligament size distribution based on suitable images generated
using the surface mesh tool of OVITO.
89

33

4.2.2 Measurement of Surface Curvature
The surface mesh generated from the entire solid phase is then used to evaluate surface curvatures
throughout the structure. The mean and Gaussian curvatures, 𝜅
U
and 𝜅
V
, are analytically calculated on
every surface mesh point following Eq. (43) and Eq. (44). The two principal curvatures, 𝜅
(
and 𝜅
)
, are then
calculated as
𝜅
1
= 𝜅
𝐻
−$|𝜅
ℎ
2
−𝜅
𝐺
|
59

𝜅
!
= 𝜅
𝐻
+$|𝜅
ℎ
2
−𝜅
𝐺
|
60
4.2.3 Simulation Details
After the generation, the BNPMG is relaxed at 50 K for 0.6 ns. The structure is then deformed
under uniaxial tensile loading along the 𝑥-direction at a constant 3×10
8
s
M(
strain rate. The temperature is
maintained at 50 K and the stress on the (𝑦 and 𝑧) perpendicular directions is kept at zero. The atomic-level
stress and strain are calculated by LAMMPS
127
based on locally averaged Virial stress
128
and OVITO based
on von Mises shear strain,
129
respectively. To resolve the normal direction throughout the surface we resort
to the local atomic data of the BNPMG system, which is separated into the surface and internal sets. The
criteria used to classify atoms in “surface” and “internal” sets are based on their local density calculated by
atomic Voronoi analysis. Surface atoms have significantly larger atomic volumes (lower number density)
compared to internal atoms. For each surface atom, we construct a set of internal atoms searching within a
cutoff distance of 6 Å. The normal direction at the position of each surface atom is then calculated as the
sum of the vectors from every internal atom to the surface atom considered.
All MD simulations used in the generation of the simulation model and to perform tensile loading
tests are performed with LAMMPS. The interatomic potential proposed by Mendelev et al.
130
is used to
describe interactions between atoms. The potential parameters were based on a set of experimental data on
34

atomic sizes, lattice parameters of different structures, cohesive energy, melting temperature, latent heat of
melting, and other properties.
130
In addition, this potential is able to predict the Cu-Zr amorphous alloy
structure in reasonable agreement with experiments. The stress is calculated based on a general formulation
of stress tensor for arbitrary many-body interaction potentials under PBCs proposed by Thompson et al.
119

All simulations use a constant integration time step of 1 fs.

Figure 11. Morphology of the BNPMG sample investigated. (a) Illustration of as-generated Cu 64Zr 36
BNPMG sample and (b) its ligament size distribution, which follows a Gaussian distribution with a mean
value of 4.4 nm. (c) The probability distribution of principal curvatures, 𝜅
(
and 𝜅
)
, from the surface of the
BNPMG sample.
35

4.3 Results
4.3.1 Analysis of Bicontinuous Topologies
The initial relaxed BNPMG structure is shown in Figure 11(a). The structure illustrated highlights
the bicontinuous morphology features, which display an open porous structure with randomly oriented
ligaments. The specific surface area is 𝑆
l
= 0.0269 Å
M(
(269 m
2
⁄cm
3
). The calculated ligament size
distribution of the sample is shown in Figure 11(b). The distribution, in the range from 0.8 nm to 10.2 nm,
is fitted with a Gaussian curve displaying a 4.44 nm mean value.
The morphology of the initial structure is quantified by calculating the statistics of principal
curvatures, 𝜅
1
and 𝜅
2
, over the BNPMG surface. The data is shown in Figure 11(c), where the curvatures
are normalized by the specific surface area 𝑆
l
, making them dimensionless. The distribution is plotted using
a bicubic interpolation of the data. Analysis of Figure 11(c) indicates that the surface is dominated by saddle
points located around the purple region, where 𝜅
1
= −𝜅
2
. The porosity level, higher than 50%, makes the
structure predominantly convex, i.e., the distribution has a higher probability near the red box region in
Figure 11(c) than near the blue box region, which is linked to a concave surface.
4.3.2 Mechanical Response to Tensile Loading
The mechanical behavior of the BNPMG sample is investigated by performing tensile loading
simulations. The resulting stress-strain curve is shown in Figure 12(a) and indicates a rather unusual
mechanical behavior for an MG. At a low density of 3.47 g/cm
3
, the yield strength and ultimate tensile
strength are about 400 MPa and 430 MPa, respectively. The sample displays a contrasting mechanical
behavior compared to typical bulk MGs at the yield point, i.e., it displays a short and flat plastic region
from e=0.06 to e=0.1 as shown in Figure 12(a). Failure occurs graciously as suggested by the gradual
stress drop in the region from e=0.1 to e=0.4. In this work, failure is defined as the point where the
material loses mechanical integrity and sharply releases the loading stress.
36

Figure 12. Tensile loading global and local stresses evolution during deformation and failure. (a) Stress-
strain curve displaying well-defined elastic, plastic, and failure regimes. The sample yield and ultimate
tensile strengths are ~400 MPa and ~430 MPa, respectively. The curve indicates an atypical ductile
behavior for an MG. (b) Illustrations of the BNPMG sample at six strain levels highlighting the local stress
distribution during the plastic deformation and failure. Atoms are colored by the component of the stress
along the loading direction averaged locally.
To rule out possible artifacts brought by the high strain rate deformation and size effects
131
during
simulations, we also consider the tensile loading mechanical behavior of two additional samples with aspect
ratios of 2:1 and 4:1. These two samples are generated by duplication of the cubic sample along the loading
direction. The corresponding curves presented in Figure 12(a) reveal that the high-aspect-ratio samples
37

have identical elastic behavior in the strain range up to e=0.06, and very similar plastic response in the
strain range 0.06< e<0.16. The stress then drops sharply for the high-aspect-ratio samples, which indicates
the fracture of the samples carried out by a cascade of ligaments rupture. By comparing the behavior of the
three samples, one can note that the long tail of the stress-strain curve of the 1:1 aspect ratio sample is an
MD artifact and not a material property. The plastic region of the nanoporous samples based on the results
of Figure 12 is estimated to be from e=0.06 to e=0.16.
In order to understand the unusual mechanical response, we investigate the distribution of local
stress during the deformation and failure processes. The atomic stress along the loading direction is
averaged locally in spherical regions of radius 11 Å. Illustrations of the deformation and failure processes
are shown in Figure 12(b), where the atoms are colored according to the value of the local stress. Six frames
are shown in Figure 12(b), corresponding to the A-F strain values indicated in Figure 12(a). Since the
surface-to-volume ratio is large for nanoscale materials such as nanowires and nanoporous structures,
surface stress may have a significant influence on the mechanical properties. For both crystalline and glassy
metallic nanowires, the presence of a stress gradient has been shown in previous studies,
132,133
i.e., the stress
state shifts from tensile on the surface to compressive in the center of the nanowires.

Figure 13. Snapshot of the sample at 0 tensile strain. Colors are assigned to atoms based on the local stress
level.
38

Figure 14. Probability distribution of the components of a unit normal vector on the surface of the sample
at three strain levels. Horizontal axes represent the component parallel to the loading (𝑥 axis), while vertical
axes represent the two perpendicular components (𝑦 and 𝑧 axes). The undeformed structure (e = 0) implies
a rather isotropic distribution, indicated by the non-biased distribution of normal components, implying a
random alignment of the ligaments. During the loading, the probability of high parallel component drops
sharply, indicating the progressive alignment of ligaments along the loading direction during the tensile
loading.
Asymmetry in the mechanical behavior of nanowires with very small diameters under tensile and
compressive loading was suggested to be caused by the significant surface stress in these nanostructures.
134–
136
Recently, this tension/compression asymmetry was also reported for nanoporous crystalline and glassy
metals.
137,138
The surface stress is expected to have a larger effect going from MG to crystalline metals since
the atoms in MG structures can rearrange their local structure in order to relax the stress more effectively
than atoms constrained to a crystalline structure. A snapshot of the as-relaxed structure of bicontinuous MG
is presented in Figure 13, indicating a gradual change from high tensile to compressive stress when moving
39

from ligament surface to internal regions. Atoms with different stress states balance each other, resulting in
overall zero stress at the beginning of the loading simulations.
At the end of the elastic region, multiple stress concentration spots can be observed. Analysis of
the structure indicates that stress is concentrated in the narrow cross-sections of the ligaments mostly
aligned with the loading direction. Due to the stochastic nature of the structure, those regions are spread
throughout the whole structure. The locally stressed ligaments develop necking and eventually fail as the
average tensile strain of the porous structure increases. The homogeneous distribution of local ligament
necking throughout the system and their gradual development and rupture results in an increase in ductility.
4.3.3 Realignment of Ligaments Towards Loading Direction
The progressive alignment of the ligaments before their deformation and rupture is an important
result of this study. To quantify the alignment of the ligaments during the tensile loading process, we turn
to the probability distribution of the surface normal direction, shown in Figure 14. We scan the surface of
the porous structure and use the normal direction at each surface point to calculate the probability
distribution of the normal direction components decomposed along the loading direction, parallel ||, and
along the two other perpendicular directions, ^1 and ^ 2. In the beginning, the ligaments are randomly
distributed, and the distribution is isotropic. As the tensile loading progresses, Figure 14(c) andFigure 14(d)
indicate that the normal direction aligns along the ^1 and ^2 directions, i.e., the probability for the parallel
component becomes concentrated around zero. These results quantify to a certain extent the significant
alignment of ligaments along with the loading direction, especially in the elastic region for strains up to
e= 0.06. This indicates that the ligaments bend significantly in the elastic region prior to deformation by
necking. Interestingly, Sun et al. also observed ligament alignment towards loading direction during tensile
loading of crystalline nanoporous gold.
10

40

4.3.4 Analysis of the Atomic Local Strain during Tensile Loading

Figure 15. Strain localization and failure of BNPMG. (a) Fraction of atoms in plastically deformed regions
- von Mises shear strain in the range, n = 0.2 - 0.7, vs. the fraction of atoms in severely deformed regions
with n > 0.7. Only 2% of atoms are in severely deformed regions, when 6% of the atoms are participating
in plastic deformation. (b) Illustrations of the BNPMG sample at six strain levels during the loading with
atoms colored by their von Mises shear strain. High strain atoms are concentrated in necking regions that
are distributed over the whole system. The delocalized necking gives rise to the ductile behavior displayed
by the BNPMG.
41

To quantify the plastic deformation process, we calculate the atomic local von Mises shear strain,
129

n, which is useful to help understand the strain localization processes that often trigger the failure in MGs.
139

In Figure 15(a) we identify strain localization in the system by plotting the fraction of atoms in severely
deformed regions, n 0.7, i.e., the fraction of atoms with n>0.7, against the fraction of atoms in plastically
deformed regions, n 0.2-0.7, i.e., the fraction of atoms with 0.2> n >0.7. One can note that the n 0.7 fraction
of atoms does not increase drastically until the overall tensile strain reaches e = 0.15, indicated by the tag
“D”. At this point, the n 0.2-0.7 fraction is 5.5% indicating that strain localization is delayed until the structure
has undergone noticeable and delocalized plastic deformation. This behavior contrasts with that typically
observed in bulk MG, where shear bands dominate the plastic deformation and the failure, which occurs
swiftly after the yield point. In comparison, bulk MGs display a very small fraction of atoms in plastically
deformed regions as compared to atoms in severely deformed regions.
139,140

To further illustrate the delocalized deformation and failure process responsible for the rather
ductile stress-strain curve shown in Figure 12(a) we display in Figure 15(b) snapshots of the BNPMG
structure, color-coded by n at different strain levels as indicated in Figure 12(a) and Figure 15(a). One can
note that delocalized necking events and other regions undergoing plastic deformation are highlighted by
the contrasting large n values of their respective atoms.
From a close examination of the deformation of single ligaments, we find that the fundamental
mechanism of necking and rupture of ligaments is viscous flow. Consistent with the trends shown in Figure
15(a), the distribution of atoms participating in deformation and severe deformation is rather delocalized
throughout the structure. The overall picture of the deformation and failure indicates a quite unusual
deformation and failure process among different classes of materials since localization and delocalization
of plastic events synergize to produce a mechanical response that is unseen for MGs. This combination of
deformation mechanisms in the BNPMG promotes the synergy between strength and ductility.
42

4.4 Discussion
4.4.1 Local Heating
Experimental data from in-situ tensile loading on Cu 47Ti 33Zr 11Ni 6Sn 2Si 1 bulk MG (BMG) has
indicated that local heating
105
can be generated at shear band events, resulting in local melting and
nanocrystallization.
141
Figure 16 reveals the local temperature profile of the BNPMG at strain e=0.12.
While in our study, the failure of the ligaments is not dominated by shear-band-like events, localized plastic
deformation by necking does occur. Nonetheless, the local temperature at localized plastic deformation
regions does not exceed 150 K. Therefore, the local temperature increase does not lead to melting or
recrystallization.

Figure 16. The local temperature in a deformed BNPMG at 𝜀 =0.12.
4.4.2 Comparison to Nanoporous Metallic Glasses with Simple-Shaped Pores
In contrast with the BNP structure considered in this work, much of the research effort in this area
has been on porous structures with simple, well-defined, geometrical pores.
39,142,143
Several factors are
reported to affect the ductility of NPMGs, including porosity, pore size, pore shape, and the arrangement
of pores. Ductility improvement is linked to changes in plastic deformation mode from single localized
shear banding to delocalized shear banding to necking-like homogeneous flow. Compared to the BNPMGs
studied here, in general, NPMGs display larger strengths and lower ductility. This is a direct consequence
43

of the differences in porosity and morphology of the porous structures. Recently, Lin et al.
137
investigated
the mechanical behavior of open-cell Cu 50Zr 50 NPMG structures, which are comparable to those considered
here, although with well-defined ligament shapes and dimensions. They found that ligaments with a low
slenderness ratio (aspect ratio) or under a combined shear/tensile stress state tend to form a shear band
instead of necking. This is consistent with our simulation results since, in our sample, shear bands can be
found in the joint of ligaments and some thick ligaments. Nonetheless, they reported a similar ductile
behavior mechanism dominated by strain delocalization and necking of ligaments.
4.4.3 Effect of Metallic Glasses Ligaments
The mechanical behavior of BNP materials encompasses contributions from both material response,
inherent to the deformation of each ligament in the structure, and structure response, coming from the
collaborative deformation mode of the ligaments. The latter also includes ligament alignment throughout
the structure. To verify the structure response, we perform additional similar tensile loading simulations on
cylindrical ligaments of 20 nm length and diameters of 3 and 7 nm. 3 nm ligament displays outstanding
ductility while 7 nm ligament fails right after the yield point. These ligaments are representative of typical
ligaments in the BNP structure considered in this work. Figure 17 demonstrates the stress-strain curves of
the individual ligaments and the BNP system. The stress for the BNPMG is scaled by the inverse of the
square of its relative density for better comparison. The stress-strain curves of the ligaments indicate the
underlying material's response present in the BNP sample rather than the structure response. The curves
indicate that the yield points for the ligaments - e= 0.06 and e= 0.07 - are in agreement with that of the
BNP system. The 7 nm ligament displays a sharp drop of the stress at the yield point, indicating the failure
of the system. In contrast, the 3 nm diameter ligament displays an extended plastic region up to e = 0.11,
delaying the failure of the system. Figure 18 indicates snapshots of the two ligaments at an engineering
strain e = 0.12 highlighting their necking failure mode. After the sharp drop of stress, the tail of the stress-
strain curves is a result of constant strain rate and low aspect ratio of the samples in the MD simulations.
44

Figure 17. Stress-strain curves for the BNPMG sample and MG ligaments of different diameters. The stress
for the BNP sample is normalized by the square of the relative density of the sample, r
)
.
Both ligaments have the same stiffness, and it is higher than the BNP system. Combining with the
significant normal distribution change shown in Figure 14 in the strain range e = 0 to e = 0.06, the slightly
elongated elastic region of the BNP sample can be partially linked to ligament reorientation, which delays
the stretching and rupture of the ligaments and results in a lower stiffness because of the bending
deformation. However, the fraction of the ligaments initially aligned with loading defines the yield point.
4.4.4 Comparison to Bicontinuous Nanoporous Crystalline Metals

Figure 18. Snapshots of the failure by necking process of the MG ligaments.
Compared to BNPMG, BNP crystalline metals have been investigated extensively.
144
Necking of
ligaments, which is linked with ductile behavior has also been reported in BNP metals at ligaments with
diameters around 100 nm.
18,145
Nonetheless, BNP metals commonly display a lack of macroscopic ductility
and fail by brittle fracture compared to bulk metals.
145
One of the pioneer studies of BNP metals indicates
45

that such material becomes extremely brittle when the relative ligament size compared to the sample size
is small.
25
Interestingly, atomistic modeling investigations suggest larger ductility of BNP metals than
indicated by experimental data.
10,146
Regardless, the failure of BNP metals is linked to the rupture of
ligaments in the structure, starting at the weakest ligament, leading to catastrophic failure by overload and
rupture of adjacent ligaments.
145
The level of ductility shown in our results for the BNPMG is comparable
to that of BNP metals reported in atomistic modeling investigations, displaying a failure strain of 0.3 or
less.
10,146
The overall strength of BNPMG presented here is higher compared to a crystalline metal
counterpart, e.g., the ultimate tensile strength reported range is 260 MPa to 320 MPa for BNP metals at
comparable ligament size and porosity.
10,138,146,147

4.4.5 Bicontinuous Nanoporous Metallic Glass vs. Collection of Metallic Glass
Nanopillars
In the absence of pores, BMGs are generally very brittle and display limited plasticity in the form
of highly localized shear-banding in both tensile and compressive loading.
125,148,149
Ductility can be
improved as well by exploring size effects in the mechanical properties of MGs.
150,151
Maybe the most
dramatic size effects are those displayed in MG nanowires and nanopillars that undergo a brittle-to-ductile
transition triggered by size reduction. This size effect is found under tensile as well as under compressive
and bending loading.
152,153
The BNPMG discussed in this work can be seen as a collection of randomly
oriented MG nanowires. However, the ligaments have irregular cross-sections, often having a narrow
diameter at their mid-length. That prevents an accurate assessment of their aspect ratios. An alternative
BNP model with cylindrical ligaments would allow such analysis.
137
Results indicate the displayed
homogeneous deformation is carried out by those ligaments, spread across the loading direction, which
delocalizes the BNPMG deformation by homogeneously distributed local necking deformation, resulting
in an overall large ductility of the bulk sample. In the framework discussed by Greer and De Hosson,
153
one
can say that the extrinsic size effect on the mechanical properties of the ligaments is dictating the intrinsic
mechanical behavior of the BNPMG. The estimated aspect ratio of the ligaments in the bicontinuous
46

structure depends closely on the porosity level. As porosity increases, the surface curvature will concentrate
less in the saddle points areas and more in the cylindrical areas than what is shown in Figure 11. Thus, the
estimated aspect ratio of the ligaments increases. An investigation of the dependence of the mechanical
properties of the BNPMG as a function of porosity is left for future work. Following previous studies on
nanoporous crystalline metals, we expect that the strength will drop dramatically with increasing porosity.
However, it is still to be determined how porosity affects the overall ductility of the BNPMG. By doing
microbending experiments on MG nanopillars,
105,153
it was found that the critical stress required to create a
shear band under bending was twice as large as that under compression/tension. This supports our results,
which indicate localized plasticity for ligaments under tension and few shear bands or necking events for
ligaments under bending, resulting in an increase of ductility in the BNPMG sample.
4.4.6 Effect of High Strain Rate Used in Molecular Dynamics Simulations
It is important to mention the possible effects caused by the high strain rate used in this work. The
strain rates typically used in MD simulations are in the range from 10
7
to 10
10
s
-1
. That is several orders of
magnitude larger than the strain rates typically used in experiments. However, it has been shown that MD
simulations using high strain rates can reproduce critical features of the mechanical behavior of BMGs and
MG nanowires and nanopillars. That includes the formation of shear bands, brittle failure under tensile
loading,
149,154,155
size induced brittle-to-ductile transition, and formation of necking.
140,156,157
Cao et al.
156
argue that the high strain rate effect in MD simulations is offset by the high cooling rates used in the
generation of the MG models. These previous studies support our results as the BNPMG presented in this
work is a bulk system consisting of a collection of MG nanowires. Although the bicontinuous design can
delocalize the stress to some degree, stress concentration still exists and dictates the failure of the structure.
As one ligament fails, the stress on the adjacent ligaments is intensified. The absence of experimental results
on BNPMG prevents us from further validating our simulation results. However, a similar cascading failure
of ligaments was observed for crystalline BNP metallic systems.
18,145
Nonetheless, the profile of local
instability and fracture of the BNP structure may have a strain rate dependence, which may affect the
47

number of shear bands and necking events observed. An investigation of the strain rate dependence is left
for a future study.
4.5 Conclusions
This study implies that ligaments initially aligned with the loading direction bear most of the load
and are responsible for the plasticity generated immediately beyond the yield point. Ligaments in other
orientations also contribute to the overall behavior of the BNPMG as they bend over ligament joints and
progressively align with the loading direction. The reorientation of ligaments occurring mostly in the elastic
region can contribute notably to the enhanced ductility of BNPMG as it occurs concomitantly to the necking
of aligned ligaments. As ligaments align with the loading direction, they are able to support more load,
which is transferred from ligaments with fully developed necking. The synergistic combination of randomly
oriented ligaments and the ductile nature inherited from their thin cross-sections results in an anomalous
ductile behavior for a nominally brittle MG composition. These findings shed light on the mechanical
behavior of BNPMGs and are expected to motivate further experimental and modeling research, which will
pave the way for practical applications of these materials.
48

Chapter 5
Structural and Compositional Effects on the Mechanical Properties
and Scaling Laws of Bicontinuous Nanoporous Metallic Glasses
5.1 Background
While monolithic MGs have well-defined mechanical properties rooted in their materials properties,
those of NPMGs also depend on their porous structure topology. Efforts have been spent in developing
scaling laws for NPMGs based on experimental results and theoretical calculations.
10,103,146,158–162
Current
scaling laws are developed empirically or based on simplified topology models that consider only porosity
as a variable. Complex functional relationships depending on multiple variables, such as porosity and
composition, can be developed using Artificial Intelligence (AI) based methods. The emergence of AI
methods is bringing a new dawn to the development of materials science.
163,164
Machine Learning (ML) and
Genetic Programming (GP) algorithms have been experiencing a resurgence in recent years within the
materials science community.
165,166
AI technologies can provide researchers with tools to overcome the
barriers between designing, synthesizing, and processing materials by accelerating simulations,
167,168

predictions of properties,
162,169–172
design of synthetic routes,
173,174
optimizations of experimental
parameters,
175,176
and enhancement of characterization methods.
177,178
As for nanoporous materials, it is
challenging to predict the mechanical properties due to their complex geometry. A suitable AI method, such
as GP, which mimics natural selection, can be employed to describe the functional relationship between
bicontinuous NPMGs mechanical properties and system variables.
Motivated by Darwin’s theory of natural selection, GP is a powerful evolutionary technique
commonly used to automatically generate programs suitable for user-defined tasks.
179
AI technologies are
currently enabling the uncovering of physical laws. Among the AI methods, GP is a promising technique
to accomplish symbolic regression allowing one to find suitable mathematical models describing data when
little knowledge of the data structure or distribution is available.
180,181
Within the GP-assisted symbolic
49

regression (GPSR) expressions are randomly generated and evolve in successive generations improving the
description of the relationships of interest.
182
GPSR has already been applied to numerous problems.
183
In
contrast to other AI methods,
184–186
one can derive physical insights from the expressions obtained with the
GP method. Chopra et al.
187
developed GP models to predict the compressive strength of concrete based on
in situ data from the literature. Cai et al.
188
used GPSR to derive heat transfer correlations, including the
functional equation and its parameters, from experimental data on heat transfer measurements, which were
used to predict the performance of thermal components. Langdon and Barrett
189
applied GP in drug
discovery by evolving simple, biologically interpretable, in silico models of human oral bioavailability.
Barmpalexis et al.
190
found a function mapping levels of four polymers to three different properties of a
pharmaceutical release tablet with the help of GPSR. Recently, Im et al.
191
applied the GP methodology to
identify governing equations in non-linear multi-physics systems.
In the present work, we use MD simulations to study the mechanical properties of CuZr BNPMGs.
We then apply GPSR to derive scaling laws describing the mechanical behavior as a function of system
variables and compare them against existing scaling relationships. The results indicate that the GPSR-
derived models are able to predict accurately both the Young’s modulus and ultimate strength as a function
of relative density and alloy composition. This work demonstrates that the GPSR is able to uncover
expressions that can predict accurately the mechanical properties and also provide physical insights into
complex systems.
5.2 Method
5.2.1 Model Generation and Morphology Characterization.
Initially, cubic CuxZr 1-x BMG simulation cells with 5 nm sides are prepared following the same
procedure described in section 4.2.1. Five different compositions, x = 0.25, 0.36, 0.5, 0.64, and 0.72 Cu are
selected to study the effect of compositions on the mechanical behavior of BNPMG. We chose this
composition range based on experimental reports on the stability of CuZr BMGs alloys.
192
PBCs are applied
along the three cartesian directions. The small BMG cubic system is then replicated multiple times along
50

the 𝑥, 𝑦, and 𝑧 directions to obtain large BMG samples. To study the effect of ligament size while
controlling the porosity, we create Cu 64Zr 36 BMG samples with different side lengths varying from 40 nm
to 80 nm. For all the other compositions, we employ a side length of 54 nm. The large systems containing
around 4,864,000 – 32,062,500 atoms are then relaxed at 800 K for 0.25 ns to minimize periodic patterns
in the sample. The final BMG models are obtained after subsequent quenching to 50 K using a -15 K/ps
cooling rate.
We use the level set method reported in Section 4.2.1 to generate a BNP structure based on the
BMG sample. Different level sets are then applied to generate samples with various relative densities. In
order to ensure the same morphology for all the samples with different side lengths, we keep the product of
the wavenumber 𝑞
K
and the side length 𝐿 constant:
𝑞
K
𝐿 =𝐶
K
61
where 𝐶
K
denotes a constant. The morphologies of the samples with different porosities are shown in Fig.
19.

Figure 19. Initial sample structures with relative density 𝜌 = 0.3, 0.4, 0.5, 0.6, and 0.7 with corresponding
surface curvature distributions. Colors in the lower plots indicate the probability density. Yellow, red, and
purple boxes indicate regions with concave, convex, and saddle surfaces. Dashed lines indicate opposite
values of the principal curvatures while the dotted lines indicate identical values.
Taking three 2-nm thick slices at random positions along each direction, we employ the AQUAMI
package
97
to evaluate the average ligament size and ligament size distribution for samples with a side length
51

of 54 nm based on suitable images generated using the surface mesh tool of OVITO.
89
The ligament size
distribution is displayed in Fig. 20 and indicates that the average ligament size varies from ~3.5 to ~6.2 nm.
5.2.2 Simulation Details
After the generation of the BNP structure, the NPMG samples are relaxed at 50 K for 0.6 ns. The
structure is then deformed under uniaxial tensile loading along the 𝑥-direction at a constant 3×10
8
s
M(
strain
rate. The temperature is maintained at 50 K and the stresses on the perpendicular directions to the loading,
i.e., 𝑦 and 𝑧-directions, are kept at zero. The atomic-level stress and strain are calculated using the same
methods as described in section 4.2.2.
All MD simulations used in the generation of the simulation models and to perform tensile loading
tests are performed with LAMMPS. The interatomic potential proposed by Mendelev et al.
130
is used to
describe interactions in the Cu-Zr binary alloy. All simulations use a constant integration time step of 1 fs.

Figure 20. Ligament size distributions for the initial structures with relative density 𝜌 = 0.3, 0.4, 0.5, 0.6,
and 0.7.
5.2.3 Genetic Programming
To study the scaling laws of BNPMGs, we employ GPSR to obtain symbolic expressions that can
describe the relation between properties, such as Young’s modulus, ultimate tensile strength (UTS), and
other variables, which in our case, are porosity and composition. The symbolic expressions are represented
by binary expression trees.
193
There are two types of tree nodes. The first type is composed of operation
52

nodes, e.g., plus, minus, multiplication, division, and exponential. The second type is composed of terminal
nodes.

Figure 21. Binary tree representations of symbolic expressions used in GP. (a)-(c) are examples of binary
tree representations. (b) and (c) are two different binary trees of equivalent symbolic expressions. Circles
and hexagons represent binary and unary operations, respectively. Squares represent terminals, i.e.,
numbers or variables.
Terminal nodes are variables and parameters. Among the first type, i.e., operation nodes, there are
still two subtypes. The first one is composed of binary operation nodes, e.g., plus and minus. The second
type is composed of unary operation nodes, e.g., exponential and sine. The binary operation nodes always
have two child nodes, while the unary nodes only have a left child node. While it is not a requirement, in
this work, we chose to have unary operation nodes that are differentiable functions. This choice is made to
allow efficient optimization of coefficients using the gradient-based Levenberg–Marquardt algorithm. In
the literature, one can find implementations of the GPSR method also using nondifferentiable or
discontinuous functions, such as the absolute function, |𝑥|, step function, 𝜃(𝑥), or sign function, sign(𝑥).
191

53

However, such functions demand less efficient methods to obtain the optimized coefficients, e.g., by using
the GP method itself.
Nodes of different types and values are randomly chosen with assigned probabilities when creating
a random expression. Figure 21(a) illustrates an example of the binary tree representation of the symbolic
expression given in Eq. (62)
2
+
,
Óexp(𝑥+𝑦)+
0.5
sin(𝑥−𝑦)
Ô 62
Each generation contains 80 symbolic expressions as its population. The parameters, if there are
any, in each expression are fitted by the Levenberg–Marquardt algorithm. Expressions in each generation
are sorted according to a score calculated based on the mean square difference between predicted and true
values and a penalty value based on the number of nodes and coefficients. Evolutionary processes are then
applied to each generation of symbolic expressions. The top 15% expressions in each generation are
regarded as “elite” and are passed directly to the next generation. The remaining expressions in the next
generation are obtained from the current generation of expressions following two evolutionary steps:
crossover between two expressions and possible mutation of the offspring expressions obtained from the
crossover step. In the crossover step, two expression trees are chosen, and each one of them is randomly
cut at a certain node resulting in a branch. The branches are then attached to the other tree at the node where
it was cut. In the mutation step, an offspring tree is cut at a certain node, and a randomly generated
subexpression is attached to that location. Similar to a natural selection process, expressions in the
population with low scores have a higher chance of crossover to have offspring than those with high scores.
Each offspring expression has a 20% chance of mutating; otherwise, it is directly passed to the next
generation. Crossover and mutation steps occur iteratively with the current generation of expressions until
the complete set of 80 expressions of the next generation is obtained. After many generations of such
evolutionary processes, the top expressions are expected to be able to predict well the scaling law of the
mechanical behavior of BNPMGs.
54

Several practical concerns need to be addressed in the GP algorithm. The first one is the generation
of exceedingly long expressions, which may lead to overfitting. To alleviate this, we add a penalty to the
number of nodes, i.e., the more nodes contained in the expression, the larger the penalty. A restriction is
also applied to the depth of the expression trees, which is limited to seven. An additional concern is the
existence of equivalent expressions, i.e., expressions with different binary tree representations yet the same
value, e.g., Eq. (63) and Eq. (64). The binary tree representations of Eqs. (63) and (64) are shown in Figs.
21(b) and (c).
2𝑥(2+3𝑦) 63
4𝑥+6𝑥𝑦 64
The GP algorithm is implemented in Python 3.7.13 with the help of NumPy 1.21.6, Sympy 1.7.1,
and SciPy 1.4.1 libraries. To avoid equivalent expressions in the population, we take advantage of the
capabilities of the Sympy library. Each generated expression is checked for any equivalence in the
population. If an equivalent expression is found in the population, the newly generated expression is
eliminated.
5.3 Results
5.3.1 Geometrical and morphological features of bicontinuous NPMG
The initial NPMG structures for 5 different porosities, 𝜌 =0.3, 0.4, 0.5, 0.6 and 0.7, are shown in
Fig. 19. The structures are self-similar to each other due to the level set method, i.e., the structures are
generated from a common topology by scaling the porosity level. These bicontinuous structures consist of
connected solid and porous phases. Figure 19 also shows the normalized curvature distribution for the
different initial BNPMG structures. The two principal curvatures, 𝜅
(
and 𝜅
)
, are calculated by fitting
bivariate surfaces represented by the Eq. (65), considering the position of each surface atom and other
surface atoms in its neighborhood, using the least square method
𝑧(𝑥,𝑦) =𝑎𝑥
)
+𝑏𝑦
)
+𝑐𝑥𝑦+𝑑𝑥+𝑒𝑦+𝑓 65
55

where (𝑥,𝑦,𝑧) are Cartesian coordinates and 𝑎, 𝑏, 𝑐, 𝑑, 𝑒, and 𝑓are parameters to be determined. Once the
parameters are found, the two principal curvatures can be analytically calculated
99
using the Eqs. (66), (67),
(68), (59) and (60):
𝑓
+
=𝑑+2𝑎𝑥+𝑐𝑦,𝑓
,
=𝑒+2𝑏𝑦+𝑐𝑥,𝑓
-
=−1,𝑓
++
=2𝑎,𝑓
,,
=2𝑏,𝑓
--
=0,𝑓
+,
=𝑐,𝑓
+-
=0,𝑓
,-
=0
66

𝜅
U
=(−1 𝑜𝑟 1)×
𝑓
+
)
(𝑓
,,
+𝑓
--
1+𝑓
,
)
(𝑓
++
+𝑓
--
)+𝑓
-
)
(𝑓
++
+𝑓
,,
1
2∙(𝑓
+
)
+𝑓
,
)
+𝑓
-
)
)
H/)
−
𝑓
+
𝑓
,
𝑓
+,
+𝑓
+
𝑓
-
𝑓
+-
+𝑓
,
𝑓
-
𝑓
,-
(𝑓
+
)
+𝑓
,
)
+𝑓
-
)
)
H/)

67
𝜅
m
=2
𝑓
+
𝑓
,
(𝑓
+-
𝑓
,-
−𝑓
+,
𝑓
--
1+𝑓
+
𝑓
-
(𝑓
+,
𝑓
,-
−𝑓
+-
𝑓
,,
1+𝑓
,
𝑓
-
(𝑓
+,
𝑓
+-
−𝑓
,-
𝑓
--
1
(𝑓
+
)
+𝑓
,
)
+𝑓
-
)
)
)
+
𝑓
+
)
(𝑓
,,
𝑓
--
−𝑓
,-
)
1+𝑓
,
)
(𝑓
++
𝑓
--
−𝑓
+-
)
)+𝑓
-
)
(𝑓
++
𝑓
,,
−𝑓
+,
)
1
(𝑓
+
)
+𝑓
,
)
+𝑓
-
)
)
)

68
It should be noted that the sign of the mean curvature 𝜅
U
shown in Eq. (67) is 1 (−1) when the
calculated surface normal is the same (opposite) as the actual normal, which can be determined following
a procedure detailed in section 4.2.1. The curvatures are then normalized by the specific area. The curvature
distribution for the 𝜌 = 0.5 structure is symmetric with respect to the dashed lines in Fig. 19(c). Most of the
surface at this porosity level is composed of saddle points. The surface has equal amounts of concave and
convex areas indicated by the yellow and red rectangular regions, respectively. The remaining plots indicate
that by increasing (decreasing) the relative density by 0.1, the surface displays more concave (convex) areas.
As the relative density reaches 𝜌 = 0.7 (𝜌 = 0.3) the structure contains mostly concave (convex) surface.
An important characteristic of NPMG is the ligament size distribution. Figure 20 shows the
calculated ligament size distributions for the five initial NPMG structures. In the calculation of the ligament
sizes, we used 2 nm-thick slices randomly cut from the structures. These cross-sections are then used to
56

evaluate the ligament size, which is done using the AQUAMI software.
97
For relative densities from 𝜌 =
0.3 to 𝜌 = 0.7 the estimated ligament size increases from 3.6 to 6.3 nm.
5.3.2 Mechanical Properties of bicontinuous CuZr NPMG
We first focus on the mechanical properties of the samples with a Cu concentration of 0.64. The
engineering stress-strain curves for Cu 64Zr 36 NPMG during tensile loading simulation of seven samples
with varying relative densities from 0.3 to 0.9 are represented in Fig. 22(a). The curves indicate a steady
drop in both the stiffness and strength of the samples as the relative density decreases. For the samples with
relative densities in the range from 𝜌 = 0.3 to 𝜌 = 0.7, the stress decreases gradually after reaching the apex,
ultimate tensile strength (UTS), and the samples fail at around the strain value e = 0.4. In contrast, the
samples with relative densities 𝜌 = 0.8 and 𝜌 = 0.9 display a distinct behavior. The UTS is much higher (>
1.5 GPa) compared to the other five samples. In addition, the stress level drops sharply after reaching the
UTS, indicating a different deformation mechanism. The stress-strain curves then stabilize at a certain level.
This is a typical result of simulations of BMGs with PBC.
194

To investigate the effect of composition on the mechanical behavior of the CuZr NPMG samples,
we consider the systems with initial dimensions 53.79 nm ´ 53.79 ´ 53.80 nm, and vary the Cu
concentration of the samples with different porosities. For each combination of porosity and composition,
a tensile loading simulation is performed to calculate the stress-strain relationship. Figure 22(b) shows the
engineering stress-strain curves obtained for samples with relative density ρ = 0.5 and CuxZr 1-x composition
𝑥 = 0.72, 0.64, 0.50, 0.36, and 0.25. After performing the complete set of tensile loading simulations, we
evaluate the Young’s modulus (E) and the UTS for each case. For the systems with relative density 𝜌 = 0.5,
it can be observed from the elastic regions that the E values increase with Cu concentration from 0.25 to
0.64. The exception is the Cu 72Zr 28 sample, which has a lower E value compared to the Cu 64Zr 36 sample.
All samples display their UTS values at around e = 0.8. The UTS values also increase with Cu concentration
from 0.25 to 0.64 to then drop in value from concentration 0.64 to 0.72. All BNPMG samples with relative
density 𝜌 = 0.5 exhibit considerable plasticity compared to their BMG counterpart.
57

Figure 22. Tensile loading engineering stress-strain relationships for the BNPMG samples investigated. (a)
Stress-strain curves for Cu 64Zr 36 samples with different relative densities. (b) Stress-strain curves for
samples with relative density 𝜌 = 0.5 and different compositions.
Nevertheless, a clear change in plastic deformation mode and increase in ductility is observed in
Fig. 22(b) when Cu concentration drops from 0.64 to 0.25. In particular, the Cu 25Zr 75 system experiences a
long plastic deformation regime displaying no signs of failure in the strain range displayed in Fig. 22(b).
Nonetheless, the ductilities are similar for Cu 72Zr 28 and Cu 64Zr 36 samples. The complete set of E and UTS
values, calculated from the mechanical behavior of the samples with different relative density–composition
combinations are displayed in Fig. 23. It can be observed that for low relative density samples, the values
of E and UTS are similar for different compositions. For relative densities different than 𝜌 = 0.5, the plots
also indicate that the values of E and UTS increase with Cu concentration from 0.25 to 0.64. The results
indicate that the E and UTS values for Cu 64Zr 36 and Cu 72Zr 28 are similar at all relative densities considered.
58

Figure 23. Young’s modulus and ultimate tensile strength for samples with different compositions as a
function of relative density.
For a system with a given relative density and composition, we can expect an effect of the system
size on the mechanical behavior since the ligament size will scale with the size of the structure. To evaluate
the possible size effects, additional simulations are performed. We construct multiple Cu 64Zr 36 samples with
different system sizes to evaluate the dependence of the mechanical behavior on the ligament size while
preserving the porosity and other morphological features.
59

Figure 24. Young’s modulus and ultimate tensile strength for Cu 64Zr 36 samples with different system size
as a function of relative density.
We consider system sizes in the range from 43.03 nm × 43.03 nm × 43.04 nm to 80.69 nm × 80.69
nm × 80.71 nm. We perform a series of tensile loading simulations considering all combinations of system
size and relative density for the fixed Cu 64Zr 36 composition. Figure 24 displays the values of E and UTS
evaluated from the stress-strain curves obtained from the tensile loading simulations. It should be noted
that for systems with small dimensions, low relative density, i.e., 𝜌 = 0.3, cannot be achieved since the
ligaments in those systems are unstable and will collapse during the relaxation phase due to their
exceedingly small size, i.e., ligament sizes smaller than 2.8 nm. The results show that neither E nor UTS
have a clear dependence on system size. The values for these two physical quantities are close for all system
sizes and do not show any clear trend with increasing or decreasing system size while keeping the porosity
and other morphological features constant.
5.3.3 Relative Density Effect on Deformation and Failure Modes
To better understand the change in overall elastic and plastic deformation modes for the samples
with different relative densities, we need to quantify the local deformation developed in the structures
60

during the tensile loading process. Figure 25 shows the von Mises atomic local strain,
129
for the
representative samples with 𝜌 = 0.4, 0.7, and 0.9 at an overall engineering strain of 0.36. One can observe
a transition in deformation and failure modes from low to high relative density. The deformation of the
samples with ρ = 0.4 and 𝜌 = 0.7 is dominated by the deformation of ligaments in the structure. Those
ligaments initially deform elastically and after exceeding their tensile strength, develop necking and fail.
The structures fail when massive failure of ligaments is triggered throughout the nanoporous structure, as
reported previously.
195
In contrast, the deformation in the sample with 𝜌 = 0.9 occurs following a distinct
mechanism. Different than the BNP samples with low relative density, the system with 𝜌 = 0.9, which has
closed pores, deforms and fails in a similar way as a typical BMG, i.e., by generation and evolution of a
critical shear band. Nonetheless, different than in a BMG, where shear transformation zones are generated
homogeneously throughout the samples during tensile loading, in the sample with 𝜌 = 0.9 the plastic
deformation onset is dominated by STZs that are generated heterogeneously at the pores of the structure.
The STZs evolve, generating an array of incipient shear bands in the nanoporous structure, which delocalize
plastic deformation, resulting in smoother stress release from the UTS point compared to that displayed by
a BMG sample. Eventually, one of the incipient shear bands becomes critical and dominates the
deformation of the structure, as displayed in Fig. 25(c). The critical shear band traverses the structure at
nearly 45 degrees with the loading direction in a similar way as a dominant shear band in a BMG. The
change in deformation and failure mode observed from the analysis of the microstructure evolution is
consistent with the stress-strain curves shown in Fig. 22(a).
61

Figure 25. Deformation and failure of Cu 64Zr 36 BNPMG samples with different relative densities. (a) and
(b) samples fail by collective necking and failure of ligaments at 𝜌 =0.4 and 0.7. (c) the sample with 𝜌 =
0.9 fails by shear banding.
5.3.4 Scaling Laws Uncovered by Genetic Programming Assisted Symbolic Regression
It is instructive to understand how the mechanical properties of nanoporous materials scale with
their structural features. In this work, GPSR is applied to find scaling relationships between materials
properties: E and UTS (𝜎
7
) and BNP variables: relative density (𝜌) and Cu concentration (𝑐). Each value
of E and UTS corresponding to a combination of relative density and composition considered in this work
is used as input data for the GPSR. The initial population of 80 expressions describing the relationships is
generated randomly. After multiple GPSR runs, we obtain various expressions that are able to describe
accurately the scaling relationships between mechanical properties and system variables. We display four
such expressions developed for E and UTS in Table 3 and 4, respectively. The corresponding root mean
square error (RMSE) and normalized root mean square error (NRMSE) are calculated and displayed aside
to quantify the accuracy of the expressions. The RMSE is calculated as 𝑅𝑀𝑆𝐸 =
n∑ (+
.
M+
/
p)
& 0
.12
!
, where 𝑁 is
the number of data point, 𝑥
$
the actual observations and 𝑥
"
Ø the estimated observations. The NRMSE is
calculated as 𝑁𝑅𝑀𝑆𝐸 =
=qrs
./0 (+
.
)M.tu (+
.
)
. The variety of equations and their diverse structures show the
ability of the GP method to uncover scaling laws with considerably independent forms. The top expressions
in the two tables are considered the best scaling relations based on their accuracy and compactness. While
62

they may not provide the ultimate lowest NRMSE values, they have a relatively simpler functional form
and a lower number of coefficients.
Expression RMSE (GPa) NRMSE
𝐸 =(𝑐(2𝑐+2)+2.5)(𝜌
)
(𝑎
K
+𝜌)−1) 0.88 0.018
𝐸 =𝑐𝜌(𝑎
K
+𝜌)(𝑎
(
+𝑎
)
𝑐) 1.24 0.032
𝐸 =𝑎
K
+𝑎
(
5
+𝑎
)
𝑐+
𝑎
H
𝜌
+𝑎
T
𝑐
G(\
3
v\
4
5
((
6
7(
8
$7$)
:
)
+𝑎
S
𝜌
)

0.45 0.0093
𝐸 =(𝑎
K
𝜌+𝑐×cosh2𝑐)(𝑎
(
+𝑎
)
𝑐+𝑎
H
𝜌
\
;
) 0.97 0.020
Table 3 Scaling laws for Young’s modulus obtained by GP.
It should be noted that in the GPSR implementation used in this work, we allow the algorithm to
use coefficients, such as 𝑎
K
,𝑎
(,
etc, that are fitted, and also constants, which are multiples of 0.5, such as
0.5, 1, 1.5, 2, etc, which are kept fixed. That allows the generation of expressions with common physical
roots such as 𝑥
)
,
+
)
, and
(
w((M+)
. The comparison between the true values of E and UTS vs. those calculated
from the first expressions shown in Table 3 and Table 4, are displayed in Figs. 26(a) and 26(b), respectively.
The orange lines are used as a reference and indicate perfect matching. Besides the values predicted with
the equations developed with GPSR, Fig. 26 also shows values predicted with existing scaling laws:
86

𝐸
𝐸
x
=𝐶
(
(𝐶
)
𝜌+𝜌
)
) 69
𝜎
7
=𝜎
x
𝐶
H
𝜌 70
Expression RMSE (GPa) NRMSE
𝜎
7
=2𝜌(𝑎
K
+𝜌)×𝑒
5
0.12 0.49
𝜎
7
=(𝑎
K
𝜌+𝑐𝜌)
).Q
0.13 0.54
𝜎
7
=b𝑎
K
+𝑎
(
𝑐𝜌
\
&
G
d(𝑎
H
+𝜌)(𝑐+𝜌) 0.11 0.46
𝜎
7
=𝜌
)
(2𝑐+𝜌)(𝑎
K
𝜌
\
2
5v\
&
Gv\
<
5
$
+𝑎
T
𝜌+𝑐) 0.11 0.46
Table 4 Scaling law for ultimate tensile strength obtained by GP.
63

5.4 Discussion
5.4.1 System Size Effect on Mechanical Properties
As we change the system size, we found that the ligament diameter has negligible influence on the
mechanical properties of NPMG for all compositions when the porosity and other morphological features
are kept fixed.

Figure 26. Comparison between existing and GP-derived scaling laws accuracy for the Young’s modulus
and ultimate tensile strength. (a) Predicted vs. simulated values for the Young’s modulus and (b) the
ultimate tensile strength of samples with different relative densities and compositions. Perfect correlation
is indicated by the reference orange line.
Similar results have been found in the study of open-cell porous Cu 50Zr 50 MG structures with 𝜌 = 0.5 by
Zhang et al.
159
The results obtained from tensile loading simulations indicate the absence of dependence
between E, the UTS, and the yield strength with ligament size. Arguably, the absence of a clear effect of
ligament size is linked to the absence of a change in deformation mode for the MG ligaments, which should
occur at a much larger critical ligament size as reported for MG nanopillars.
140,196
These results for MG
64

samples contrast with the behavior predicted for crystalline nanoporous materials. Sun et al.
10
investigated
the effect of ligament size on the mechanical properties of BNP gold and found negligible changes in the
values of E and the yield strength. Nonetheless, they found that UTS values decrease with increasing
ligament size.
5.4.2 Composition Effect on Mechanical Properties
Changes in the composition typically lead to changes in the mechanical properties of materials.
Our results indicate that the values of E and UTS increase with Cu concentration in the concentration range
from 𝑥 = 0.25 to 0.64. However, there is no clear change in E and UTS values for concentrations beyond 𝑥
= 0.64. Similar results have been found in the study of Cu 50Zr 50 and Cu 64Zr 36 BMGs,
197
in which the Cu 50Zr 50
samples display lower values of E, yield stress, and UTS. In addition, higher Cu concentration significantly
decreases the observed ductility of bicontinuous NPMGs, as observed in Fig. 22(b). That is in agreement
with observed behavior in BMGs.
197
Previous investigations indicate that a higher fraction of atomic
icosahedra motifs in the amorphous structure is linked to the higher strength and lower plasticity displayed
by samples with higher Cu concentration.
198,199
Ward et al.
192
found that an increase of Cu fraction in CuZr
MG leads to an increase of the fraction of Cu centered <0,0,12,0> full icosahedra until it reaches a maximum
value at 0.7 Cu fraction, followed by a slight decrease with further increasing of Cu fraction. A higher
density of icosahedra in the structure may lead to a network of interconnected icosahedra motifs. The
network then behaves as a backbone that is resistant to stress-induced shear transformations during
deformation. Imran et al.
200
studied the mechanical properties of Cu 25Zr 75, Cu 50Zr 50, and Cu 75Zr 25 MGs
using indentation tests and concluded that samples with higher Cu concentrations display higher hardness
and resistance force at the same load depth. Lee et al.
201
suggested that the connectivity of the icosahedra
network reaches 98% at 0.65 Cu fraction. Our results agree with these results and suggest that such
interconnected icosahedra motif is achieved at 𝑥 = 0.64.
192,201

65

5.4.3 Existing Scaling Laws with Relative Density
It is interesting to compare the scaling laws obtained by the GP method with existing scaling laws.
After thoroughly studying the mechanical properties of open-cell foams with micron or larger characteristic
length scale, Gibson and Ashby proposed the well-known scaling law
202
for Young’s modulus based on
Timoshenko beam theory shown in Eq. (71),
𝐸
𝐸
x
=𝐶
T
𝜌
)
71
where 𝐶
T
is a coefficient to be determined. In the study of Cu 50Zr 50 BNPMGs, Zhang et al.
158
found the
relation between Young’s modulus and relative density to be in good agreement with Eq. (71). From our
results, after converting Eq. (71) to 𝐸 =𝐸
x
𝐶
T
𝜌
)
and finding the optimized value for 𝐶
T
, considering all
data, the obtained NRMSE value is about 2%.
While this accuracy is satisfactory, the scaling law proposed by Sun et al.
10
shown in Eq. (69)
describes better how the mechanical properties scale with the relative density of BNPMGs, i.e., it results in
a NRMSE value of 0.77%. The proposed scaling law shown in Eq. (69) combines the effect of porosity
202

as well as the effect of ligament deformation.
203
Sun et al.
10
used the cubic array unit cell model established
by Gibson and Ashby, which can account for the effect of ligaments.
202
The relative density 𝜌 was assumed
to be related to the ligament aspect ratio by 𝜌 ∝(𝑑/𝑙)
)
, where d and l are the ligament diameter and length,
respectively. In addition, the normalized Young’s modulus was also assumed to follow the same
relationship
s
s
=
∝(𝑑/𝑙)
)
.
10

Thus,
s
s
=
∝𝜌. The linear term was then added to Eq. (71) to obtain the proposed scaling law shown
in Eq. (69). We use our data values for Young’s modulus at different relative densities of MGs of different
Cu concentrations and Eq. (69) to find the optimized values for the coefficients 𝐶
(
and 𝐶
)
. The comparison
between the simulation results and values obtained from Eq. (69) is shown in Fig. 26(a). We can see that
this scaling law fits well the simulation data, displaying an NRMSE of 0.77%. The results show that the
66

relationship provided in Eq. (69) further improves the agreement with simulation results from that obtained
using Eq. (71).
Equation (69) was first proposed to describe the scaling relation for nanoporous gold.
10
Later it was
also applied to study nanoporous aluminum by Winter et al.,
204
who also found it to describe well the
experimental data. In the original work of Sun et al.,
10
Eq. (69) was tested considering relative densities in
the range from 0.24 to 0.36, which is significantly narrower than in our study. Our results show that Eq.
(69) is also able to accurately predict the Young’s modulus of CuZr NPMG considering a wide range of
relative densities and also compositions. The fact that the same scaling law applies for both nanoporous
crystalline metals and MGs reveals that their mechanical behavior is intrinsically similar, at least in the
elastic region.
There are alternative scaling laws proposed in the literature for the mechanical behavior of
nanoporous materials.
24,103,205
Pia et al.
205
studied the effective Young’s modulus of nanoporous gold
obtained from selective dissolution of Ag from an Ag 70Au 30 alloy, which was followed by thermal annealing.
The ligament size of the nanoporous gold samples varied from 5 to 200 nm. They simplified the geometry
of the nanoporous structure to that of a representative cubic unit cell and derived a scaling law of effective
Young’s modulus as a function of the shape factors of the ligaments, as well as the Young’s modulus and
Poisson's ratio of bulk gold. Mangipudi et al.
103
investigated the stiffness of nanoporous gold considering
three independent topologies: (i) 3D reconstruction from tomography of experimental nanoporous gold
samples (ii) spinodal decomposed structure obtained from phase-field simulations and (iii) a gyroid
structure. They claimed that while the scaling law of Young’s modulus for all topologies can be well
described by Eq. (73), the pre-factor 𝐶
T
is linearly dependent on the genus of the particular nanoporous
structure. Badwe et al.
24
studied the tensile properties of dealloyed nanoporous gold. The samples were
thermally annealed to coarsen the ligaments. The Young's modulus was found to obey a power law.
However, the calculated exponent was larger than that predicted by the Gibson-Ashby scaling law shown
in Eq. (71). They further generalized the scaling law of Young’s modulus proposing the law shown in Eq.
(72):
67

𝐸
𝐸
x
=𝐶
(
(𝐶
)
𝜌+𝜌
y
)+𝐶
H
72
Among the scaling relationships proposed in the literature, the modified Gibson and Ashby law provided
in Eq. (69) is the most effective and succinct at describing relationsips in a large range of relative densities.
As for the UTS, Sun et al.
10
proposed a linear relationship between the UTS and the relative density
as shown in Eq. (70). It was demonstrated that Eq. (70) can describe well the scaling relation between
porosity and UTS for BNP gold.
10
However, as illustrated in Fig. 8(b), the actual values of the UTS of CuZr
bicontinuous NPMGs deviate considerably from the predictions of Eq. (70) with the value of the coefficient
𝐶
H
optimized to the data. The corresponding NRMSE is found to be 15%. A similar significant deviation
in the predictions of UTS was also reported for Cu 50Zr 50 bicontinuous NPMGs.
158
This contrast in results
for crystalline and amorphous samples implies that the scaling law for UTS of crystalline and amorphous
materials is intrinsically distinct, which is expected considering the contrasting deformation mechanisms
displayed by crystalline and amorphous metals.
5.4.4 Comparison of GPSR Uncovered and Existing Scaling Laws
The GPSR method is able to generate functional expressions with different forms that can
accurately describe the scaling relationships between mechanical properties and system variables as shown
in Tables 3 and 4. For simplicity, only the top two expressions in Tables 3 and 4 are used in the calculation
of the predictions shown in Figs. 26(a) and 26(b), respectively, as GP scaling laws for E and UTS. Those
two expressions are able to describe well the values of E and UTS of BNP CuZr MGs, with NRMSE 1.8%
and 5% for E and UTS, respectively. Compared to the predictions made with Eq. (70), which is the current
scaling law for nanoporous gold, the predictions of the GP scaling law for UTS display a much better
agreement with the simulation data. The corresponding NRMSE for the GP scaling law at 5% is much
lower than that of the current scaling law at 15%. The GP method is able to uncover the nonlinear
relationship between UTS and relative density and assign appropriate expressions to it.
As for the predictions of the Young’s modulus, the results indicate that both the current scaling law
and the GP-derived relationship predictions are in good agreement with the simulation data, with
68

corresponding NRMSE of 0.77% and 1.8%, respectively. Even though the agreement is good in both cases
the GP scaling law is slightly less accurate than the predictions made with Eq. (69). However, it should be
noted that the scaling law derived with the GPSR method is universal and applicable to all compositions.
In contrast, the current scaling law requires knowledge of the 𝐸
x
value, which is used as an additional
parameter to define the specific scaling law at each composition. Instead, the GP method is capable of
uncovering the dependence of 𝐸
x
on composition, which is typically not a simple increasing or decreasing
relationship, and incorporating it in the universal scaling law predicted. For instance, both E and UTS of
bulk CuxZr 1-x MG increase from x = 0.25 to 0.64, and keep the same value in the range 0.64 < x < 0.72.
While the accuracy of the predictions using the GP-derived universal scaling law is good, it can be further
improved if one allows more complicated expressions with more fitting coefficients to describe the
nonlinear relationships, e.g., the third expression in Table 1. An important aspect of the GP method is that
it is also able to retrieve the scaling law proposed based on physical insights, e.g., the second expression in
Table 1 has the same form as the current scaling law shown in Eq. (69).
5.5 Conclusions
In this work, we have investigated the mechanical properties of BNP CuZr MGs with different
relative densities, system sizes, and compositions during tensile loading using MD simulations. With
samples with self-similar structures, the relative density was varied to tune the surface morphology and
ligament size. We found that the system size of NPMG, in the range considered, has little impact on the
mechanical properties, while the relative density and the Cu concentration cause strong effects. A brittle-
to-ductile transition occurs at relative density 0.7. NPMGs with large relative density, e.g., 0.9, fail by
developing a critical shear band, which is nucleated at large pores, while those with small relative density
fail by progressive necking of ligaments. We employed the GP method to uncover universal scaling laws
and compared their predictions with that of existing scaling laws. The scaling law based on Gibson and
Ashby's theory considering the effect of ligaments, can be applied well to the Young’s modulus of BNP
CuZr MGs. However, the existing linear scaling law for UTS of nanoporous structure is not able to offer
69

accurate predictions for the NPMGs studied in this work. Relationships uncovered by GPSP are able to find
succinct expressions to describe accurately the dependence of mechanical properties on relative density and
Cu concentration with physical insights.

70

Chapter 6
Pore-Size Dependence of Permeability in Bicontinuous Nanoporous
Media
6.1. Background
Fluid flow through porous media has been drawing experimental and theoretical attention for
several decades.
47,49,206–209
The understanding of the different aspects intrinsic to fluid flow in confined
space has progressed from the application of the kinetic theory to fluid transport through pores
210–212
and
analysis of experimental results
47,212
in the light of simulation predictions.
208,213–215
A major challenge in the
study of flow in nanoporous media is the description of its complex geometry, which involves many factors,
such as porosity, pore size, pore shape, and pore distribution. Among these, pore size is a crucial factor. It
has been shown that the fluid flow can behave in distinct ways in fine-grained nanopores and coarse-grained
pores due to the increasing effect of fluid-wall interactions in nanopores.
51

Simulation techniques are particularly promising for the understanding of flow behavior in
nanoporous media. For example, Ghadiri et al.
44
applied computational fluid dynamics to simulate the water
desalination process using nanoporous membranes, reporting experimentally validated results. Razavi et
al.
216
studied carbon dioxide mass transport from a gas mixture through nanoporous membranes equipped
with an absorbent solution and activators of absorption. It was found that an increase in flow rate and
concentration of activators enhanced CO 2 absorption. Pan et al.
217
simulated two-phases flow in a complex
pore network using the Lattice Boltzmann Method (LBM). Zhao et al.
218
simulated liquid flow through a
single nanopore and nanoporous media also with LBM. Atomistic simulations were instrumental in many
studies. For instance, Takaba et al.
219
and Sokhan et al.
47
used MD simulations to describe water flow
properties through a nanopore. Adiga et al.
220
studied fluid gating through polymer-grafted nanofluidic
channels. Strong and Eaves
212
used Gaussian dynamics to study the dynamics of water in porous two-
71

dimensional crystals. Wang et al.
221
studied the pressure-driven flow of water in CNT nanoporous
membranes.
Modeling methods can be very useful, yet they have their limitations. While continuum scale
simulations might not be able to capture the physical and chemical behavior of fluid particles in nanoporous
media,
50,51
atomistic simulations are often limited by their small length and time scales. As a compromise,
the DPD method, a mesoscale model originally introduced by Hoogerbrugge and Koelman
222
can bridge
the length and time scales gap between continuum and atomistic simulations. The DPD formulation models
the dynamics of atom clusters instead of single atoms, while preserving important atomistic-like
interactions. The particle-particle interaction is described by a soft potential that allows overlapping
between DPD particles. In the DPD method, interactions include conservative, stochastic, and dissipative
components. DPD has been previously employed in the study of colloid suspension flow
223–226
and flow in
nanoporous media.
222,227,228
As a variant of the original DPD model, many-body DPD (mDPD) features an
attractive term and a density-dependent repulsive force term.
229
mDPD has a flexible parameterization that
allows the modeling of different fluids and is particularly suitable for the modeling of multi-fluid systems.
mDPD has been successfully applied in the study of liquid-vapor interface and surface tension,
214
single-
phase unsaturated flows in porous media,
230
interfacial tension between two fluids,
51
and fluids in
capillaries.
51
Recently, a graphics processing unit (GPU) accelerated mDPD package has been developed
by Xia et al.
231
enabling significantly better performance than its central processing unit (CPU) counterpart.
While the capabilities of fluid flow simulation tools, such as mDPD, have significantly evolved, long-
standing challenges remain, such as the understanding of the effect of nanoscale porosity in the fluid flow
in complex nanoporous structures.
An important aspect of the mDPD simulations is the interaction between fluid particles and the
nanoporous media walls. Because of the soft nature of the mDPD pairwise interaction, there is a chance
that fluid particles can penetrate the solid phase, which is usually also represented by DPD particles. This
is an unrealistic phenomenon and should not occur during the simulations. Several efforts have been put
72

into designing proper model boundary conditions to avoid such a problem.
232–235
A particularly simple
boundary condition was proposed by Xia et al.
51
using a denser configuration for the solid walls. With the
use of such a method, one can prevent the fluid from diffusing into the solid phase while producing a non-
slip boundary condition at the fluid-solid interface.
The topic of pore size effect on flow through nanoporous structures has been investigated in several
previous atomistic simulation studies. Sun et al.
236
studied the mechanisms of molecular gas permeation
through nanoporous graphene membranes while Hossain and Kim
237
studied liquid transport through
similar membranes. Yuan et al.
238
studied gas permeation through graphene nanopores. Glavatskiy and
Bhatia
239
investigated the dependence of the thermodynamic resistance on the radius of a nanopore and the
thickness of a pore wall, considering the transport of carbon dioxide and methane through carbon nanotubes.
Chen et al.
240
investigated the transport behavior of water molecules through carbon nanotubes. Liu and
Li
241
studied the interactions of the fundamental parameters, including binding energies, temperature, and
driving force, and their effects on the flow motion in nanoscale Poiseuille flows. Takaba et al.
219

investigated the transport properties of pressure-driven fluid flow in thin nanoporous membranes. They
found that the fluid flux and velocity distribution in the membrane pores were directly dependent on the
pore size and the interaction between fluid particles and pore walls. A linear correlation was established
between the density of the fluid in the membrane pores and the deviation of the flux estimated from the
Hagen-Poiseuille flow. Even though these previous studies of pore size effect are instructive and insightful,
they are all based on relatively simple geometries, e.g., a 2-dimensional circle or a cylindrical tube.
In this section, we focus on the effect of pore size on the fluid flow in realistic nanoporous structures.
We employ mDPD simulations to describe accurately fluid-pore interactions and BNP structures to
represent nanoporous rock structures.
126
The BNP models provide a stochastic description of the
morphology and pore size distribution while allowing for a direct investigation of the effects of average
pore size. The results show the dependency of permeability on pore size and provide physical insights into
such phenomena. The permeability describes how well a material allows liquids or gases to pass through it.
73

When the flux increases linearly with the applied pressure difference, the permeability of a material is
defined by Darcy’s law:
𝑞 =−
𝑘
𝜇
∇𝑝 73
where 𝑞 is the flux, ∇𝑝 is the pressure difference, 𝜇 is the dynamic viscosity of the fluid and k is the
permeability.
6.2 Method
6.2.1 Bicontinuous Nanoporous Model and its Tortuosity
To study the permeation process through nanoporous media, a system is prepared by filling a BNP
structure with fluid particles. Illustrations of the unit BNP structure as well as of the combined systems are
shown in Fig. 27. We then quantify the pore size distribution and tortuosity of the BNP system. In Fig. 27(c)
we show the pore size distribution for the S80 sample, calculated with the AQUAMI software.
97
The
tortuosity is quantified by applying the Dijkstra shortest path algorithm.
242

To calculate the tortuosity, we considered the fluid mDPD particles occupying the porous phase of
the model. The set of mDPD particles and the connectivity of every particle to their nearest neighbor
particles is used to create a graph, which is employed in the search algorithm. The particles with position
𝑥 < 0.01 𝐿, where L is the length from source to sink, are considered departure positions. Then for every
other particle, we applied the Dijkstra shortest path algorithm
242
to find the shortest path distance to all
departure positions. The tortuosity is then defined as τ=
2
,>
2
, where 𝐿
Zz
is the smallest among all the
shortest paths to the departure particles. The resulting tortuosity for all particles is shown in Fig. 27(b). Two
fluid reservoirs are placed adjacent to the nanoporous media, where the fluid is confined by two mobile
pistons. By the action of the two pistons, the fluid is made to move between the reservoirs by percolation
through the nanoporous structure.
74

Figure 27. DPD simulation models. (a) Unit nanoporous media morphology used in the construction of the
DPD models. (b)-(e) Four models with increasing square cross-sections 20, 40, 60, and 80 𝜎 wide. Blue,
gray, and red-colored particles represent DPD particles in the fluid, nanoporous media, and pistons.
6.2.2 Configurations of Solid and Liquid Phases
The solid phase present in the pistons and the nanoporous media is generated from a parent bulk
system composed of DPD particles at a density of 7.5
[
{
<
in an initially FCC crystal structure. This bulk
75

system was relaxed for 10,000 timesteps (Dt) (691 ps), where Dt = 0.01 𝜏 = 0.0691 ps, allowing the system
to reach thermodynamic equilibrium. Equilibrium is reached when all intensive quantities, including
temperature, chemical potential, and pressure reach a steady-state value in all regions of the system. In
order to focus on the fluid properties, the DPD particles in the nanoporous structure are made immobile,
while the DPD particles in the pistons move as a rigid body. The thickness of the pistons is 10 𝜎 (8.52 nm).
The side of the system square cross-section varies from 20 𝜎 to 80 𝜎 (17.04 nm to 68.16 nm) across the
different systems. For the nanoporous media, we start from a cubic solid bulk and then apply the method
proposed in Chapter 3 to generate a BNP structure template.
126
The wavenumber employed in the method
is chosen to be |𝒒| =
)O×T
2
, where 𝐿 is the edge of the cubic system. The initial Fibonacci grid contains
3000 waves. The level set threshold defining the porosity level is set as 𝑓(𝒓) < 0.385. The resultant
porosity of the bicontinuous structure is 𝜌 =0.65. The samples we prepared are morphologically consistent
to nanoporous crystalline metals,
84
non-crystalline metals,
40
polymers,
243
and silica
244
obtained from
selective dealloying or corrosion of a spinodally decomposed solution. Our bicontinuous samples have
similar features to that of natural nanoporous rocks, as they both feature stochastically distributed nanopores
and tortuous flow paths.
245–247
Nonetheless, nanoporous rock samples are known to display a complex
multiscale heterogeneous pore size distribution. Our well-defined bicontinuous model is useful to
characterize the effect of average pore sizes in the flow properties, ruling out other geometrical features.
Other simplified models, such as single cylindrical pores,
248
or systems with circular and square-shaped
pore channels
50
are also widely used to understand the flow properties in nanoporous systems. To rule out
the possible dependence of the results on the specific nanoporous morphologies used we use the same
nanoporous morphology template in all systems. We generate systems with 𝐿 equal to 20 𝜎, 40 𝜎, 60 𝜎,
and 80 𝜎 (17.04 nm, 34.08 nm, 51.12 nm, and 68.16 nm), and the corresponding average pore sizes are 4
𝜎, 8 𝜎, 12 𝜎, and 16 𝜎 (3.41 nm, 6.82 nm, 10.22 nm, and 13.63 nm). The initial length of the nanoporous
systems is fixed at 80 𝜎 (68.16 nm) allowing direct comparison of the flow properties among the different
nanoporous models. Thus, the systems with 𝐿 = 20 𝜎, 40 𝜎, and 80 𝜎 (17.04 nm, 34.08 nm,68.16 nm) have
76

4, 2, and 1 concatenating cubic unit nanoporous media, respectively. For the case of 𝐿 = 60 𝜎 (51.12 nm),
the system contains 1
(
H
cubic unit nanoporous media. With the solid phases assembled, the fluid particles
are added to the systems. The solid-liquid and liquid-liquid interactions are both described by mDPD
interactions. The parameters used for the interactions are 𝛾 = 4.5,𝐵 =25,𝑟
F
= 0.75, and 𝑟
5
=1.0. In our
simulations, 𝑟
F
is chosen to be 0.75 𝑟
5
. The choice of such 𝑟
F
and 𝑟
5
combination is made based on the fact
that mDPD simulations using a single cutoff value are unable to describe a stable liquid-vacuum
interface.
249
The use of a repulsive interaction with a cutoff 𝑟
F
<𝑟
5
can address this problem.
229
Several
values of 𝑟
F
, including 𝑟
F
= 0.5, 0.65, 0.75, 0.86 𝑟
5
are proven to be able to retrieve the empirical equation
of state of a vapor-liquid-coexisting system
𝑝 =𝜌𝑘
E
𝑇+𝛼𝐴𝜌
)
+2𝛼𝐵𝑟
F
T
𝜌
)
(𝜌−𝑐) 74
where 𝛼 and 𝑐 are constants.
229
The different choices of 𝑟
F
can be compensated from the correction constant
c. The choice of 𝑟
F
= 0.75 𝑟
5
was also made in several previous investigations, providing a reasonable
equation of state.
229,231,250,251
We use 𝐴 = −40 and −35 for the liquid-liquid and solid-liquid interactions,
respectively. The choice of the parameters follows the work of Xia et al.
51
on multiphase flow in shale rocks.
The choice of parameters provides good wettability and reasonable surface tension. Solid-solid interactions
are turned off since pistons are treated as rigid bodies and the nanoporous media is undeformable and
immobile. The initial density of the fluid is 6.544
[
{
<
, before the application of pressure. It should be noted
that the solid phase density, 7.5
[
{
<
, is slightly higher than the fluid density. This is set to prevent the
penetration of fluid particles into the solid phase. The force on the piston is calculated as:
𝑭=II𝒇(𝑟
$%
)
[
$;(
y
%;(
+I𝒈
%
y
%;(
75
where 𝒇(𝑟
$%
) stands for the force from a fluid particle 𝑖 to the piston particle 𝑗, and 𝒈
%
stands for the body
force added on piston particle 𝑗 when external pressure is applied to the piston. While the fluid is under
77

hydrostatic pressure the pistons receive only a pressure component along the flow direction, as shown in
Figs. 27(d)-(g).
DPD value Real value
𝜎 8.52 Å
𝑚 9.42×10
M)*
kg
𝜏 6.91 ps
𝜀 1.43×10
M)(
J
𝜌 =6.54
𝑚
𝜎
H
997 kg∙m
MH

𝑝 =7.65
𝜀
𝜎
H
17.67 MPa
𝜇 =6.25
𝑚
𝜎𝜏
0.1 mPa∙s
Table 5 Conversion between DPD values and real values.
After the composite systems are assembled, an initial relaxation step is performed where the left
piston is fixed, please see Fig. 28, while zero pressure is applied to the right piston allowing it to move
freely and reach an equilibrium position. The relaxation is done over 80,000 Dt, which is 800 𝜏. To apply
pressure difference and drive the fluid advective flow, we added a −
K.(}
{
force magnitude on each particle
in the left piston and 0.3, 0.4, 0.5, and 0.6
}
{
on each particle in the right piston, corresponding to a pressure
𝑃
(
= 7.65
}
{
<
and 𝑃
)
=22.97, 30.62, 38.29, 45.94
}
{
<
. The equivalent pressure difference generated in the
fluid reservoirs is then 15.31, 22.97, 30.62 and 38.29
}
{
<
. We use LAMMPS
127
to run all simulations of fluid
flow and OVITO
89
to perform visualizations. As a reference, we provide in Table 5 conversion data between
Lennard Jones and real values.
78

6.3 Results
6.3.1 Self-Closed Void During Relaxation
During the setup process, the liquid particles are initially introduced in an FCC arrangement. At
the first stages of the relaxation, the liquid particles move toward thermodynamic equilibrium in the
presence of the nanoporous solid phase with complex geometries. Due to the immediate wetting process,
several voids then appear in the nanoporous portion, as shown in Fig. 28(a). As relaxation progresses, the
liquid particles spontaneously fill up the voids due to the surface tension, which is well described by the
attractive term in the mDPD model. After the fluid voids are filled up and the fluid density is homogeneous,
the right piston moves to a stable position. The final relaxed distance between the two pistons for different
samples varies from 215 𝜎 to 225 𝜎 (183.18 nm to 191.70 nm).

Figure 28. Illustrations of the S40 model during the initial relaxation process. (a)-(c) Configurations at
different times indicating the nanovoid fill-up process. The colors follow the scheme used in Figure 27.
It should be noted that the fluid is unable to fill up the voids automatically when 𝐿 = 20 𝜎 (17.04 nm). This
is due to the large fraction of fluid particles that interact with the nanoporous walls in the nanoporous
structure, which constrains to some degree the flow of fluid particles. As the pore size, in this case, is
relatively very small, the fluid particles are unable to rearrange themselves enough to fill up the voids. Thus,
for this case, we add an extra force of 0.2
}
{
to the right piston during the first 400 𝜏 of relaxation. After that,
the simulation of fluid flow proceeds for additional 800 𝜏 without any extra force.
79

6.3.2 Pressure Distribution and Variation in the System
To induce fluid flow through the nanoporous models, various pressure levels are applied to the
right piston, while a constant pressure level is applied to the left piston. Due to the non-equilibrium nature
of the simulation, pressure oscillation is observed during the initial steps of the simulation. To compare the
pressure level in different regions of the sample, we decompose the system into 3 domains. A nanoporous
domain, a source domain from the rightmost part of the nanoporous to the leftmost part of the right piston,
and a sink domain from the rightmost part of the left piston to the leftmost part of the nanoporous media.
To illustrate the pressure evolution in the system, we show in Fig. 29 the pressure as a function of
simulation time for all samples when the pressure difference applied is 38.29
}
{
<
. It is observed that the
highest-pressure peak in the source is around 65 - 75
}
{
<
regardless of the pore size and sample size. The
pressure oscillation gradually is dissipated with the simulation time as the flow reaches a steady state. The
oscillation of the pressure in the sink region displays a much lower amplitude than in the source. For the
small pore size sample, S20, the oscillation in the sink pressure, even in the first stages, is negligible. The
steady average pressure for source and sink, after the decay of the oscillation, is the same for all the samples
with different pore sizes. In addition, the equilibrated pressure levels are equal to the externally applied
pressure to the right and left pistons, respectively.
It is important to understand how the local fluid pressure changes during the flow through the
nanoporous structure. In order to evaluate the local pressure distribution, we first calculate the hydrostatic
pressure on each atom, using the virial stress, we then average the value within a radius of 3 𝜎. The resulting
local pressure distribution is shown in Fig. 30. For simplicity, we only display the case with the pressure
difference of 38.29
}
{
<
and all the snapshots are taken after the systems achieve a steady-state flow. One can
observe that the local pressure in the source and sink regions have well-defined values besides
thermodynamic fluctuations. In contrast, the local pressure in the nanoporous media drops continuously
from the source to the sink reservoirs. To better quantify the pressure change in the nanoporous media, we
80

split the system into 10 bins and calculate the average atomic hydrostatic pressure in each bin. The
calculated pressure distribution is shown in Fig. 30(e), which indicates that the local pressure drops linearly
in the nanoporous media, regardless of pore and sample size.

Figure 29. Pressure evolution in the source and sink regions for different samples during the fluid flow
process.
6.3.3 Motion of Pistons
The motion of the pistons during the fluid flow process is monitored by recording the position of
each piston as a function of simulation time. Due to the different pore sizes and different permeation, the
difference in pressure applied to the two pistons generates distinct evolution of positions of the two pistons,
as shown in Fig. 31(a). At the beginning of the simulation, the position of the right piston significantly
oscillates while overall moving towards the nanoporous structure. In contrast, the position of the left piston
oscillates to a much lower amplitude.
81

Figure 30. Local pressure distribution in the systems during the fluid flow process. (a)-(d) Pressure
distribution in the S20, S40, S60, and S80 models when about half the sink fluid has flowed through the
nanoporous structure. (e) The plot of the pressure along the model length.
For both pistons, the oscillation is considerably dissipated by 250 t when the motion becomes
smooth. One can note that the moving speed of the pistons is significantly reduced on reduction of the pore
size, implying a strong pore size dependence on the permeability. The distance between the right and left
pistons is plotted in Fig. 31(b). This is due to the different pressure in the sink and the source and the
resulting different local densities in these regions. As more fluid particles flow from the source into the
sink, the region with lower density in the sink increases in detriment to the region with relatively higher
density in the source. The result is a gradual increase in the separation between pistons.
82

Figure 31. Position of the pistons during the simulation. (a) Position of the front and back piston and (b)
distance between the two pistons during the simulation time.
Consistent with the results shown in Fig. 30, the asymptotic rate of separation of the two pistons suggests
a decrease in permeability as porous sizes are reduced. For instance, the orange curve in Fig. 31 for the S20
model shows an almost negligible increase in the distance between the pistons within the time frame
considered.
6.3.4 Local Velocities
The local velocity distribution along the horizontal (x-) direction in the different models is shown
in Fig. 32. At the position of each particle, the local velocity is obtained by averaging the 𝑥 component of
velocity over the particle and its neighbors within a 2𝜎 cutoff.
83

Figure 32. Fluid velocity distribution in the different models. (a)-(d) Profile of local fluid velocity through
the different models. (e) Fluid velocity probability densities for the four models.
The illustrations are based on configurations taken at different times for each system when a steady-
state flow is achieved. Due to the large difference in particle velocity among the different systems, different
color scale bars are used. One can observe the inhomogeneous flow throughout the nanoporous media, with
some preferred channels displaying much higher velocity than others in large-pore-size samples. While the
velocity contrast is quite strong for S60 and S80, the contrast for S20 and S40 is weak. To understand this,
we collect all the local particle velocities and normalize them to obtain the velocity probability density
shown in Fig. 32(e). The velocity value at the probability peak is directly related to the permeability of the
84

nanostructure. The results show that the model S80 shows predominantly pressure-driven motion with an
average local velocity peaked at ~0.25
{
~
. In contrast, in the S20 model, the local velocity probability is
peaked at 0.02
{
~
, which is much smaller than that in large samples, namely S60 and S80. As the overall
flow rate is smaller than the magnitude of the particles random motion, the flow displays a homogenous
profile in the nanoporous media. The skewing seen from S20 to S80 is caused by the increasingly
inhomogeneous flow in the nanoporous structures at the large average pore.
The diffusion of the fluid particles under a strong gradient of pressure is further quantified by
calculating the mean square displacement (MSD). The time evolution of the MSD of fluid particles in the
nanoporous and sink regions of the samples with different pore sizes at 30.62
}
{
<
pressure difference is
shown in Fig. 33. The MSD accounts for the thermal diffusion of the particles and the drift along 𝑥-direction
caused by the different pressures applied to the pistons. The results indicate that the diffusion in both the
sink and the nanoporous structure increases with the pore size. This is a direct consequence of the larger
flow rate in samples with larger pore sizes. It should be noted that the larger MSD in the nanoporous region
of the samples with larger pore size is not only caused by the higher drifting rate along the 𝑥-direction since
the displacement components along the other two directions are also stronger. Indeed, the particles flowing
along the 𝑥-direction will be forced to change direction to follow the complex topology of the nanoporous
structure. The second observation is that the diffusion in the nanoporous region is stronger than that in the
sink region for all the samples except for S20. For S40, S60, and S80 samples, many particles move faster
because of the limiting number of channels available for them to go through in the nanoporous region.
However, in the S20 sample, the nanoporous topology is the finest, and the pore size is the smallest of all
the samples. A consequence of that is that fluid particles may follow a path mostly aligned with the x flow
direction. In addition, a large fraction of fluid particles is pinned by the surface effect of the nanoporous
solid phase. As a result, the MSD in the nanoporous region is smaller than that in the sink region for the
S20 sample. It should be noted that the mDPD model employed in this work is not able to quantify
accurately both diffusivity and viscosity of the fluid at the same time. However, the qualitative observations
85

demonstrate the impact of average pore size on fluid diffusive behavior. A modified mDPD model recently
proposed by Rao et al.,
251
is able to quantitatively describe diffusivity and viscosity concurrently, which is
appropriate for future in-depth work on fluid diffusion.

Figure 33. MSD of fluid particles in the nanoporous and sink regions for all samples at 30.62
}
{
<
pressure
difference.
Finally, in order to quantify the permeability in the different models, we evaluate the rate of
particles entering the sink region after 300 𝜏 after a steady-state flow is achieved. The specific time, 300 𝜏,
is chosen because the oscillation of the position of the left piston is negligible thereafter, as shown in Fig.
31(a), and the increase rate of the number of sink fluid particles follows a well-defined linear relationship.
It is observed that the density of particles at a pressure of 15.31
}
{
<
, which is the pressure in the sink region,
is 6.88
[
{
<
. Thus, the rate of particles entering the sink can be calculated as 6.88×𝑣
]z
×𝐿
)
, where 𝑣
]z
is the
velocity of the left piston and 𝐿 is the side of the cross-section of the system. The calculated fluxes per unit
area for all pressure differences and pore sizes considered are shown in Fig. 34. We find that the fluxes, for
all samples, increase linearly with the applied pressure difference in agreement with Darcy’s law Eq. (73).
To calculate the dynamic viscosity of the fluid, we employed the Müller-Plathe method.
252
We first
create a simulation box with dimensions 𝐿
+
×𝐿
,
×𝐿
-
= 80 × 40 × 40 𝜎
H
(68.16 × 34.08 × 34.08 nm
H
).
The box is then filled with the same mDPD particles as used in the permeation simulations. The system is
divided into 20 bins, from top to bottom, along the 𝑧-direction. At every 10 𝜏, momentum is changed in the
first bin and the middle (11
th
) bin, i.e., 20 particles with the most positive velocity along the x-direction in
86

the first bin will interchange their velocities with the 20 particles with the most negative velocity along the
x-direction in the 11
th
bin. With time evolution, the system forms a steady-state shear velocity profile. The
dynamic viscosity is then calculated:
𝜇 =−
𝑃
+
2𝑡𝐴
b
𝜕𝑣
+
𝜕𝑧
d
M(

76
where 𝑃
+
is the total momentum exchange along 𝑥-direction, 𝑡 is the total simulation time, and 𝐴 is the
cross-section area. The first factor,

!

, can be considered the momentum exchange rate per unit area, which
is equivalent to the shear stress. The ratio between the shear stress to the velocity slope defines the dynamic
viscosity. The factor 2 in the denominator comes from the periodicity of the system, such that the flux is
going in two directions, causing the effective momentum exchange to be reduced by a factor of 2. It should
be noted that the dynamic viscosity is directly related to the density and indirectly related to the pressure of
the system. In our non-equilibrium simulations, the fluid mDPD particles will experience pressure to
different levels causing corresponding density variations. The results indicate that the calculated value of
the dynamic viscosity increases from 4.55
[
{~
to 6.25
[
{~
when the system density increases from 6.544
[
{
<

to 6.977
[
{
<
. The average density in the nanoporous media varies from 6.88 to 7.02
[
{
<
for the applied
pressure difference in the range of 15.31 to 38.29
}
{
<
. For simplicity, we approximate the dynamic viscosity
for all the scenarios as 6.25
[
{~
.
We finally calculate the permeability for each sample according to Eq. (73). The relation between
permeability and pore size is shown in Fig. 35. An extra point is added at (0, 0), accounting for the limit
when pore size reaches zero and no fluid particle is able to flow from the source to the sink because of the
solid obstacle. The simulation data lies perfectly on a parabolic curve. In this section, we analyze the results
obtained from the novel mDPD non-equilibrium simulations of fluid flow through complex stochastic BNP
media with different average pore sizes, highlighting both non-steady and steady states. More discussion
about the dependence of permeability on pore size can be found in the following section.
87

Figure 34. Particle flux per unit area under different pressure differences for different samples.
6.4 Discussion
6.4.1 Pressure Evolution
It is instructive to consider the pressure evolution in the source and sink regions for the different
models shown in Fig. 29. As soon as the external pressure is applied to the pistons at the beginning of the
simulations, the pressure level in the source and the sink swiftly increases and oscillates with time until a
steady-state pressure level is achieved. As can be observed in Fig. 29 the amplitude of the pressure
oscillations in the source and sink are different due to the different external pressure applied to the right
and left pistons. As time evolves, the pressure oscillations in both source and sink are attenuated by
dissipative effects related to the viscous flow, represented in the simulations by the mDPD dissipative term.
We will focus on the pressure oscillations in the source reservoir since they are stronger and better defined
than those in the sink. One can note that in the S20 model, the dissipation of the pressure oscillation is
relatively slower than in the other three models, which have the same nanoporous topology and porosity
level. To understand that, we should consider the interaction of fluid particles with the solid particles of the
nanoporous media. When fluid particles are in contact with the nanoporous media, their interaction acts to
reduce the velocity of the fluid particles, effectively damping the initial pressure oscillation. Considering
this argument, it is expected that the system with a smaller pore size should have the fastest dissipation of
the oscillations since the smaller size of the pores in the nanoporous media will lead to a higher specific
area and, as a consequence, bring a relatively larger portion of fluid particles in contact with the nanoporous
88

media surface. Nonetheless, the S20 model, which is the sample with the smallest pore size, displays the
strongest oscillation and the slowest dissipation. This contrast can be understood since the permeability of
the S20 sample, as shown in Fig. 35, is much lower than the other samples. As a consequence, the dissipative
term, which depends directly on the relative velocity between particles, will be smaller for the S20 model
than for the other systems. The frequency of the pressure oscillation is also worth discussing. A direct
observation is that the pressure level in the source reaches the first peak at the same time for all the samples.
Initially, the piston starts to move due to the applied pressure. The movement creates an initial pressure
wave in the system. The speed of the wave is equal to the speed of the sound in the system since there are
no obstacles.
253
The speed of sound 𝑐 can be calculated by
𝑐 = ¡
𝐾
𝜌
77
where 𝐾 and 𝜌 are the bulk modulus and the density of the fluid, respectively. As these two quantities and
the length between the right end of the nanoporous media and the right piston are the same for all the
systems, the first peak of the pressure level in the source is expected to occur at the same time for all systems.
One can also note that the pressure level at the peaks for the samples has a trend: smaller pore size induces
a higher pressure level. Indeed, the sample with smaller pore sizes has much lower permeability. That leads
to much less fluid entering the nanoporous media during the pressure oscillation process and thus less
dissipation of the source pressure wave and generation of a stronger reflected wave from the nanoporous
media surface. As a result, the pressure level of the first peak increases with decreasing pore size. It can
also be observed that the frequency of the pressure oscillation is about the same for each system over time.
This is because, after the first peak, the system behaves similarly to a spring, whose frequency is not
dependent on the length of the source or the sink. However, the frequency for different systems slightly
differs from each other, and the complex geometries of the nanoporous media prevent us from establishing
a clear trend in the pore size.
89

Figure 35. Permeability dependence on pore size compared with predictions from the corrected Hagen-
Poiseuille Law.
6.4.2 Boundary Condition at the Solid-Fluid Interface
As presented in the background section, there are several approaches to setting the boundary
conditions between solid and fluid phases. In this work, we applied the boundary condition as proposed by
Xia et al.,
51
using a denser configuration for the solid walls. This boundary condition is effective at
preventing the penetration of fluid particles into the solid walls of the nanoporous media while not freezing
the particles in close contact with the walls. It is argued that, in the nanometer-sized pores, the surface slip
might be an important physical phenomenon to observe.
47,213,218
However, by looking at Fig. 32, we
observed that the velocity at the proximity of the solid walls is nearly zero. Therefore, the boundary
condition implemented effectively promotes a non-slip fluid-solid interface.
6.4.3 Pore Size Effect on Permeability Derived from Hagen-Poiseuille Law
The contrast in fluid velocity as a function of pore size is an important result of this work. The fluid
velocity profiles produced by the action of the pressure difference between the source and the sink are quite
different among the models. As shown in Figs. 31, 32, and 34, the overall flow velocity in the S80 model,
at a comparatively similar flow time, is almost 10 times larger than that in the S20 model, indicating that
the pore size has a significant impact on the overall flow velocity. One can note that the contrast of flow
velocity is stronger in samples with large pore sizes, i.e., samples S60 and S80, while it is challenging to
distinguish it in samples with small pore sizes, especially in the S20 sample. To understand this result, we
90

turn to Fig. 32(e), which shows the probabilities of particle velocities in the different models. It can be
observed that the velocity distribution for S80 has a peak around 0.25, while the peak for S20 is only 0.02.
Besides, the probability distribution for the S80 model has a strong bias and long tail to the right-hand side
of the distribution, while the distribution is almost symmetric with respect to the peak for S20. This means
that in the sample with small pore sizes, the advective fluid flow, driven by the applied pressure difference,
is small compared to the convective fluid flow component, driven by the isotropic thermodynamic self-
diffusion of the fluid particles. In contrast, in the S60 and S80 models, the convective fluid flow velocity,
indicated by the probability distribution peak, is comparable to the advective flow, indicated by the width-
at-half-maximum of the probability distribution, making the blockages and channels in the nanoporous
media clearly delineated.
The nanoporous media can be seen as a collection of tubes with complex geometries. To better
understand the effect of nanoporous media on the fluid flow, we study the flow properties in single
cylindrical channels with the same diameters as the average pore size in the nanoporous media. The sample
we designed is shown in Fig. 36(a). Cylindrical channels with diameters 𝑑 equal to 4, 8, 12, and 16 𝜎 (3.41,
6.82, 10.22, and 13.63 nm) and a uniform length 80 𝜎 (68.16 nm) are made, and models are prepared
following the same setup as for the nanoporous media as shown in Fig. 36(a). In the figure, the view shows
a cut displaying only half of the system along the direction perpendicular to the page. Two pistons are
located at the left- and right-hand sides of the channels. The system is then filled with liquid particles. We
apply two different pressure differences, 15.31
}
{
<
and 30.62
}
{
<
on each system and observe the flow
induced through the channel. For each channel, the radius is divided into bins of size 0.5 𝜎 (0.43 nm).
Particle velocities are averaged in each bin over 30 simulation frames in the steady-state regime. The
velocity profiles along the radius direction for different samples are plotted in Fig. 36(b). The profiles
indicate parabolic shapes with particle velocities displaying a maximum at the center of the channel and
tending to zero at the surface of the channel. The latter indicates that the boundary condition resulting from
the use of higher density mDPD particles in the solid phase generates an effective non-slip interface.
91

Figure 36. Fluid flow in cylindrical channels of different diameters. (a) Channel nanoporous model where
the view is a sliced system at mid-length. (b) Velocity profiles of fluid flow in channels with different radii
from 2 to 8 𝜎 under two pressure differences, 15.31 and 30.62
}
{
<
. The theoretical prediction based on the
Hagen-Poiseuille law is also plotted.
It is interesting to compare the calculated channel velocity profiles shown in Fig. 36 against
theoretical predictions. For the rectilinear flow of Newtonian fluids flowing through a cylindrical tube, the
equation of motion can be written as:
254

𝜕
𝜕𝑟
'𝑟𝜇
𝜕𝑢
𝑥
𝜕𝑟
( =
𝜕𝑝
𝜕𝑥
𝑟 69
where 𝑟 is the radial coordinate, 𝑢
+
(𝑟) is the fluid velocity, and 𝑝 is the local pressure. In this case, the fluid
flow is along the x-direction. Assuming non-slip boundary condition and a maximum velocity in the axis
of the tube,
𝑢
𝑥
(𝑟 = 𝑅) =
𝜕𝑢
𝑥
𝜕𝑟
+
𝑟=0
= 0 70
The velocity profile in the tube can be given as:
𝑢
𝑥
(𝑟) =
𝑅
2
4𝜇
'1−,
𝑟
𝑅
-
2
(
𝜕𝑝
𝜕𝑥

71
As the pressure gradient is constant in the channel,
z
+
= 15.31
}
{
<
or 30.62
}
{
<
for the two pressure
differences applied. Using the dynamic viscosity estimated previously, 6.25
[
{~
, we can use Eq. (71) to plot
92

the theoretical velocity distribution for the different channels. The results are shown in the solid curves
plotted in Fig. 36(b), which show an excellent agreement between simulation results and theoretical
predictions.
With the reference analysis of the flow in the cylindrical channels, we can now better analyze the
flow in the nanoporous media. It is worthwhile to compare the nanoporous media to the cylindrical tubes
to understand the calculated relationship between pore size and permeability. The derivative of the pressure
can be assumed to be a well-defined gradient in the nanoporous media for the cylindrical channel. This is
demonstrated in Fig. 30(e). The integration of Eq. (71) over the entire cross-section leads to the equation
for volumetric flow as:
𝑄 = 2𝜋𝑟
=
K
𝑅
)
4𝜇
b1−i
𝑟
𝑅
j
)
d
∆𝑝
𝑙
d𝑟 =
𝜋∆𝑝𝑅
T
8𝜇𝑙
72
which is also known as the Hagen-Poiseuille law.
254

The resulting volume flux is given by:
𝑞 =
𝑄
𝜋𝑅
)
=
∆𝑝𝑅
)
8𝜇𝑙
73
Substituting 𝑞 into Eq. (73), the permeability can be expressed as:
𝑘 =
𝑅
)
8

74
As the permeability in the simulation is evaluated taking the whole cross-section area, the porosity 𝜙 must
have a negative effect on the permeability. Finally, the permeability of the nanoporous system should have
the form:
𝑘 =
𝜙
8
b
𝐶
yz
𝐷
2
d
)

75
where 𝐷 is the average pore size of the nanoporous media, and 𝐶
yz
is a constant related to the nanoporous
morphology. The permeability as a function of pore size calculated from Eq. (75) with the 𝐶
yz
constant set
at the optimal value of 0.85 is plotted in Fig. 35 along the permeability values obtained from the simulations.
93

The model shows excellent agreement with the simulation results, indicating that a model derived from
cylindrical tubes can be used to successfully estimate the pore-size-dependency of the permeability of
nanoporous media with complex topology. It should be noted that 𝐶
yz
is a constant smaller than 1, which
indicates that the nanoporous structure will create extra resistance compared to an all-cylindrical-tube
model due to its complex geometry. It is interesting to compare Eq. (75) to the model obtained by Koponen
et al.
255
In that model, the permeability of a porous structure is expressed as:
𝑘′ =
𝜙
H
𝑐𝜏
)
𝑆
)

76
where 𝜏 is the tortuosity, 𝑆 is the specific area and 𝑐 is the Kozeny coefficient.
256
One can also invoke the
unit cell model established by Gibson and Ashby,
202
where we can relate the relative density 𝜌 to the size
of the ligaments by 𝜌 ∝i
F
]
j
)
, where 𝑑 is the ligament diameter.
248
Assuming that in our BNP sample the
relation between pore size and porosity should be:
𝜙 ∝b
𝐷
𝑙
d
)

77
If we combine Eq. (77) and Eq. (75), we can obtain a similar expression to Eq. (76) yet following
a quadratic relationship of the permeability with the porosity: 𝑘 ∝𝜙
)
. The shape factor in our expression
could include the effects of tortuosity and specific area. Our model emphasizes that the relationship between
permeability and average pore size is parabolic, while other geometrical factors are fixed.
In the current study, we use BNP models with complex morphological features that are
representative of the complex topology of shale rocks. The samples we create are strictly identical in terms
of geometrical features other than the average pore size. Our results indicate that the dependence of
permeability on pore size for such complex BNP structures is akin to the dependence found for simple
cylindrical channels. Our method to generate stochastic BNP structures is able to control different geometric
features benefiting future studies, e.g., the effects of other geometrical features, such as porosity and genus,
on the flow properties. A combination of all these understandings will then provide us with a framework to
predict the hydraulic properties of nanoporous media.
94

6.5 Conclusions
In this work, we used mDPD to investigate the fluid flow process through BNP media. Simulations
considered average pore size in the range from 4 𝜎 to 16 𝜎, where 𝜎 = 8.52 Å. Geometrical features such
as the nanoporous topology and porosity level (𝜌 = 0.65) are kept constant while the focus is placed on the
effect of average nanopore size on permeability. Simulations of flow on nanosized cylindrical channels
indicate that a no-slip boundary condition can be realized by using denser mDPD walls and that the Hagen-
Poiseuille law can describe well the fluid flow law if a suitable shape factor is added to account for the
resistance of the tortuous nature of the nanoporous media. Future work should clarify the effects of other
important geometrical features, such as porosity and genus to allow the prediction of hydraulic behavior in
nanoporous media with arbitrary feature sizes and shapes.

95

Chapter 7
Conclusions
This thesis focuses intensively on BNP materials and their properties. It starts with the generation
of the atomistic model of a BNP structure. After the development of an efficient method to generate BNP
structures on an arbitrary non-cubic cell with PBCs, we produced CuZr MGs with BNP structures and use
MD simulations to study the mechanical properties of BNPMGs, which is a novel class of materials that
have been successfully synthesized experimentally. Finally, we generate BNP media and study the
properties of fluid flow going through them using the mDPD method. The following is a summary of each
part of the work.
i) We created an efficient method to generate periodic BNP structures in general parallelepiped
cells of different length sizes and shapes. A level set value can be adjusted to have structures with different
relative densities. To verify the method, we have calculated mean curvatures, Gaussian curvatures, and
genus in periodic cubic, rectangular prism, and general parallelepiped cells and found an excellent
agreement with theoretical predictions. Compared to the traditionally used phased field simulation of the
spinodal decomposition process, this method provides an efficient alternative for the generation of BNP
models for atomistic simulations.
ii) A bulk of novel results and insights have been obtained from the study of the mechanical
properties of CuZr BNPMG using MD simulations. The study of CuZr BNPMG with 0.45 relative density
suggests that ligaments initially aligned with the loading direction bear most of the load and are responsible
for the plasticity generated immediately beyond the yield point. Ligaments in other orientations also
contribute to the overall behavior of the BNPMG as they bend over ligament joints and progressively align
with the loading direction. The reorientation of ligaments occurring mostly in the elastic region can
contribute notably to the enhanced ductility of BNPMG as it occurs concurrently with the necking of aligned
ligaments. As ligaments align with the loading direction, they are able to enhance the load support, which
96

is transferred from ligaments with fully developed necking. The synergistic combination of randomly
oriented ligaments and the ductile nature inherited from their thin cross-sections results in an anomalous
ductile behavior for a nominally brittle MG composition. We varied the relative density, composition, and
system size of CuZr BNPMG to then investigate the dependence of mechanical properties on those variables.
While the samples have self-similar structures, the relative density dictates the surface morphology and
ligament size. We found that the system size of BNPMG has little impact on the mechanical properties.
However, the relative density and Cu concentration have significant effects. A brittle-to-ductile transition
occurs at a relative density 0.7. BNPMGs with large relative density, e.g., 0.9, fail by developing a critical
shear band, which is nucleated at large pores. While those with small relative density fail by progressive
necking of ligaments. We then employed the GP method to uncover universal scaling laws of the Young’s
modulus and UTS as a function of relative density and composition and compared their predictions with
that of existing laws. The scaling law based on Gibson and Ashby's theory considering the effect of
ligaments, can be applied well to the Young’s modulus of CuZr NPMGs. However, the existing linear
scaling law for UTS is not able to accurately predict the simulation results of this work. Relationships
uncovered by GPSR are able to find succinct expressions to describe accurately the dependence of the
mechanical properties on relative density and Cu concentration. These findings shed light on the mechanical
behavior of NPMGs. The successful utilization of GPSR in this work demonstrates the great potential of
AI technologies in the study of NPMGs. The results are expected to motivate further experimental and
modeling research, which will pave the way for practical applications of these materials.
iii) Finally, We conduct mDPD simulations to investigate the fluid flow process through BNP
media. The work considered average pore size in the range from 4 𝜎 to 16 𝜎 , where 𝜎 = 8.52 Å .
Geometrical features such as the nanoporous topology and porosity level (𝜌 = 0.65) are kept constant while
the focus was placed on the effect of average nanopore size on permeability. A no-slip boundary condition
was realized by using denser mDPD walls based on a close examination of the flow in nanosized cylindrical
97

channels. It was found that the Hagen-Poiseuille law can describe well the fluid flow with an additional
shape factor that is responsible for the resistance of the tortuous nature of the nanoporous media.

98

References
1. Fratzl, P. & Weinkamer, R. Nature’s hierarchical materials. Prog. Mater. Sci. 52, 1263–1334 (2007).
2. Currey, J. D. Hierarchies in Biomineral Structures. Science (80-. ). 309, 253–254 (2005).
3. Rintoul, M. D. et al. Structure and transport properties of a porous magnetic gel via x-ray
microtomography. Phys. Rev. E 54, 2663–2669 (1996).
4. Joo, S. H. et al. Ordered nanoporous arrays of carbon supporting high dispersions of platinum
nanoparticles. Nature 412, 169–172 (2001).
5. You, T., Niwa, O., Tomita, M. & Hirono, S. Characterization of Platinum Nanoparticle-Embedded
Carbon Film Electrode and Its Detection of Hydrogen Peroxide. Anal. Chem. 75, 2080–2085 (2003).
6. BOND, G. C. & THOMPSON, D. T. Catalysis by Gold. Catal. Rev. 41, 319–388 (1999).
7. Cheng, I. C. & Hodge, A. M. Morphology, Oxidation, and Mechanical Behavior of Nanoporous Cu
Foams. Adv. Eng. Mater. 14, 219–226 (2012).
8. Cheng, I. C. & Hodge, A. M. High temperature morphology and stability of nanoporous Ag foams.
J. Porous Mater. 21, 467–474 (2014).
9. Biener, J., Hodge, A. M., Hamza, A. V., Hsiung, L. M. & Satcher, J. H. Nanoporous Au: A high
yield strength material. J. Appl. Phys. 97, 024301 (2005).
10. Sun, X.-Y., Xu, G.-K., Li, X., Feng, X.-Q. & Gao, H. Mechanical properties and scaling laws of
nanoporous gold. J. Appl. Phys. 113, 023505 (2013).
11. Farkas, D., Caro, A., Bringa, E. & Crowson, D. Mechanical response of nanoporous gold. Acta
mater. 61, 3249–3256 (2013).
12. Jin, H.-J. et al. Deforming nanoporous metal: Role of lattice coherency. Acta Mater. 57, 2665–2672
(2009).
13. Mameka, N., Wang, K., Markmann, J., Lilleodden, E. T. & Weissmüller, J. Nanoporous Gold—
Testing Macro-scale Samples to Probe Small-scale Mechanical Behavior. Mater. Res. Lett. 4, 27–
36 (2016).
14. Yeh, W. J. & Chava, S. Fabrication of metallic nanoporous films by dealloying. J. Vac. Sci. Technol.
B Microelectron. Nanom. Struct. 27, 923 (2009).
15. Erlebacher, J., Aziz, M. J., Karma, A., Dimitrov, N. & Sieradzki, K. Evolution of nanoporosity in
dealloying. Nature 410, 450–453 (2001).
16. Juarez, T. & Hodge, A. M. Synthesis of Nanoporous Gold Tubes. Adv. Eng. Mater. 18, 65–69 (2016).
17. Chen, L.-Y., Yu, J.-S., Fujita, T. & Chen, M.-W. Nanoporous Copper with Tunable Nanoporosity
for SERS Applications. Adv. Funct. Mater. 19, 1221–1226 (2009).
18. Hodge, A. M., Hayes, J. R., Caro, J. A., Biener, J. & Hamza, A. V. Characterization and Mechanical
Behavior of Nanoporous Gold. Adv. Eng. Mater. 8, 853–857 (2006).
19. Qian, L. H., Yan, X. Q., Fujita, T., Inoue, A. & Chen, M. W. Surface enhanced Raman scattering of
nanoporous gold: Smaller pore sizes stronger enhancements. Appl. Phys. Lett. 90, 153120 (2007).
20. Dan, Z., Qin, F., Sugawara, Y., Muto, I. & Hara, N. Fabrication of nanoporous copper by dealloying
amorphous binary Ti–Cu alloys in hydrofluoric acid solutions. Intermetallics 29, 14–20 (2012).

99

21. Zhao, C., Qi, Z., Wang, X. & Zhang, Z. Fabrication and characterization of monolithic nanoporous
copper through chemical dealloying of Mg–Cu alloys. Corros. Sci. 51, 2120–2125 (2009).
22. Sun, L., Chien, C.-L. & Searson, P. C. Fabrication of Nanoporous Nickel by Electrochemical
Dealloying. Chem. Mater. 16, 3125–3129 (2004).
23. Jin, T. et al. Nanoporous Copper Metal Catalyst in Click Chemistry: Nanoporosity-Dependent
Activity without Supports and Bases. Adv. Synth. Catal. 353, 3095–3100 (2011).
24. Badwe, N., Chen, X. & Sieradzki, K. Mechanical properties of nanoporous gold in tension. Acta
Mater. 129, 251–258 (2017).
25. Li, R. & Sieradzki, K. Ductile-brittle transition in random porous Au. Phys. Rev. Lett. 68, 1168–
1171 (1992).
26. Dursun, A., Pugh, D. V. & Corcoran, S. G. Dealloying of Ag-Au Alloys in Halide-Containing
Electrolytes. J. Electrochem. Soc. 150, B355 (2003).
27. Chen, L. Y., Fujita, T. & Chen, M. W. Biofunctionalized nanoporous gold for electrochemical
biosensors. Electrochim. Acta 67, 1–5 (2012).
28. Biener, J., Hodge, A. M. & Hamza, A. V. Microscopic failure behavior of nanoporous gold. Appl.
Phys. Lett. 87, 121908 (2005).
29. Ngô, B. N. D. et al. Anomalous compliance and early yielding of nanoporous gold. Acta Mater. 93,
144–155 (2015).
30. CROWSON, D., FARKAS, D. & CORCORAN, S. Geometric relaxation of nanoporous metals: The
role of surface relaxation. Scr. Mater. 56, 919–922 (2007).
31. Griffiths, E., Bargmann, S. & Reddy, B. D. Elastic behaviour at the nanoscale of innovative
composites of nanoporous gold and polymer. Extrem. Mech. Lett. 17, 16–23 (2017).
32. Guillotte, M., Godet, J. & Pizzagalli, L. A fully molecular dynamics-based method for modeling
nanoporous gold. Comput. Mater. Sci. 161, 135–142 (2019).
33. Cahn, J. W. Phase separation by spinodal decomposition in isotropic systems. J. Chem. Phys. 42,
93–99 (1965).
34. KLEMENT, W., WILLENS, R. H. & DUWEZ, P. Non-crystalline Structure in Solidified Gold–
Silicon Alloys. Nature 187, 869–870 (1960).
35. Yavari, A. R. et al. The glass transition of bulk metallic glasses studied by real-time diffraction in
transmission using high-energy synchrotron radiation. Mater. Sci. Eng. A 375–377, 709–712 (2004).
36. Shan, Z. W. et al. Plastic flow and failure resistance of metallic glass: Insight from in situ
compression of nanopillars. Phys. Rev. B 77, 155419 (2008).
37. SCHUH, C., HUFNAGEL, T. & RAMAMURTY, U. Mechanical behavior of amorphous alloys.
Acta Mater. 55, 4067–4109 (2007).
38. Greer, A. L. & Ma, E. Bulk Metallic Glasses: At the Cutting Edge of Metals Research. MRS Bull.
32, 611–619 (2007).
39. Şopu, D. et al. Structure-property relationships in nanoporous metallic glasses. Acta Mater. 106,
199–207 (2016).
40. Jiao, W. et al. Tunable Nanoporous Metallic Glasses Fabricated by Selective Phase Dissolution and
Passivation for Ultrafast Hydrogen Uptake. Chem. Mater. 29, 4478–4483 (2017).

100

41. Larsson, M. Global Energy Transformation. (Palgrave Macmillan UK, 2009).
doi:10.1057/9780230244092
42. Capuano, L. International Energy Outlook 2020 ( IEO2020 ) United States milestones in meeting
global energy consumption. US Energy Information Administration 2020, (2020).
43. Ren, Y. & Stein, D. Slip-enhanced electrokinetic energy conversion in nanofluidic channels.
Nanotechnology 19, 195707 (2008).
44. Ghadiri, M., Fakhri, S. & Shirazian, S. Modeling and CFD Simulation of Water Desalination Using
Nanoporous Membrane Contactors. Ind. Eng. Chem. Res. 52, 3490–3498 (2013).
45. Schlapbach, L. & Züttel, A. Hydrogen-storage materials for mobile applications. Nature 414, 353–
358 (2001).
46. Ally, J., Molla, S. & Mostowfi, F. Condensation in Nanoporous Packed Beds. Langmuir 32, 4494–
4499 (2016).
47. Sokhan, V. P., Nicholson, D. & Quirke, N. Fluid flow in nanopores: An examination of
hydrodynamic boundary conditions. J. Chem. Phys. 115, 3878–3887 (2001).
48. Kanellopoulos, N. K. Recent Advances in Gas Separation by Microporous Ceramic Membranes,
Volume 6. (2000).
49. Churaev, N. V., Lyklema, J. & Churaeva, M. N. Liquid and Vapor Flows in Porous Bodies.
(Routledge, 2018). doi:10.1201/9780203748794
50. ZUO, H., DENG, S. & LI, H. Limitations of Lattice Boltzmann Modeling of Micro‐Flows in
Complex Nanopores. Acta Geol. Sin. - English Ed. 93, 1808–1822 (2019).
51. Xia, Y. et al. Many-body dissipative particle dynamics modeling of fluid flow in fine-grained
nanoporous shales. Phys. Fluids 29, (2017).
52. Chen, I. J., Yin, D. & MacKerell, A. D. Combinedab initio/empirical approach for optimization of
Lennard-Jones parameters for polar-neutral compounds. J. Comput. Chem. 23, 199–213 (2002).
53. Daw, M. S. & Baskes, M. I. Embedded-atom method: Derivation and application to impurities,
surfaces, and other defects in metals. Phys. Rev. B 29, 6443–6453 (1984).
54. Rota, G.-C. Handbook of stochastic methods. Adv. Math. (N. Y). 55, 101 (1985).
55. Liu, M. B., Liu, G. R., Zhou, L. W. & Chang, J. Z. Dissipative Particle Dynamics (DPD): An
Overview and Recent Developments. Arch. Comput. Methods Eng. 22, 529–556 (2015).
56. Groot, R. D. & Warren, P. B. Dissipative particle dynamics: Bridging the gap between atomistic
and mesoscopic simulation. J. Chem. Phys. 107, 4423–4435 (1997).
57. Balk, T. J., Eberl, C., Sun, Y., Hemker, K. J. & Gianola, D. S. Tensile and compressive
microspecimen testing of bulk nanoporous gold. JOM 61, 26–31 (2009).
58. Biener, J. et al. Nanoporous Plasmonic Metamaterials. Adv. Mater. 20, 1211–1217 (2008).
59. Cheng, I.-C. & Hodge, A. M. Strength scale behavior of nanoporous Ag, Pd and Cu foams. Scr.
Mater. 69, 295–298 (2013).
60. Crowson, D. A., Farkas, D. & Corcoran, S. G. Mechanical stability of nanoporous metals with small
ligament sizes. Scr. Mater. 61, 497–499 (2009).
61. Hodge, A. M. et al. Scaling equation for yield strength of nanoporous open-cell foams. Acta Mater.
55, 1343–1349 (2007).

101

62. Jin, H.-J. et al. Nanoporous Au−Pt Alloys As Large Strain Electrochemical Actuators. Nano Lett.
10, 187–194 (2010).
63. Kramer, D., Viswanath, R. N. & Weissmüller, J. Surface-Stress Induced Macroscopic Bending of
Nanoporous Gold Cantilevers. Nano Lett. 4, 793–796 (2004).
64. Biener, J. et al. Surface-chemistry-driven actuation in nanoporous gold. Nat. Mater. 8, 47–51 (2009).
65. Snyder, J., Fujita, T., Chen, M. W. & Erlebacher, J. Oxygen reduction in nanoporous metal–ionic
liquid composite electrocatalysts. Nat. Mater. 9, 904–907 (2010).
66. Ding, Y., Chen, M. & Erlebacher, J. Metallic Mesoporous Nanocomposites for Electrocatalysis. J.
Am. Chem. Soc. 126, 6876–6877 (2004).
67. Wittstock, A., Zielasek, V., Biener, J., Friend, C. M. & Bäumer, M. Nanoporous Gold Catalysts for
Selective Gas-Phase Oxidative Coupling of Methanol at Low Temperature. Science (80-. ). 327,
319–322 (2010).
68. Qian, L. H., Inoue, A. & Chen, M. W. Large surface enhanced Raman scattering enhancements from
fracture surfaces of nanoporous gold. Appl. Phys. Lett. 92, 093113 (2008).
69. Lang, X. Y., Chen, L. Y., Guan, P. F., Fujita, T. & Chen, M. W. Geometric effect on surface
enhanced Raman scattering of nanoporous gold: Improving Raman scattering by tailoring ligament
and nanopore ratios. Appl. Phys. Lett. 94, 213109 (2009).
70. Hanks, D. F. et al. Nanoporous membrane device for ultra high heat flux thermal management.
Microsystems Nanoeng. 4, 1 (2018).
71. Kunugi, T., Muko, K. & Shibahara, M. Ultrahigh heat transfer enhancement using nano-porous layer.
Superlattices Microstruct. 35, 531–542 (2004).
72. Wadley, H. Fabrication and structural performance of periodic cellular metal sandwich structures.
Compos. Sci. Technol. 63, 2331–2343 (2003).
73. Evans, A. G. et al. Concepts for enhanced energy absorption using hollow micro-lattices. Int. J.
Impact Eng. 37, 947–959 (2010).
74. Surani, F. B., Kong, X. & Qiao, Y. Two-staged sorption isotherm of a nanoporous energy absorption
system. Appl. Phys. Lett. 87, 251906 (2005).
75. Biener, J. et al. Size Effects on the Mechanical Behavior of Nanoporous Au. Nano Lett. 6, 2379–
2382 (2006).
76. Ruestes, C. J., Farkas, D., Caro, A. & Bringa, E. M. Hardening under compression in Au foams.
Acta Mater. 108, 1–7 (2016).
77. Ruestes, C. J., Schwen, D., Millán, E. N., Aparicio, E. & Bringa, E. M. Mechanical properties of Au
foams under nanoindentation. Comput. Mater. Sci. 147, 154–167 (2018).
78. Beets, N. & Farkas, D. Mechanical Response of Au Foams of Varying Porosity from Atomistic
Simulations. JOM 70, 2185–2191 (2018).
79. Xiang, M. et al. Shock responses of nanoporous aluminum by molecular dynamics simulations. Int.
J. Plast. 97, 24–45 (2017).
80. Jian, W. R., Li, B., Wang, L., Yao, X. H. & Luo, S. N. Shock response of open-cell nanoporous Cu
foams: Effects of porosity and specific surface area. J. Appl. Phys. 118, 165902 (2015).
81. Xian, Y., Li, J., Wu, R. & Xia, R. Softening of nanocrystalline nanoporous platinum: A molecular
dynamics simulation. Comput. Mater. Sci. 143, 163–169 (2018).
102

82. JIANG, Z., CHEN, W. & BURKHART, C. Efficient 3D porous microstructure reconstruction via
Gaussian random field and hybrid optimization. J. Microsc. 252, 135–148 (2013).
83. Torquato, S. & Haslach, H. Random Heterogeneous Materials: Microstructure and Macroscopic
Properties. Appl. Mech. Rev. (2002). doi:10.1115/1.1483342
84. Soyarslan, C., Bargmann, S., Pradas, M. & Weissmüller, J. 3D stochastic bicontinuous
microstructures: Generation, topology and elasticity. Acta Mater. 149, 326–340 (2018).
85. Zhu, B. et al. Short Review on Porous Metal Membranes—Fabrication, Commercial Products, and
Applications. Membranes (Basel). 8, 83 (2018).
86. Sun, X.-Y., Li, Q., Gu, Y. & Feng, X.-Q. Mechanical properties of bioinspired bicontinuous
nanocomposites. Comput. Mater. Sci. 80, 71–78 (2013).
87. Huang, J. et al. Gradient nanoporous gold: a novel surface-enhanced Raman scattering substrate.
RSC Adv. 7, 15747–15753 (2017).
88. Lian, L., Yao, Y., Liu, Y. & Fang, X. A segmental dealloying for fabricating the gradient nanoporous
metal materials. J. Porous Mater. 24, 211–215 (2017).
89. Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO–the Open
Visualization Tool. Model. Simul. Mater. Sci. Eng. 18, 015012 (2010).
90. Naylor, M. Golden, , and π Flowers: A Spiral Story. Math. Mag. 75, 163–172 (2002).
91. Swinbank, R. & James Purser, R. Fibonacci grids: A novel approach to global modelling. Q. J. R.
Meteorol. Soc. 132, 1769–1793 (2006).
92. Saff, E. B. & Kuijlaars, A. B. J. Distributing many points on a sphere. Math. Intell. 19, 5–11 (1997).
93. Chukkapalli, G., Karpik, S. R. & Ethier, C. R. A Scheme for Generating Unstructured Grids on
Spheres with Application to Parallel Computation. J. Comput. Phys. 149, 114–127 (1999).
94. JOHANSSON, B. & LUKACS, A. FIBONACCI OPTICAL LATTICE: A FRACTAL
BIOMODULATION DEVICE. Fractals 27, 1950039 (2019).
95. Coates, S., Smerdon, J. A., McGrath, R. & Sharma, H. R. A molecular overlayer with the Fibonacci
square grid structure. Nat. Commun. 9, 3435 (2018)
96. Qin, L. & Xu, Y. Fibonacci helps to evacuate from a convex region in a grid network. J. Comb.
Optim. 34, 398–413 (2017).
97. Stuckner, J., Frei, K., McCue, I., Demkowicz, M. J. & Murayama, M. AQUAMI: An open source
Python package and GUI for the automatic quantitative analysis of morphologically complex
multiphase materials. Comput. Mater. Sci. 139, 320–329 (2017).
98. Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO-the Open
Visualization Tool. Model. Simul. Mater. Sci. Eng. 18, (2010).
99. Albin, E., Knikker, R., Xin, S., Paschereit, C. O. & D’Angelo, Y. Computational Assessment of
Curvatures and Principal Directions of Implicit Surfaces from 3D Scalar Data. in Lecture Notes in
Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in
Bioinformatics) 10521 LNCS, 1–22 (2017).
100. Teubner, M. Level Surfaces of Gaussian Random Fields and Microemulsions. Europhys. Lett. 14,
403–408 (1991).
101. Fujita, T., Qian, L.-H., Inoke, K., Erlebacher, J. & Chen, M.-W. Three-dimensional morphology of
nanoporous gold. Appl. Phys. Lett. 92, 251902 (2008).
103

102. Mangipudi, K. R., Radisch, V., Holzer, L. & Volkert, C. A. A FIB-nanotomography method for
accurate 3D reconstruction of open nanoporous structures. Ultramicroscopy 163, 38–47 (2016).
103. Mangipudi, K. R., Epler, E. & Volkert, C. A. Topology-dependent scaling laws for the stiffness and
strength of nanoporous gold. Acta Mater. 119, 115–122 (2016).
104. Chen, Y. K. et al. Morphological and topological analysis of coarsened nanoporous gold by x-ray
nanotomography. Appl. Phys. Lett. 96, 043122 (2010).
105. Chen, C. Q., Pei, Y. T. & De Hosson, J. T. M. Effects of size on the mechanical response of metallic
glasses investigated through in situ TEM bending and compression experiments. Acta Mater. 58,
189–200 (2010).
106. Wu, X. & Zhu, Y. Heterogeneous materials: a new class of materials with unprecedented mechanical
properties. Mater. Res. Lett. 5, 527–532 (2017).
107. Lu, K. Making strong nanomaterials ductile with gradients. Science (80-. ). 345, 1455–1456 (2014).
108. Wei, Y. et al. Evading the strength–ductility trade-off dilemma in steel through gradient hierarchical
nanotwins. Nat. Commun. 5, 3580 (2014).
109. Wu, X., Jiang, P., Chen, L., Yuan, F. & Zhu, Y. T. Extraordinary strain hardening by gradient
structure. Proc. Natl. Acad. Sci. 111, 7197–7201 (2014).
110. Ning, Y., Yao, Z., Guo, H. & Fu, M. W. Structural-gradient-materials produced by gradient
temperature heat treatment for dual-property turbine disc. J. Alloys Compd. 557, 27–33 (2013).
111. Fang, T. H., Li, W. L., Tao, N. R. & Lu, K. Revealing Extraordinary Intrinsic Tensile Plasticity in
Gradient Nano-Grained Copper. Science (80-. ). 331, 1587–1590 (2011).
112. Russew, K. & Stojanova, L. Glassy Metals. Glassy Metals (Springer Berlin Heidelberg, 2016).
doi:10.1007/978-3-662-47882-0
113. Chen, M. A brief overview of bulk metallic glasses. NPG Asia Mater. 3, 82–90 (2011).
114. Chen, M. Mechanical Behavior of Metallic Glasses: Microscopic Understanding of Strength and
Ductility. Annu. Rev. Mater. Res. 38, 445–469 (2008).
115. Wang, W. H. Bulk Metallic Glasses with Functional Physical Properties. Adv. Mater. 21, 4524–
4544 (2009).
116. Ding, S. et al. Combinatorial development of bulk metallic glasses. Nat. Mater. 13, 494–500 (2014).
117. Liu, P. S. & Liang, K. M. Functional materials of porous metals made by P/M, electroplating and
some other techniques. J. Mater. Sci. 36, 5059–5072 (2001).
118. Brothers, A. H. & Dunand, D. C. Amorphous metal foams. Scr. Mater. 54, 513–520 (2006).
119. Williams, D. E., Newman, R. C., Song, Q. & Kelly, R. G. Passivity breakdown and pitting corrosion
of binary alloys. Nature 350, 216–219 (1991).
120. Erlebacher, J. & Sieradzki, K. Pattern formation during dealloying. Scr. Mater. 49, 991–996 (2003).
121. Chen, X. H. et al. A Porous Bulk Metallic Glass with Unidirectional Opening Pores. Electrochem.
Solid-State Lett. 10, E21 (2007).
122. Rösler, J. & Mukherji, D. Design of nanoporous superalloy membranes for functional applications.
Adv. Eng. Mater. 5, 916–918 (2003).
123. Gebert, A., Kündig, A. A., Schultz, L. & Hono, K. Selective electrochemical dissolution in two-
phase La-Zr-Al-Cu-Ni metallic glass. Scr. Mater. (2004). doi:10.1016/j.scriptamat.2004.07.021
104

124. Wada, T., Inoue, A. & Greer, A. L. Enhancement of room-temperature plasticity in a bulk metallic
glass by finely dispersed porosity. Appl. Phys. Lett. 86, 251907 (2005).
125. Adibi, S., Branicio, P. S. & Joshi, S. P. Suppression of Shear Banding and Transition to Necking
and Homogeneous Flow in Nanoglass Nanopillars. Sci. Rep. 5, 15611 (2015).
126. Liu, C. & Branicio, P. S. Efficient generation of non-cubic stochastic periodic bicontinuous
nanoporous structures. Comput. Mater. Sci. 169, 109101 (2019).
127. Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 117,
1–19 (1995).
128. Thompson, A. P., Plimpton, S. J. & Mattson, W. General formulation of pressure and stress tensor
for arbitrary many-body interaction potentials under periodic boundary conditions. J. Chem. Phys.
(2009). doi:10.1063/1.3245303
129. Shimizu, F., Ogata, S. & Li, J. Theory of Shear Banding in Metallic Glasses and Molecular
Dynamics Calculations. Mater. Trans. 48, 2923–2927 (2007).
130. Mendelev, M. I., Sordelet, D. J. & Kramer, M. J. Using atomistic computer simulations to analyze
x-ray diffraction data from metallic glasses. J. Appl. Phys. 102, 043501 (2007).
131. Sha, Z. D. et al. Effect of aspect ratio on the mechanical properties of metallic glasses. Scr. Mater.
93, 36–39 (2014).
132. Wu, H. A. Molecular dynamics study on mechanics of metal nanowire. Mech. Res. Commun. 33, 9–
16 (2006).
133. Zhang, Q., Li, Q. K. & Li, M. Internal stress and its effect on mechanical strength of metallic glass
nanowires. Acta Mater. 91, 174–182 (2015).
134. Zhang, Q., Li, Q. K. & Li, M. Key factors affecting mechanical behavior of metallic glass nanowires.
Sci. Rep. 7, 1–8 (2017).
135. Diao, J., Gall, K. & Dunn, M. L. Yield strength asymmetry in metal nanowires. Nano Lett. 4, 1863–
1867 (2004).
136. Zhou, X., Zhou, H., Li, X. & Chen, C. Size effects on tensile and compressive strengths in metallic
glass nanowires. J. Mech. Phys. Solids 84, 130–144 (2015).
137. Lin, W., Teng, Y., Sha, Z., Yuan, S. & Branicio, P. S. Mechanical properties of nanoporous metallic
glasses: Insights from large-scale atomic simulations. Int. J. Plast. 10.1016/j.ijplas.2019.102657
(2020). doi:10.1016/j.ijplas.2019.102657
138. Beets, N., Farkas, D. & Corcoran, S. Deformation mechanisms and scaling relations in the
mechanical response of nano-porous Au. Acta Mater. 165, 626–637 (2019).
139. Adibi, S. et al. A transition from localized shear banding to homogeneous superplastic flow in
nanoglass. Appl. Phys. Lett. 103, 1–6 (2013).
140. Adibi, S. et al. Surface roughness imparts tensile ductility to nanoscale metallic glasses. Extrem.
Mech. Lett. 5, 88–95 (2015).
141. Matthews, D. T. A., Ocelík, V., Bronsveld, P. M. & De Hosson, J. T. M. An electron microscopy
appraisal of tensile fracture in metallic glasses. Acta Mater. 56, 1762–1773 (2008).
142. Song, H. Y. et al. Atomic simulations of plastic deformation behavior of Cu 50 Zr 50 metallic glass.
J. Non. Cryst. Solids 471, 312–321 (2017).
143. Zhou, X., Wang, L. & Chen, C. Q. Strengthening mechanisms in nanoporous metallic glasses.
Comput. Mater. Sci. 155, 151–158 (2018).
105

144. Weissmüller, J., Newman, R. C., Jin, H., Hodge, A. M. & Kysar, J. W. Nanoporous Metals by Alloy
Corrosion: Formation and Mechanical Properties. MRS Bull. 34, 577–586 (2009).
145. Biener, J., Hamza, A. V. & Hodge, A. M. Deformation Behavior of Nanoporous Metals. in Micro
and Nano Mechanical Testing of Materials and Devices 118–135 (Springer US, 2008).
doi:10.1007/978-0-387-78701-5_6
146. Li, J. et al. Microstructure-sensitive mechanical properties of nanoporous gold: a molecular
dynamics study. Model. Simul. Mater. Sci. Eng. 26, 075003 (2018).
147. Kiani, A. & Fard, E. N. Fabrication of palladium coated nanoporous gold film electrode via
underpotential deposition and spontaneous metal replacement: A low palladium loading electrode
with electrocatalytic activity. Electrochim. Acta 54, 7254–7259 (2009).
148. Adibi, S., Branicio, P. S. & Ballarini, R. Compromising high strength and ductility in nanoglass-
metallic glass nanolaminates. RSC Adv. 6, 13548–13553 (2016).
149. Greer, A. L., Cheng, Y. Q. & Ma, E. Shear bands in metallic glasses. Mater. Sci. Eng. R Reports 74,
71–132 (2013).
150. Ghidelli, M. et al. Extrinsic mechanical size effects in thin ZrNi metallic glass films. Acta Mater.
90, 232–241 (2015).
151. Liu, M. C. et al. Is the compression of tapered micro- and nanopillar samples a legitimate technique
for the identification of deformation mode change in metallic glasses? Scr. Mater. 66, 817–820
(2012).
152. Chen, C. Q. et al. Intrinsic size effects in the mechanical response of taper-free nanopillars of
metallic glass. Phys. Rev. B - Condens. Matter Mater. Phys. 83, 4–7 (2011).
153. Greer, J. R. & De Hosson, J. T. M. Plasticity in small-sized metallic systems: Intrinsic versus
extrinsic size effect. Prog. Mater. Sci. 56, 654–724 (2011).
154. Şopu, D., Scudino, S., Bian, X. L., Gammer, C. & Eckert, J. Atomic-scale origin of shear band
multiplication in heterogeneous metallic glasses. Scr. Mater. 178, 57–61 (2020).
155. Sarac, B., Klusemann, B., Xiao, T. & Bargmann, S. Materials by design: An experimental and
computational investigation on the microanatomy arrangement of porous metallic glasses. Acta
Mater. 77, 411–422 (2014).
156. Wu, Z., Zhang, Y. W., Jhon, M. H., Gao, H. & Srolovitz, D. J. Nanowire failure: Long = brittle and
short = ductile. Nano Lett. 12, 910–914 (2012).
157. Yang, G. J., Xu, B., Kong, L. T., Li, J. F. & Zhao, S. Size effects in Cu 50 Zr 50 metallic glass films
revealed by molecular dynamics simulations. J. Alloys Compd. 688, 88–95 (2016).
158. Zhang, Y. et al. Characterization of the deformation behaviors under uniaxial stress for bicontinuous
nanoporous amorphous alloys. Phys. Chem. Chem. Phys. 24, 1099–1112 (2022).
159. Zhang, Y. et al. Atomistic investigation on the mechanical properties of 3D nanoporous metallic
glasses under uniaxial tension and compression. Mater. Today Commun. 27, 102460 (2021).
160. Zhang, Z. et al. Deformation behavior of a nanoporous metallic glass at room temperature. Int. J.
Plast. 152, 103232 (2022).
161. Li, J. et al. Structurally ordered nanoporous Pt–Co alloys with enhanced mechanical behaviors in
tension. Microporous Mesoporous Mater. 295, 109955 (2020).
162. Schmidt, M. & Lipson, H. Distilling Free-Form Natural Laws from Experimental Data. Science
(80-. ). 324, 81–85 (2009).
106

163. Kline, R. Cybernetics, Automata Studies, and the Dartmouth Conference on Artificial Intelligence.
IEEE Ann. Hist. Comput. 33, 5–16 (2011).
164. Sha, W. et al. Artificial Intelligence to Power the Future of Materials Science and Engineering. Adv.
Intell. Syst. 2, 1900143 (2020).
165. Agrawal, A. & Choudhary, A. Perspective: Materials informatics and big data: Realization of the
“fourth paradigm” of science in materials science. APL Mater. 4, 053208 (2016).
166. Vasudevan, R. K. et al. Materials science in the artificial intelligence age: high-throughput library
generation, machine learning, and a pathway from correlations to the underpinning physics. MRS
Commun. 9, 821–838 (2019).
167. Panapitiya, G. et al. Machine-Learning Prediction of CO Adsorption in Thiolated, Ag-Alloyed Au
Nanoclusters. J. Am. Chem. Soc. 140, 17508–17514 (2018).
168. Artrith, N., Hiller, B. & Behler, J. Neural network potentials for metals and oxides - First
applications to copper clusters at zinc oxide. Phys. status solidi 250, 1191–1203 (2013).
169. Sahu, H., Rao, W., Troisi, A. & Ma, H. Toward Predicting Efficiency of Organic Solar Cells via
Machine Learning and Improved Descriptors. Adv. Energy Mater. 8, 1801032 (2018).
170. Lu, S. et al. Accelerated discovery of stable lead-free hybrid organic-inorganic perovskites via
machine learning. Nat. Commun. 9, 3405 (2018).
171. Li, Z., Xu, Q., Sun, Q., Hou, Z. & Yin, W.-J. Thermodynamic Stability Landscape of Halide Double
Perovskites via High-Throughput Computing and Machine Learning. Adv. Funct. Mater. 29,
1807280 (2019).
172. Salmenjoki, H., Alava, M. J. & Laurson, L. Machine learning plastic deformation of crystals. Nat.
Commun. 9, 5307 (2018).
173. Kim, E. et al. Materials Synthesis Insights from Scientific Literature via Text Extraction and
Machine Learning. Chem. Mater. 29, 9436–9444 (2017).
174. Schneider, N., Lowe, D. M., Sayle, R. A. & Landrum, G. A. Development of a Novel Fingerprint
for Chemical Reactions and Its Application to Large-Scale Reaction Classification and Similarity.
J. Chem. Inf. Model. 55, 39–53 (2015).
175. Zhang, H., Moon, S. K. & Ngo, T. H. Hybrid Machine Learning Method to Determine the Optimal
Operating Process Window in Aerosol Jet 3D Printing. ACS Appl. Mater. Interfaces 11, 17994–
18003 (2019).
176. Lin, B., Hedrick, J. L., Park, N. H. & Waymouth, R. M. Programmable High-Throughput Platform
for the Rapid and Scalable Synthesis of Polyester and Polycarbonate Libraries. J. Am. Chem. Soc.
141, 8921–8927 (2019).
177. Gonzalez, T. F. Handbook of Approximation Algorithms and Metaheuristics. Handbook of
Approximation Algorithms and Metaheuristics (Chapman and Hall/CRC, 2007).
doi:10.1201/9781420010749
178. Li, W., Field, K. G. & Morgan, D. Automated defect analysis in electron microscopic images. npj
Comput. Mater. 4, 36 (2018).
179. Koza, J. Genetic programming as a means for programming computers by natural selection. Stat.
Comput. 4, 87–112 (1994).
180. Wang, Y., Wagner, N. & Rondinelli, J. M. Symbolic regression in materials science. MRS Commun.
9, 793–805 (2019).
107

181. Koza, J. R. et al. Genetic Programming IV: Routine Human-Competitive Machine Intelligence.
(2003).
182. Zhong, J., Feng, L., Cai, W. & Ong, Y.-S. Multifactorial Genetic Programming for Symbolic
Regression Problems. IEEE Trans. Syst. Man, Cybern. Syst. 50, 4492–4505 (2020).
183. Sastry, K. N. Genetic Algorithms and Genetic Programming for Multiscale Modeling: Applications
in Materials Science and Chemistry and Advances in Scalability. (2007).
184. Quaranta, G., Lacarbonara, W. & Masri, S. F. A review on computational intelligence for
identification of nonlinear dynamical systems. Nonlinear Dyn. 99, 1709–1761 (2020).
185. Prieto, A. et al. Neural networks: An overview of early research, current frameworks and new
challenges. Neurocomputing 214, 242–268 (2016).
186. Ghasemi, F., Mehridehnavi, A., Pérez-Garrido, A. & Pérez-Sánchez, H. Neural network and deep-
learning algorithms used in QSAR studies: merits and drawbacks. Drug Discov. Today 23, 1784–
1790 (2018).
187. Chopra, P., Sharma, R. K. & Kumar, M. Prediction of Compressive Strength of Concrete Using
Artificial Neural Network and Genetic Programming. Adv. Mater. Sci. Eng. 2016, 1–10 (2016).
188. Cai, W., Pacheco-Vega, A., Sen, M. & Yang, K. T. Heat transfer correlations by symbolic regression.
Int. J. Heat Mass Transf. 49, 4352–4359 (2006).
189. Langdon, W. B. & Barrett, S. J. Genetic Programming in Data Mining for Drug Discovery. in
Evolutionary Computation in Data Mining 211–235 (Springer-Verlag, 2006). doi:10.1007/3-540-
32358-9_10
190. Barmpalexis, P., Kachrimanis, K., Tsakonas, A. & Georgarakis, E. Symbolic regression via genetic
programming in the optimization of a controlled release pharmaceutical formulation. Chemom.
Intell. Lab. Syst. 107, 75–82 (2011).
191. Im, J., Rizzo, C. B., de Barros, F. P. J. & Masri, S. F. Application of genetic programming for model-
free identification of nonlinear multi-physics systems. Nonlinear Dyn. 104, 1781–1800 (2021).
192. Ward, L., Miracle, D., Windl, W., Senkov, O. N. & Flores, K. Structural evolution and kinetics in
Cu-Zr metallic liquids from molecular dynamics simulations. Phys. Rev. B 88, 134205 (2013).
193. Forouzan, R. F. G. B. A. Data Structures: A Pseudocode Approach with C. (Cengage Learning,
2005).
194. Albe, K., Ritter, Y. & Şopu, D. Enhancing the plasticity of metallic glasses: Shear band formation,
nanocomposites and nanoglasses investigated by molecular dynamics simulations. Mech. Mater. 67,
94–103 (2013).
195. Liu, C., Yuan, S. & Branicio, P. S. Bicontinuous nanoporous design induced homogenization of
strain localization in metallic glasses. Scr. Mater. 192, 67–72 (2021).
196. Wang, X. L. et al. Plasticity of a scandium-based nanoglass. Scr. Mater. 98, 40–43 (2015).
197. Peng, C. X. et al. Deformation behavior of designed dual-phase CuZr metallic glasses. Mater. Des.
168, 107662 (2019).
198. Yang, M. H., Li, J. H. & Liu, B. X. Proposed correlation of structure network inherited from
producing techniques and deformation behavior for Ni-Ti-Mo metallic glasses via atomistic
simulations. Sci. Rep. 6, 29722 (2016).

108

199. Wu, Z. W., Li, M. Z., Wang, W. H. & Liu, K. X. Correlation between structural relaxation and
connectivity of icosahedral clusters in CuZr metallic glass-forming liquids. Phys. Rev. B 88, 054202
(2013).
200. Imran, M., Hussain, F., Rashid, M., Cai, Y. & Ahmad, S. A. Mechanical behavior of Cu—Zr bulk
metallic glasses (BMGs): A molecular dynamics approach. Chinese Phys. B 22, 096101 (2013).
201. Lee, M., Lee, C. M., Lee, K. R., Ma, E. & Lee, J. C. Networked interpenetrating connections of
icosahedra: Effects on shear transformations in metallic glass. Acta Mater. 59, 159–170 (2011).
202. Gibson, L. J. & Ashby, M. F. Cellular Solids. Cellular Solids: Structure and Properties, Second
Edition (Cambridge University Press, 1997). doi:10.1017/CBO9781139878326
203. Wang, G.-F. & Feng, X.-Q. Surface effects on buckling of nanowires under uniaxial compression.
Appl. Phys. Lett. 94, 141913 (2009).
204. Winter, N., Becton, M., Zhang, L. & Wang, X. Failure Mechanisms and Scaling Laws of
Nanoporous Aluminum: A Computational Study. Adv. Eng. Mater. 18, 632–642 (2016).
205. Pia, G., Carta, M. & Delogu, F. Nanoporous Au foams: Variation of effective Young’s modulus
with ligament size. Scr. Mater. 144, 22–26 (2018).
206. Kikkinides, E. S., Kainourgiakis, M. E. & Kanellopoulos, N. K. Simulation of Gas Transport in a
“Network of Micropores”. The effect of Pore Structure on Transport Properties. in Recent Advances
in Gas Separation by Microporous Ceramic Membranes (ed. Kanellopoulos, N. K.) 297–322
(Elsevier, 2000). doi:10.1016/S0927-5193(00)80013-3
207. Bhatia, S. K. Modeling Pure Gas Permeation in Nanoporous Materials and Membranes. Langmuir
26, 8373–8385 (2010).
208. Drahushuk, L. W. & Strano, M. S. Mechanisms of Gas Permeation through Single Layer Graphene
Membranes. Langmuir 28, 16671–16678 (2012).
209. Sun, C., Zhou, R., Zhao, Z. & Bai, B. Extending the Classical Continuum Theory to Describe Water
Flow through Two-Dimensional Nanopores. Langmuir 37, 6158–6167 (2021).
210. Pozhar, L. A. & Gubbins, K. E. Transport theory of dense, strongly inhomogeneous fluids. J. Chem.
Phys. 99, 8970–8996 (1993).
211. Davis, H. T. Kinetic theory of inhomogeneous fluid: Tracer diffusion. J. Chem. Phys. 86, 1474–
1477 (1987).
212. Strong, S. E. & Eaves, J. D. The dynamics of water in porous two-dimensional crystals. J. Phys.
Chem. B 121, 189–207 (2017).
213. Yoshioka, T. et al. Molecular dynamics simulation study on the mechanisms of liquid-phase
permeation in nanopores. Sep. Purif. Technol. 220, 259–267 (2019).
214. Ghoufi, A. & Malfreyt, P. Mesoscale modeling of the water liquid-vapor interface: A surface tension
calculation. Phys. Rev. E 83, 051601 (2011).
215. Billemont, P., Coasne, B. & De Weireld, G. An Experimental and Molecular Simulation Study of
the Adsorption of Carbon Dioxide and Methane in Nanoporous Carbons in the Presence of Water.
Langmuir 27, 1015–1024 (2011).
216. Razavi, S. M. R., Razavi, S. M. J., Miri, T. & Shirazian, S. CFD simulation of CO2 capture from
gas mixtures in nanoporous membranes by solution of 2-amino-2-methyl-1-propanol and piperazine.
Int. J. Greenh. Gas Control 15, 142–149 (2013).

109

217. Pan, C., Hilpert, M. & Miller, C. T. Lattice-Boltzmann simulation of two-phase flow in porous
media. Water Resour. Res. 40, 1–14 (2004).
218. Zhao, J. et al. Lattice Boltzmann simulation of liquid flow in nanoporous media. Int. J. Heat Mass
Transf. 125, 1131–1143 (2018).
219. Takaba, H., Onumata, Y. & Nakao, S. Molecular simulation of pressure-driven fluid flow in
nanoporous membranes. J. Chem. Phys. 127, 054703 (2007).
220. Adiga, S. P. & Brenner, D. W. Flow control through polymer-grafted smart nanofluidic channels:
Molecular dynamics simulations. Nano Lett. 5, 2509–2514 (2005).
221. Wang, L., Dumont, R. S. & Dickson, J. M. Nonequilibrium molecular dynamics simulation of
pressure-driven water transport through modified CNT membranes. J. Chem. Phys. 138, (2013).
222. Hoogerbrugge, P. J. & Koelman, J. M. V. A. Simulating Microscopic Hydrodynamic Phenomena
with Dissipative Particle Dynamics. Europhys. Lett. 19, 155–160 (1992).
223. Pan, W. & Tartakovsky, A. M. Dissipative particle dynamics model for colloid transport in porous
media. Adv. Water Resour. 58, 41–48 (2013).
224. Dzwinel, W. & Yuen, D. A. A Two-Level, Discrete-Particle Approach for Simulating Ordered
Colloidal Structures. J. Colloid Interface Sci. 225, 179–190 (2000).
225. Bolintineanu, D. S. et al. Particle dynamics modeling methods for colloid suspensions. Comput.
Part. Mech. 1, 321–356 (2014).
226. Abu-Nada, E. Modeling of Various Heat Transfer Problems Using Dissipative Particle Dynamics.
Numer. Heat Transf. Part A Appl. 58, 660–679 (2010).
227. Adiga, S. P. & Brenner, D. W. Stimuli-Responsive Polymer Brushes for Flow Control through
Nanopores. J. Funct. Biomater. 3, 239–256 (2012).
228. Kaçar, G. Dissipative particle dynamics simulation parameters and interactions of a hydrogel. J.
Turkish Chem. Soc. Sect. A Chem. 5, 19–28 (2018).
229. Warren, P. B. Vapor-liquid coexistence in many-body dissipative particle dynamics. Phys. Rev. E
68, 066702 (2003).
230. Liu, M., Meakin, P. & Huang, H. Dissipative particle dynamics simulation of pore-scale multiphase
fluid flow. Water Resour. Res. 43, 1–14 (2007).
231. Xia, Y. et al. A GPU-accelerated package for simulation of flow in nanoporous source rocks with
many-body dissipative particle dynamics. Comput. Phys. Commun. 247, 106874 (2020).
232. Revenga, M., Zúñiga, I. & Español, P. Boundary conditions in dissipative particle dynamics.
Comput. Phys. Commun. 121–122, 309–311 (1999).
233. Pivkin, I. V. & Karniadakis, G. E. A new method to impose no-slip boundary conditions in
dissipative particle dynamics. J. Comput. Phys. 207, 114–128 (2005).
234. Henrich, B., Cupelli, C., Moseler, M. & Santer, M. An adhesive DPD wall model for dynamic
wetting. Europhys. Lett. 80, 60004 (2007).
235. Li, Z., Bian, X., Tang, Y.-H. & Karniadakis, G. E. A dissipative particle dynamics method for
arbitrarily complex geometries. J. Comput. Phys. 355, 534–547 (2018).
236. Sun, C. et al. Mechanisms of Molecular Permeation through Nanoporous Graphene Membranes.
Langmuir 30, 675–682 (2014).

110

237. Hossain, J. Al & Kim, B. Scale Effect on Simple Liquid Transport through a Nanoporous Graphene
Membrane. Langmuir 37, 6498–6509 (2021).
238. Yuan, Z., Misra, R. P., Rajan, A. G., Strano, M. S. & Blankschtein, D. Analytical Prediction of Gas
Permeation through Graphene Nanopores of Varying Sizes: Understanding Transitions across
Multiple Transport Regimes. ACS Nano 13, 11809–11824 (2019).
239. Glavatskiy, K. S. & Bhatia, S. K. Effect of pore size on the interfacial resistance of a porous
membrane. J. Memb. Sci. 524, 738–745 (2017).
240. Liu, J., Li, S. & Chen, Y. A fast and practical method to pack spheres for mesh generation. Acta
Mech. Sin. 24, 439–447 (2008).
241. Liu, C. & Li, Z. Flow regimes and parameter dependence in nanochannel flows. Phys. Rev. E 80,
036302 (2009).
242. Pardo-Alonso, S., Vicente, J., Solórzano, E., Rodriguez-Perez, M. Á. & Lehmhus, D. Geometrical
Tortuosity 3D Calculations in Infiltrated Aluminium Cellular Materials. Procedia Mater. Sci. 4,
145–150 (2014).
243. Krause, B., Sijbesma, H. J. P., Münüklü, P., van der Vegt, N. F. A. & Wessling, M. Bicontinuous
Nanoporous Polymers by Carbon Dioxide Foaming. Macromolecules 34, 8792–8801 (2001).
244. Yang, X. et al. Preparation and lithium-storage performance of carbon/silica composite with a
unique porous bicontinuous nanostructure. Carbon N. Y. 77, 275–280 (2014).
245. Raeini, A. Q., Blunt, M. J. & Bijeljic, B. Direct simulations of two-phase flow on micro-CT images
of porous media and upscaling of pore-scale forces. Adv. Water Resour. 74, 116–126 (2014).
246. Zhao, J. et al. Study of Gas Flow Characteristics in Tight Porous Media with a Microscale Lattice
Boltzmann Model. Sci. Rep. 6, 32393 (2016).
247. Goral, J. et al. Confinement Effect on Porosity and Permeability of Shales. Sci. Rep. 10, 49 (2020).
248. Roy, S., Raju, R., Chuang, H. F., Cruden, B. A. & Meyyappan, M. Modeling gas flow through
microchannels and nanopores. J. Appl. Phys. 93, 4870–4879 (2003).
249. Pagonabarraga, I. & Frenkel, D. Dissipative particle dynamics for interacting systems. J. Chem.
Phys. 115, 5015–5026 (2001).
250. Chen-Wiegart, Y. C. K. et al. Tortuosity characterization of 3D microstructure at nano-scale for
energy storage and conversion materials. J. Power Sources 249, 349–356 (2014).
251. Rao, Q. et al. A modified many-body dissipative particle dynamics model for mesoscopic fluid
simulation: methodology, calibration, and application for hydrocarbon and water. Mol. Simul. 47,
363–375 (2021).
252. Müller-Plathe, F. Reversing the perturbation in nonequilibrium molecular dynamics: An easy way
to calculate the shear viscosity of fluids. Phys. Rev. E 59, 4894–4898 (1999).
253. Mosland, E. N., Lohne, K. D., Ystad, B. & Hallanger, A. Pressure Wave Velocity in Fluid-Filled
Pipes with and without Deposits in the Low-Frequency Range. J. Hydraul. Eng. 144, 04018064
(2018).
254. Sutera, S. P. & Skalak, R. The History of Poiseuille’s Law. Annu. Rev. Fluid Mech. 25, 1–20 (1993).
255. Koponen, A., Kataja, M. & Timonen, J. Permeability and effective porosity of porous media. Phys.
Rev. E 56, 3319–3325 (1997).
256. Kruczek, B. Encyclopedia of Membranes. Encyclopedia of Membranes (Springer Berlin Heidelberg,
2020). doi:10.1007/978-3-642-40872-4

Abstract (if available)

Abstract Nanoporous materials have been drawing increasing attention in the last decades due to their lightweight and large surface area to volume ratio. The materials have proved their usefulness in many applications, such as catalysts, sensors, and lightweight structural design. This research is focused on the mechanical and hydraulic properties of nanoporous materials, which are investigated using a combination of multi-scale simulation methods. First, we propose a direct method to generate three-dimensional (3D) stochastic bicontinuous nanoporous (BNP) microstructures on general periodic parallelepiped cells, providing flexibility to the design of nanoporous simulation models in different length scales. The method combines the ideas of Cahn with an adjustable Fibonacci grid to ensure the generation of a stochastic structure, while strictly enforcing periodicity in an arbitrary parallelepiped cell. The topological and morphological features, including ligament size, mean and Gaussian curvatures, and genus, are investigated to validate the method. The proposed method offers an efficient framework for the generation of 3D stochastic periodic bicontinuous microstructures. The method is particularly convenient in the design of adjustable nanoporous models for atomistic simulations. With the support of the novel method, we investigated the mechanical properties of BNP metallic glasses (BNPMGs), which were recently experimentally synthesized. We conduct molecular dynamics (MD) simulations of tensile loading deformation and failure of Cu(64)Zr(36) BNPMG with 55% porosity and 4.4 nm ligament size. Results indicate an anomalous mechanical behavior featuring delocalized plastic deformation preceding ductile failure. The deformation follows two mechanisms: i) Necking of ligaments aligned with the loading direction and ii) progressive alignment of randomly oriented ligaments. Failure occurs at 0.16 strain, followed by a massive rupture of ligaments. This work indicates that a BNP design is able to effectively delocalize strain localization in a metallic glass (MG) due to a combination of size effects on the ductility of MGs, resulting in nano ligaments necking and progressive asynchronous alignment of ligaments. We further investigate the effect of porosity levels and compositions on the mechanical behavior of such materials. MD simulations are employed to study the mechanical properties of nanoporous Cu(x)Zr(1-x) MGs with five different compositions, x = 0.28, 0.36, 0.5, 0.64, and 0.72, and porosity in the range from 0.1 to 0.7. Results from tensile loading simulations indicate a strong dependence of the Young’s modulus (E) and ultimate tensile strength (UTS) on porosity and composition. By increasing the porosity from 0.1 to 0.7, the topology of the nanoporous MG shifts from closed cell to open-cell bicontinuous. The change in nanoporous topology enables a brittle-to-ductile transition in deformation and failure mechanisms from a single critical shear band to necking and rupture of ligaments. Genetic Programming (GP) is employed to find scaling laws for E and UTS as a function of porosity and composition. A comparison of the GP-derived scaling laws against existing relationships shows that the GP method is able to uncover expressions that can predict accurately both the values of E and UTS in the whole range of porosity and compositions considered. Other than mechanical properties, a good understanding of the hydraulic behavior of nanoporous media is necessary for its wide application field, including electrokinetic energy conversion devices, water desalination using nanosized membranes, and extraction of hydrocarbons from ultra-tight shales. To address this, we employed Many-body Dissipative Particle Dynamics simulations to shed light on the fluid flow process through BNP media, which are representative models for a wide class of nanoporous materials. BNP models are constructed considering a defined morphology and porosity level and varying pore sizes. Fluid flow is generated through the nanoporous media by the action of two pistons added to the ends of the system. Simulation results, obtained by imposing different pressure differences along with the nanoporous model, indicate a linear pressure drop within the nanoporous structure, regardless of pore size. Regardless of the complexities and different scales of the porous media considered, the steady-state fluid flow through the nanoporous models is proportional to the pressure gradient applied, in agreement with Darcy’s law. The calculated pore-size dependence of permeability is well described by the Hagen-Poiseuille law, derived for single tubes, considering a single shape correction factor that accounts for the flow resistance due to the complex nanoporous morphology.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Atomistic modeling of the mechanical properties of metallic glasses, nanoglasses, and their composites

PDF

Synthesis, characterization, and mechanical properties of nanoporous foams

PDF

Large-scale molecular dynamics simulations of nano-structured materials

PDF

Phase change heterostructures for electronic and photonic applications

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Development and characterization of hierarchical cellular structures

PDF

Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras

PDF

Mechanical behavior and deformation of nanotwinned metallic alloys

PDF

Multi-scale quantum dynamics and machine learning for next generation energy applications

PDF

Mechanical behavior and microstructure optimization of ultrafine-grained aluminum alloys and nanocomposites

PDF

Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

PDF

Simulation and machine learning at exascale

PDF

Reactivation of multiple faults in oilfields with injection and production

$Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model$

PDF

Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model

PDF

Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow

PDF

Microstructural evolution and mechanical properties of a copper-zirconium alloy processed by severe plastic deformation

PDF

Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

PDF

Expanding the synthesis space of 3D nano- and micro-architected lattice materials

PDF

Exploring properties of silicon-carbide nanotubes and their composites with polymers

Asset Metadata

Creator Liu, Chang (author)

Core Title Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Materials Science

Degree Conferral Date 2022-08

Publication Date 07/22/2022

Defense Date 05/17/2022

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag bicontinuous,dissipative dynamics,genetic programing,Hagen-Poiseuille law,Materials Science,mechanical properties,metallic glass,molecular dynamics,nanoporous,OAI-PMH Harvest,permeability

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Branicio, Paulo (committee chair), Hodge, Andrea (committee member), Jha, Birendra (committee member), Nakano, Aiichiro (committee member), Vashishta, Priya (committee member)

Creator Email chang.liu8bcz745@gmail.com,liu396@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC111375213

Unique identifier UC111375213

Legacy Identifier etd-LiuChang-10928

Document Type Dissertation

Format application/pdf (imt)

Rights Liu, Chang

Type texts

Source 20220728-usctheses-batch-962 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu