- Email: [email protected]
Computationally efficient simulation method for conductivity modeling of 2D-based conductors
Computational Materials Science 161 (2019) 364–370
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Computationally efficient simulation method for conductivity modeling of 2D-based conductors
T
⁎
Leo Rizzia,c, , Andreas Zienerta,b, Jörg Schustera,b, Martin Köhnec, Stefan E. Schulza,b a
Faculty of Electrical Engineering and Information Technology, TU Chemnitz, Reichenhainer Str. 70, 09126 Chemnitz, Germany Fraunhofer ENAS, Technologie-Campus 3, 09126 Chemnitz, Germany c Robert Bosch GmbH, Robert-Bosch-Campus 1, 71272 Renningen, Germany b
A R T I C LE I N FO
A B S T R A C T
Keywords: 2D Nanocomposites Electrical conductivity Network simulation Nanostructured materials
Macroscopic materials made of two-dimensional components such as flakes of graphene or transition metal dichalcogenides represent a material class with great potential for large-scale applications. Depending on the structure, they can inherit the exceptional properties of the nanoscale building blocks while developing new features on the macroscopic scale. Supported by theoretical considerations and finite element analysis, we developed a network simulation method to model 2D-based electrical conductors. Here, we systematically explain the technical and methodological details of our approach, using the example of graphene-based conductor materials. Apart from the raw material properties, we discuss the importance of homogeneity and internal structure of the material. Our findings are supported by finite element analysis. We demonstrate the application of our method by studying the intricate interaction of several material parameters and the resulting effect on the macroscopic network. Finally, we provide guidelines for adapting our method to different physical situations.
1. Introduction The discovery of graphene in 2004 [1] and the Nobel Prize that followed in 2010 initiated intensive research on two-dimensional (2D) materials [2]. Apart from the elemental allotropes such as graphene or silicene [3], a variety of multi-atomic 2D compositions is readily available [4]: Examples include transition metal dichalcogenides [5] or hexagonal boron nitride [6]. However, among the 2D materials, graphene is still one of the most promising candidates for electronics. The atomically thin carbon layer exhibits extraordinary mechanical, electrical, and optical properties. Especially its high thermal [7] and electrical conductivity [8] as well as the enormous current carrying capacity [9] are ideal for conductor applications. Nevertheless, despite extensive research on graphene, the first large-scale application is yet to come. One approach to exploit the benefits of graphene is the synthesis of graphene-based conductor materials (GCMs) in the form of films [10,11] or fibers [12,13]. They consist of individual graphene flakes that are assembled in a layered structure to form a bulk material. The goal of this approach is to make the record electrical conductivity of graphene accessible on a macroscopic level while at the same time profiting from its high mechanical strength at a very low density. GCMs represent a material class that bridges the gap between the nano- and
⁎
the macroscale. The properties depend largely on the graphene flakes used in the production process. At the same time, new properties arise when the nanoscopic objects are combined and arranged to create a macroscopic compound. There is lots of experimental work on GCMs, but the lack of theoretical descriptions and simulations can impair technical progress. In a previous publication [14], we proposed a network simulation approach and used it to analyze the electrical characteristics of GCMs by tuning the properties of the nanoscale components. While we concentrated on the physical results in our previous publication, we are going to present the methodological details in this work. We provide a step-by-step explanation from the geometric modeling of a GCM to the computation of the overall electrical conductivity. We explain the advantages of our method and compare it to finite element analysis. As our approach is very general, our model is able to describe any macroscopic conductor composed of 2D flakes. 2. Computation method 2.1. Geometry generation At the lowest level, our simulation of a GCM starts with the generation of the input geometry. We create 2D layers filled with graphene
Corresponding author at: Robert Bosch GmbH, Robert-Bosch-Campus 1, 71272 Renningen, Germany. E-mail address: [email protected] (L. Rizzi).
https://doi.org/10.1016/j.commatsci.2019.02.022 Received 21 December 2018; Received in revised form 6 February 2019; Accepted 18 February 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Fig. 1. Geometry generation. (a) Randomly positioned points are generated in a 2D layer. (b) Points are connected and excess points are discarded as described in the main text. (c) Polygons are generated from connected points and shrunk to fit the desired packing density. Extremely thin polygons are divided in the middle as indicated by the red arrow.
Our geometry algorithm has several important advantages and offers control over various material parameters. The generated polygons are of random shapes and sizes and model real graphene flakes as closely as practically possible. The number N of points considered for each polygon defines the maximum number of vertices and thus, the complexity of the shapes. It can be adjusted as required. We keep control over the flake size distribution by manipulating the upper and lower flake size limits as well as the perimeters of the polygons. Another parameter that alters the average flake size is the number of total points per layer M . For a given N , more points per layer result in smaller polygons. The packing density and spatial flake distribution are accessible by choosing the appropriate downscaling of the polygons. With all these fine tuning possibilities, the algorithm allows us to consider boundary conditions set by production methods or the raw material itself while retaining a realistic degree of randomness. For our simulations, we model GCMs with a minimum of 20–30 layers and a total of 30.000–40.000 graphene flakes. At this size, statistical fluctuation is minimized while preserving fast computation times (see Section 3.1 for more details). A section of a finalized 3D GCM geometry is shown in Fig. 2. We implemented our method in Python 2.7 where we make extensive use of the Shapely library [17] to handle geometric objects.
Fig. 2. 3D section of a GCM. A finite thickness is assigned to the 2D layers and they are stacked to form a 3D object.
flakes that are represented by polygons. A finite thickness is assigned to the layers and they are stacked on top of each other to model a threedimensional (3D) object. The dimensions of the layers and their total number are determined by the desired dimensions of the 3D GCM. The geometry algorithm is illustrated in Fig. 1: As a first step, each layer is filled with M randomly positioned points (Fig. 1a). Subsequently, the individual points are selected one by one in random order. At the first selected point P , we choose the N nearest neighboring points and form the convex hull. Thus, a polygon is created. From the remaining points, a new one is selected to form a new polygon. The procedure is repeated until all points are considered either as selected points or as neighboring points. Due to the randomness of the process, some points might not be incorporated into polygons, and voids might be created within a layer. Our algorithm identifies and removes leftover points, and fills the voids with new polygons so that we obtain a completely filled layer (Fig. 1b). At this stage, we refine the geometry to fit realistic boundary conditions: During the production of highly conductive GCMs, extremely small graphene flakes are often undesirable and are thus removed [15,16]. On the other hand, even the largest graphene flakes are limited in size due to the flake dimensions of the raw material and due to production methods (e.g. sonication is known to damage graphene flakes). To account for these aspects, our algorithm removes all the flakes with an area below a certain threshold and splits polygons that exceed a certain upper limit. We avoid extremely long and thin shapes by also restricting the perimeters of the polygons. These measures slightly reduce the packing density p in each layer but it remains close to 1. To finalize the geometry, we separate the polygons by scaling them down (Fig. 1c). We have two main options to adjust the overall packing density. We can either continue to scale down the polygons until the desired packing density is obtained or we eliminate randomly selected flakes from each layer. The first approach results in a homogeneous spatial distribution of the polygons while the second one corresponds to the simulation of spatially clustered graphene flakes. We are interested in densely packed GCMs (p > 0.70) as they are the best electrical conductors. Compared to foams or other porous structures, the greater amount of current-carrying material leads to higher conductivities. Because of the high packing density, we can model graphene flakes as rigid. It can be assumed that the flakes sufficiently support each other to prevent sagging.
2.2. Resistor network generation Once the 3D geometry is completed, we transform it into a graph with nodes and edges. We iterate through the individual layers of the 3D object and identify the overlaps between flakes of adjacent layers. Edge-to-edge contact within a single layer is not considered as we expect the contribution to current flow to be negligible. For a substantial contribution, a strong edge-to-edge connection of a large number of flakes per layer would be required. Due to variable flake shapes with random orientation, such an arrangement is highly improbable, as it would correspond to a spontaneous arrangement with ideal packing. Overlaps between adjacent layers on the other hand are inevitable. They cover large areas of the graphene flakes, so that many more atoms than at the flake edges can participate in electronic transport. Hence, tunneling from flake to flake is facilitated. Since the overlapping surfaces are crucial for the formation of current paths, we choose them as the nodes for our graph. At the geometric center of every overlap, a node is created in each of the two flakes involved. Nodes that share the same graphene flake or the same overlap are subsequently connected by an edge. In this way, a 3D graph is obtained (Fig. 3). As the next step, we assign electrical conductances to the edges based on geometric specifications. At the overlaps, we compute the conductance Gout from one flake to another according to
Gout = σout·
Aout l out
(1)
with the electrical conductivity σout between two graphene flakes, also referred to as out-of-plane conductivity, the distance l out between adjacent layers and the area of the overlap Aout . The situation is more complex for the in-plane conductances Gin since many nodes can be 365
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Fig. 3. Graph construction. (a) Top view of two overlapping graphene layers. (b) Central section from (a). Overlaps are identified; nodes are generated and connected. (c) Side view of four connected layers.
interconnected within a single flake. Current is not limited to the flake area directly between overlaps but it can potentially spread through the whole graphene flake. Thus, overlaps are not point-like nodes but objects with a finite surface area. Due to multiple possible overlaps, an accurate description of current flow would require additional methods such as finite element calculations. For large structures, this would result in unacceptable computation times. Nevertheless, the conductance values Gin should be dependent on geometric specifications, namely the total flake area, the distance between connected overlaps, the areas, and the orientations of the overlaps. Current flow from or to large overlaps is expected to be more efficient due to a smaller spreading resistance. To solve this dilemma, we chose to combine several simulation methods. We set up our networks by considering an empirical parameter that is determined by finite element modeling of medium-sized structures. Note that the following considerations are necessary when transitioning from 1D-based to 2D-based materials. In nanowire networks or CNT-based macro-materials for example, inplane conductances can be modeled as 1D, which facilitates the network setup significantly [18,19]. As a starting value for the in-plane conductivity between two overlaps, we compute the electrical conductance of a trapezoid as depicted in Fig. 4. The height of the trapezoid corresponds to the connection line between the node positions, i.e. the centers of the overlaps. The base edges are the intersections of the overlaps and the perpendiculars to the connection line at the two nodes. With the base edges − being the two terminals, the conductance Gin of such a trapezoidal resistor is given as
The trapezoidal description considers several important geometrical aspects: The distance between the overlaps, the part of the flake where the majority of the current is expected to flow, the area of the overlaps and their orientation towards each other. The conductances are normalized for every individual flake such that the total current within a given flake does not exceed the theoretical maximum −
Gin = Gin
tgr Δb ⎛ ⎛ b1 ⎞ ⎞ ⎜log ⎟ l in ⎝ ⎝ b2 ⎠ ⎠ ⎜
⎟
(2)
with the in-plane electrical conductivity σin , the thickness of a graphene sheet tgr , the distance between the two nodes, i.e. the height of the trapezoid l in , the lengths of the baselines b1 and b2 (b1 > b2 ) and their difference Δb = b1 − b2 > 0 . For b1 = b2 = b, Eq. (2) becomes −
Gin = σin·
(4)
Gin is the normalized conductance between two overlaps, Aflake signifies the area of the flake, and the Ai correspond to the areas of the individual trapezoids in the given flake. A normalization is necessary because for a large number of overlaps in a single flake, the sum of the areas of the trapezoids ∑i Ai can be much greater than Aflake . Without normalization to Aflake , the conductances could become unphysically large and allow more current flow than what is physically achievable. For a small number of overlaps on the other hand, the normalization to Aflake area guarantees that current is not limited to the region between two overlaps. During normalization, each trapezoid is equally weighted, irrespective of the amount of current actually flowing through it. This introduces a source of error that needs to be taken into account. An alternative approach that considers current flow requires a self-consistency loop where each connection is re-weighted after the electrical network has been set up and solved. However, we find that this procedure leads to substantially longer computation times and still unsatisfactory precision, which is why we opt for a different method: We introduce an empirical parameter k to compensate for the approximations in our model. We initially determine this parameter by generating test structures and comparing multiple methods. The computation of k only has to be performed once and will be explained comprehensively in the following section. Before that, we demonstrate how the electrical network is constructed and how it is solved: We set up the conductance matrix (Gij ) of the electrical network according to nodal analysis. On the main diagonal, for every node i , we write the sum of all conductances of the edges that are connected to i
−1
−
Gin = σin·
Aflake ∑i Ai
tgr b l in
Gii =
(3)
∑j Gij,
i≠j
(5)
The off-diagonal elements of the conductance matrix are the negative conductances of the edges between the nodes i and j : −Gij ∀ i ≠ j . If nodes i and j are not connected, Gij is zero. To compute the total electrical conductivity of the GCM, we simulate a current through the structure. We define a border region at the edges of the structure where we connect all flakes to current sources (Fig. 5). At the points where current is injected, the current density is kept constant. This corresponds to homogeneous current injection at the side faces of the GCM. Furthermore, current is conserved, i.e. the same absolute amount of current Itot that enters one side leaves the structure on the other side. Apart from the border regions, the external current at every node is zero, in accordance with Kirchhoff’s current law. For a total of n nodes, we obtain the matrix form of the current-
Fig. 4. Trapezoid used for the calculation of the in-plane conductance. 366
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
The thickness h is given as the total number of 2D layers times the average layer thickness. The width w corresponds to the vertical axis of the sample in Fig. 5. 2.3. Adjusting the empirical parameter k As shown in our previous publication [14], the total conductivity σtot of an idealized GCM is ultimately limited by the packing density p of the GCM and the in-plane conductivity σin of the individual graphene flakes according to
σtot ≤ p ·σin ≡ σmax
The empirical parameter k has to be chosen such that this fundamental upper limit is respected. However, since the limit value p ·σin was derived for an idealized GCM, we need to investigate σmax for structures that are more realistic. We generate geometries as described above with varying packing densities in the relevant range for GCMs p > 0.70 . Subsequently, we compute the limit value of the total connetwork . From the same geometries, ductivity within our network model σmax we randomly select sections where we perform finite element analysis FEA (FEA) and compute the limit value of the total conductivity σmax . For FEA, we use the commercial software COMSOL Multiphysics® [20]. FEA is computationally far more expensive than our network approach, which is why we have to limit FEA to exemplary sections of the complete geometries. By considering multiple sections per structure (Fig. 6) however, we minimize the effects of statistical fluctuation. The results are compared to the limit value of the full structure computed with our network network approach σmax . FEA of various structures investigated with The conductivity limits σmax FEA are shown in Fig. 7a. As expected from theoretical prediction, the packing density and the in-plane conductivity are decisive factors for the limit of the total electrical conductivity. Apart from that, the homogeneity of the spatial distribution of graphene flakes plays a major role. Most notably, our FEA shows that for a GCM with homogeneously distributed graphene flakes, the limit value of the electrical conductivity corresponds to the limit value of the idealized GCM. This holds true for arbitrarily shaped and oriented flakes. Correspondingly, network we obtain the largest values σmax within our network model for homogeneous structures. Our FEA shows that inhomogeneities such as holes and clusters within a layer substantially reduce the total electrical conductivity, even if the overall packing density and the in-plane conductivity of the graphene flakes remain the same. From a physical point of view, these findings are coherent. Inhomogeneities lead to constrictions in the current path and therefore longer effective pathways, which in turn cause a smaller electrical conductivity. As a prime example of an inhomogeneous spatial flake distribution, we model a checkerboard structure (Inset Fig. 7). Note that this geometry is regular but it is not homogeneous in terms of spatial distribution. In a homogeneous
Fig. 5. Top view of a GCM. Current enters or leaves the structure through the yellow flakes at the sides. All the flakes that intersect with the red lines are considered for voltage measurement.
voltage-equation
⎛ G11 −G12 ⋯ −G1n ⎞ ⎛ Φ1 ⎞ ⎛ I1 ⎞ ⎜ −G21 G22 ⋯ −G2n ⎟·⎜ Φ2 ⎟ = ⎜ I2 ⎟ ⋯ ⋯ ⋯ ⎟ ⎜ ⋯ ⎟ ⎜⋯⎟ ⎜⎜ ⋯ ⎟⎜ ⎟ ⎜ ⎟ ⎝ −Gn1 −Gn2 ⋯ Gnn ⎠ ⎝ Φn ⎠ ⎝ In ⎠
(6)
with the conductances Gij , the electric potential Φi at node i and the injected current values Ii at node i . The conductance matrix is singular and we are free to choose a reference node where the potential is set to zero. We subsequently solve the equation for the potential vector so that we obtain the electric potential at every overlap in the structure. To compute the global voltage drop, we define two probe surfaces parallel to the side faces of the GCM where current is injected (Fig. 5). They are located outside the border region and separated by a distance d that is at least five times as large as the largest graphene flakes. The probe voltage is computed by averaging all the potentials at the overlaps that intersect with the respective probe surface. The potential difference V between the two probes corresponds to the voltage drop measured in a four-point probes configuration. For GCM simulations with tens of thousands of graphene flakes and sufficiently large d , the statistical error of this approach is negligible. We chose this form of measurement, as we are interested in the pure material parameters and not in contact resistances or voltage probe effects. However, our code could be used to simulate surface voltage probes of arbitrary sizes as well. With the empirical parameter k , the probe distance d , the thickness h and the width w of the GCM, the total electrical conductivity σtot can be determined according to
σtot = k ·
Itot hw · V d
(8)
(7)
Fig. 6. The central structure is generated by our geometry algorithm and evaluated with our network model. Randomly chosen sections are investigated with FEA. The color code indicates the homogeneous voltage drop across the sections. 367
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Fig. 7. (a) Limit of the total electrical conductivity versus overall packing density. The symbols mark different degrees of spatial homogeneity. The green line is the theoretical maximum. The red line is the linear fit for a checkerboard structure, shown as the inset. (b) Empirical parameter k versus packing density.
structure, the probability to find a flake at a certain position is constant across the whole 3D object. In a checkerboard structure on the other hand, the regular arrangement of square flakes leads to periodic constrictions across all layers. This reduces the overall electrical conductivity significantly. From our FEA results, we determine the empirical parameter k : As reference value, we choose the limit of the electrical conductivity of extremely homogeneous structures. For those FEA cases, σmax is equal to the theoretical maximum p ·σin . We choose k such that the same holds true for the network model and we obtain FEA network σmax = k ·σmax = p ·σin . The value of k is dependent on the overall packing density but does not vary with flake size or other microscopic parameters within the boundaries described in the geometry section. It ranges from 2.25 at p = 0.70 to 2.12 at p = 0.99. The full dependence is shown in Fig. 7b. Once the empirical parameter is adjusted, computations such as the ones for our previously published results can be performed. No further comparisons between methods are required and the high computational speed of our network simulation can be exploited. Fig. 8. Total conductivity of a GCM as a function of out-of-plane conductivity. Three characteristic regimes are indicated.
3. Discussion 3.1. Convergence and robustness of the method
plane conductivity. For this exemplary study, we investigated a system of 0.8 packing density, 25 layers and 39,000 flakes in total. Three regimes can be identified: For small σout , the overall conductivity is suppressed (Regime I), for very large σout , saturation is achieved (Regime III). The largest effects take place in Regime II, the transitional regime. Over a comparatively small range of out-of-plane conductivity values, the total conductivity moves from negligibly small to the saturation value. In Eq. (8), we already introduced the saturation value as the upper limit for the total conductivity of a GCM σmax . It is defined as the product of the packing density p and the in-plane conductivity σin . However, this simple mathematical dependence only holds true for highly conductive out-of-plane connections – a state that is achieved by large flake overlaps and high out-of-plane conductivities. For smaller out-of-plane conductivities, the impact of both p and σin decreases. This effect is visualized in Fig. 9. Fig. 9a shows the total conductivity of a GCM with respect to the in-plane conductivity of the individual graphene flakes. The different curves correspond to different out-of-plane conductivities σout . In general, a larger in-plane conductivity is always beneficial for the total conductivity of the macro-material. But the smaller σout , the more the overall conductivity levels off, even for the largest values of σin. The physical interpretation is straightforward: The connections between individual flakes become less and less conductive. At some point, even infinitely large in-plane conductivities cannot compensate small out-of-plane conductivities. The four upper curves in Fig. 9a are located in the transition regime as indicated in Fig. 8. Here, a
As our simulation method relies on random geometric elements, a minimum number of flakes has to be included in the simulation to obtain statistically reliable results. Furthermore, to treat a simulated GCM as one material and apply Eq. (7), the dimensions of the modeled GCM have to exceed a certain threshold. In extensive parameter studies, we analyzed the necessary conditions for successful modeling with our network method. In brief, a minimum number of 20–30 layers is required for convergence, depending on the structural details. This corresponds to a thickness of 6–10 nm, which is easily fulfilled by real GCMs. The length and width of the simulated GCM should be at least 25 times as large as the average flake length. For statistical robustness, a total of 30.000 flakes across all layers is required. The most important results of the convergence tests are illustrated in more detail in the Supplementary Material. If extremely inhomogeneous materials are simulated, the system size should be 1.5 times as large. With these parameters, the statistical fluctuation of the resulting total electrical conductivity can be limited to less than 2%. 3.2. Application of the method In this section, we use our method to study the interplay between several material parameters in a GCM. We focus on a critical regime in the parameter space that we did not analyze in our previous work. In Fig. 8, the total conductivity of a GCM is plotted as a function of out-of368
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Fig. 9. Total conductivity of a GCM as a function of (a) in-plane conductivity and (b) packing density for various σout . The symbols represent simulated data points while the lines indicate the trend. For (a), p = 0.8 and for (b), σin = 100 MS/m .
small increment in out-of-plane conductivity has a strong impact on the total conductivity. Fig. 9b shows the total conductivity of a GCM as a function of packing density for different values of σout . Higher packing density always corresponds to more current-carrying material and thus a higher total conductivity. But a minimum out-of-plane conductivity is required in order to exploit the material to full extent. Analogous to the in-plane-conductivity, the influence of the packing density decreases for smaller σout . Mathematically, this can be deduced from the slope of the curves, which decreases continuously with decreasing σout . We conclude that when fabricating a high-performance GCM, the material parameters have to be adjusted carefully. As different influencing factors compete, they overlap or even dominate each other. Experimental results are thus difficult to interpret and our simulation helps to understand the underlying physical reasons. Further comprehensive parameter studies and a comparison to experimental data can be found in our previous publication [14].
In conclusion, we present the methodological details of a modeling approach for the electrical conductivity of 2D-based conductors. We generate randomized geometries where we consider a variety of microscopic material properties such as in-plane and out-of-plane conductivity, packing density, or flake size. These parameters can be adjusted to fit the desired materials and the experimental conditions. Starting from the geometric structure of a material, we create a 3D resistor network that allows us to compute the overall electrical conductivity. By varying the microscopic parameters, their influence on the macroscopic electrical conductivity can be analyzed. We demonstrate this by investigating the complex interplay of packing density, in-plane and out-of-plane conductivity in a GCM. The method is fast and efficient compared to other simulative approaches such that large systems with several ten thousand flakes can be investigated.
3.3. Conclusions and outlook
Leo Rizzi: Conceptualization, Methodology, Software, Validation, Investigation, Writing - original draft. Andreas Zienert: Methodology, Writing - review & editing. Jörg Schuster: Conceptualization, Methodology, Resources, Writing - review & editing. Martin Köhne: Conceptualization, Project administration. Stefan E. Schulz: Conceptualization, Resources, Supervision.
CRediT authorship contribution statement
The strength and speed of our modeling method for 2D-based conductors come from its simplicity. By using simple geometric objects and a variable parameter set, we can adjust the computation to various physical situations. For pure graphene-only GCMs, the modeling is performed as described above. The conductivities should be varied according to statistical distributions across the whole GCM. Multilayers and defects are considered by increasing the thickness and adjusting the conductivity of the respective flakes. In principle, any geometry can be used as a starting point. To include dopants, the polygonal objects have to be regarded as graphene platelets decorated with non-carbon atoms or molecules. In-plane and out-of-plane conductivities need to be modified as well as the thickness of the platelets. Again, statistical distributions should be applied to decide which flakes are decorated with dopants and which are not. In our previous publication, we used this procedure to compare experimental data to our simulations and we obtained excellent agreement [14]. Going beyond graphene-based materials, the method can be used to model any bulk material composed of quasi-2D elements or similar ultrathin components. Examples might include conductive polymer nanocomposites or composites based on transition metal dichalcogenides. In our study, we concentrated on GCMs since they appear to be the most promising materials for conductor applications as of right now. To the best of our knowledge, no other simulation method is able to model similar systems more efficiently. Compared to FEA, our approach requires less time and is less memory intensive at system sizes that are ten times as large or more. Due to that, our approach enables the usage of statistical methods to reduce the fluctuation of the results.
Acknowledgment The authors would like to thank W. Pieper and F. Teichert for fruitful discussions as well as E. Lorenz for technical support. Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. Declaration of interest None. Data availability The raw data required to reproduce these findings cannot be shared at this time due to technical limitations. The processed data required to reproduce these findings cannot be shared at this time due to time limitations. However, all the raw input parameters required to reproduce our results are presented in this article. The processed data is fully presented in the form of graphs and plots. 369
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Appendix A. Supplementary material
Phys. Lett. 91 (2007) 163513. [10] S. Stankovich, D.A. Dikin, G.H.B. Dommett, K.M. Kohlhaas, E.J. Zimney, E.A. Stach, R.D. Piner, S.T. Nguyen, R.S. Ruoff, Graphene-based composite materials, Nature 442 (2006) 282. [11] S. Stankovich, D.A. Dikin, R.D. Piner, K.A. Kohlhaas, A. Kleinhammes, Y. Jia, Y. Wu, S.T. Nguyen, R.S. Ruoff, Synthesis of graphene-based nanosheets via chemical reduction of exfoliated graphite oxide, Carbon 45 (2007) 1558–1565. [12] Z. Xu, C. Gao, Graphene chiral liquid crystals and macroscopic assembled fibres, Nat. Commun. 2 (2011) 571. [13] Z. Xu, C. Gao, Graphene fiber: a new trend in carbon fibers, Mater. Today 18 (2015) 480–492. [14] L. Rizzi, A. Zienert, J. Schuster, M. Köhne, S.E. Schulz, Electrical conductivity modeling of graphene-based conductor materials, ACS Appl. Mater. Interfaces 10 (2018) 43088–43094. [15] L. Peng, Z. Xu, Z. Liu, Y. Guo, P. Li, C. Gao, Ultrahigh thermal conductive yet superflexible graphene films, Adv. Mater. 29 (2017) 1700589. [16] Z. Wang, B. Mao, Q. Wang, J. Yu, J. Dai, R. Song, Z. Pu, D. He, Z. Wu, S. Mu, Ultrahigh conductive copper/large flake size graphene heterostructure thin-film with remarkable electromagnetic interference shielding effectiveness, Small 14 (2018) 1704332. [17] S. Gillies et al., Shapely: manipulation and analysis of geometric objects, 2007. [18] C. O'Callaghan, C. Gomes da Rocha, H.G. Manning, J.J. Boland, M.S. Ferreira, Effective medium theory for the conductivity of disordered metallic nanowire networks, Phys. Chem. Chem. Phys. 18 (2016) 27564–27571. [19] F. Dalmas, R. Dendievel, L. Chazeau, J.-Y Cavaillé, C. Gauthier, Carbon nanotubefilled polymer composites. Numerical simulation of electrical conductivity in threedimensional entangled fibrous networks, Acta Mater. 54 (2006) 2923–2931. [20] COMSOL Multiphysics® v. 5.3, COMSOL AB, Stockholm, Sweden. www.comsol. com.
Supplementary data to this article can be found online at https:// doi.org/10.1016/j.commatsci.2019.02.022. References [1] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Electric field effect in atomically thin carbon films, Science 306 (2004) 666–669. [2] K.S. Novoselov, V.I. Fal′ko, L. Colombo, P.R. Gellert, M.G. Schwab, K. Kim, A roadmap for graphene, Nature 490 (2012) 192. [3] P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M.C. Asensio, A. Resta, B. Ealet, G. Le Lay, Silicene: compelling experimental evidence for graphene-like two-dimensional silicon, Phys. Rev. Lett. 108 (2012) 155501. [4] P. Miró, M. Audiffred, T. Heine, An atlas of two-dimensional materials, Chem. Soc. Rev. 43 (2014) 6537–6554. [5] Q.H. Wang, K. Kalantar-Zadeh, A. Kis, J.N. Coleman, M.S. Strano, Electronics and optoelectronics of two-dimensional transition metal dichalcogenides, Nat. Nanotechnol. 7 (2012) 699. [6] L. Song, L. Ci, H. Lu, P.B. Sorokin, C. Jin, J. Ni, A.G. Kvashnin, D.G. Kvashnin, J. Lou, B.I. Yakobson, P.M. Ajayan, Large scale growth and characterization of atomic hexagonal boron nitride layers, Nano Lett. 10 (2010) 3209–3215. [7] A.A. Balandin, Thermal properties of graphene and nanostructured carbon materials, Nat. Mater. 10 (2011) 569. [8] J.-H. Chen, C. Jang, S. Xiao, M. Ishigami, M.S. Fuhrer, Intrinsic and extrinsic performance limits of graphene devices on SiO2, Nat. Nanotechnol. 3 (2008) 206. [9] J. Moser, A. Barreiro, A. Bachtold, Current-induced cleaning of graphene, Appl.
370
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Computationally efficient simulation method for conductivity modeling of 2D-based conductors
T
⁎
Leo Rizzia,c, , Andreas Zienerta,b, Jörg Schustera,b, Martin Köhnec, Stefan E. Schulza,b a
Faculty of Electrical Engineering and Information Technology, TU Chemnitz, Reichenhainer Str. 70, 09126 Chemnitz, Germany Fraunhofer ENAS, Technologie-Campus 3, 09126 Chemnitz, Germany c Robert Bosch GmbH, Robert-Bosch-Campus 1, 71272 Renningen, Germany b
A R T I C LE I N FO
A B S T R A C T
Keywords: 2D Nanocomposites Electrical conductivity Network simulation Nanostructured materials
Macroscopic materials made of two-dimensional components such as flakes of graphene or transition metal dichalcogenides represent a material class with great potential for large-scale applications. Depending on the structure, they can inherit the exceptional properties of the nanoscale building blocks while developing new features on the macroscopic scale. Supported by theoretical considerations and finite element analysis, we developed a network simulation method to model 2D-based electrical conductors. Here, we systematically explain the technical and methodological details of our approach, using the example of graphene-based conductor materials. Apart from the raw material properties, we discuss the importance of homogeneity and internal structure of the material. Our findings are supported by finite element analysis. We demonstrate the application of our method by studying the intricate interaction of several material parameters and the resulting effect on the macroscopic network. Finally, we provide guidelines for adapting our method to different physical situations.
1. Introduction The discovery of graphene in 2004 [1] and the Nobel Prize that followed in 2010 initiated intensive research on two-dimensional (2D) materials [2]. Apart from the elemental allotropes such as graphene or silicene [3], a variety of multi-atomic 2D compositions is readily available [4]: Examples include transition metal dichalcogenides [5] or hexagonal boron nitride [6]. However, among the 2D materials, graphene is still one of the most promising candidates for electronics. The atomically thin carbon layer exhibits extraordinary mechanical, electrical, and optical properties. Especially its high thermal [7] and electrical conductivity [8] as well as the enormous current carrying capacity [9] are ideal for conductor applications. Nevertheless, despite extensive research on graphene, the first large-scale application is yet to come. One approach to exploit the benefits of graphene is the synthesis of graphene-based conductor materials (GCMs) in the form of films [10,11] or fibers [12,13]. They consist of individual graphene flakes that are assembled in a layered structure to form a bulk material. The goal of this approach is to make the record electrical conductivity of graphene accessible on a macroscopic level while at the same time profiting from its high mechanical strength at a very low density. GCMs represent a material class that bridges the gap between the nano- and
⁎
the macroscale. The properties depend largely on the graphene flakes used in the production process. At the same time, new properties arise when the nanoscopic objects are combined and arranged to create a macroscopic compound. There is lots of experimental work on GCMs, but the lack of theoretical descriptions and simulations can impair technical progress. In a previous publication [14], we proposed a network simulation approach and used it to analyze the electrical characteristics of GCMs by tuning the properties of the nanoscale components. While we concentrated on the physical results in our previous publication, we are going to present the methodological details in this work. We provide a step-by-step explanation from the geometric modeling of a GCM to the computation of the overall electrical conductivity. We explain the advantages of our method and compare it to finite element analysis. As our approach is very general, our model is able to describe any macroscopic conductor composed of 2D flakes. 2. Computation method 2.1. Geometry generation At the lowest level, our simulation of a GCM starts with the generation of the input geometry. We create 2D layers filled with graphene
Corresponding author at: Robert Bosch GmbH, Robert-Bosch-Campus 1, 71272 Renningen, Germany. E-mail address: [email protected] (L. Rizzi).
https://doi.org/10.1016/j.commatsci.2019.02.022 Received 21 December 2018; Received in revised form 6 February 2019; Accepted 18 February 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Fig. 1. Geometry generation. (a) Randomly positioned points are generated in a 2D layer. (b) Points are connected and excess points are discarded as described in the main text. (c) Polygons are generated from connected points and shrunk to fit the desired packing density. Extremely thin polygons are divided in the middle as indicated by the red arrow.
Our geometry algorithm has several important advantages and offers control over various material parameters. The generated polygons are of random shapes and sizes and model real graphene flakes as closely as practically possible. The number N of points considered for each polygon defines the maximum number of vertices and thus, the complexity of the shapes. It can be adjusted as required. We keep control over the flake size distribution by manipulating the upper and lower flake size limits as well as the perimeters of the polygons. Another parameter that alters the average flake size is the number of total points per layer M . For a given N , more points per layer result in smaller polygons. The packing density and spatial flake distribution are accessible by choosing the appropriate downscaling of the polygons. With all these fine tuning possibilities, the algorithm allows us to consider boundary conditions set by production methods or the raw material itself while retaining a realistic degree of randomness. For our simulations, we model GCMs with a minimum of 20–30 layers and a total of 30.000–40.000 graphene flakes. At this size, statistical fluctuation is minimized while preserving fast computation times (see Section 3.1 for more details). A section of a finalized 3D GCM geometry is shown in Fig. 2. We implemented our method in Python 2.7 where we make extensive use of the Shapely library [17] to handle geometric objects.
Fig. 2. 3D section of a GCM. A finite thickness is assigned to the 2D layers and they are stacked to form a 3D object.
flakes that are represented by polygons. A finite thickness is assigned to the layers and they are stacked on top of each other to model a threedimensional (3D) object. The dimensions of the layers and their total number are determined by the desired dimensions of the 3D GCM. The geometry algorithm is illustrated in Fig. 1: As a first step, each layer is filled with M randomly positioned points (Fig. 1a). Subsequently, the individual points are selected one by one in random order. At the first selected point P , we choose the N nearest neighboring points and form the convex hull. Thus, a polygon is created. From the remaining points, a new one is selected to form a new polygon. The procedure is repeated until all points are considered either as selected points or as neighboring points. Due to the randomness of the process, some points might not be incorporated into polygons, and voids might be created within a layer. Our algorithm identifies and removes leftover points, and fills the voids with new polygons so that we obtain a completely filled layer (Fig. 1b). At this stage, we refine the geometry to fit realistic boundary conditions: During the production of highly conductive GCMs, extremely small graphene flakes are often undesirable and are thus removed [15,16]. On the other hand, even the largest graphene flakes are limited in size due to the flake dimensions of the raw material and due to production methods (e.g. sonication is known to damage graphene flakes). To account for these aspects, our algorithm removes all the flakes with an area below a certain threshold and splits polygons that exceed a certain upper limit. We avoid extremely long and thin shapes by also restricting the perimeters of the polygons. These measures slightly reduce the packing density p in each layer but it remains close to 1. To finalize the geometry, we separate the polygons by scaling them down (Fig. 1c). We have two main options to adjust the overall packing density. We can either continue to scale down the polygons until the desired packing density is obtained or we eliminate randomly selected flakes from each layer. The first approach results in a homogeneous spatial distribution of the polygons while the second one corresponds to the simulation of spatially clustered graphene flakes. We are interested in densely packed GCMs (p > 0.70) as they are the best electrical conductors. Compared to foams or other porous structures, the greater amount of current-carrying material leads to higher conductivities. Because of the high packing density, we can model graphene flakes as rigid. It can be assumed that the flakes sufficiently support each other to prevent sagging.
2.2. Resistor network generation Once the 3D geometry is completed, we transform it into a graph with nodes and edges. We iterate through the individual layers of the 3D object and identify the overlaps between flakes of adjacent layers. Edge-to-edge contact within a single layer is not considered as we expect the contribution to current flow to be negligible. For a substantial contribution, a strong edge-to-edge connection of a large number of flakes per layer would be required. Due to variable flake shapes with random orientation, such an arrangement is highly improbable, as it would correspond to a spontaneous arrangement with ideal packing. Overlaps between adjacent layers on the other hand are inevitable. They cover large areas of the graphene flakes, so that many more atoms than at the flake edges can participate in electronic transport. Hence, tunneling from flake to flake is facilitated. Since the overlapping surfaces are crucial for the formation of current paths, we choose them as the nodes for our graph. At the geometric center of every overlap, a node is created in each of the two flakes involved. Nodes that share the same graphene flake or the same overlap are subsequently connected by an edge. In this way, a 3D graph is obtained (Fig. 3). As the next step, we assign electrical conductances to the edges based on geometric specifications. At the overlaps, we compute the conductance Gout from one flake to another according to
Gout = σout·
Aout l out
(1)
with the electrical conductivity σout between two graphene flakes, also referred to as out-of-plane conductivity, the distance l out between adjacent layers and the area of the overlap Aout . The situation is more complex for the in-plane conductances Gin since many nodes can be 365
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Fig. 3. Graph construction. (a) Top view of two overlapping graphene layers. (b) Central section from (a). Overlaps are identified; nodes are generated and connected. (c) Side view of four connected layers.
interconnected within a single flake. Current is not limited to the flake area directly between overlaps but it can potentially spread through the whole graphene flake. Thus, overlaps are not point-like nodes but objects with a finite surface area. Due to multiple possible overlaps, an accurate description of current flow would require additional methods such as finite element calculations. For large structures, this would result in unacceptable computation times. Nevertheless, the conductance values Gin should be dependent on geometric specifications, namely the total flake area, the distance between connected overlaps, the areas, and the orientations of the overlaps. Current flow from or to large overlaps is expected to be more efficient due to a smaller spreading resistance. To solve this dilemma, we chose to combine several simulation methods. We set up our networks by considering an empirical parameter that is determined by finite element modeling of medium-sized structures. Note that the following considerations are necessary when transitioning from 1D-based to 2D-based materials. In nanowire networks or CNT-based macro-materials for example, inplane conductances can be modeled as 1D, which facilitates the network setup significantly [18,19]. As a starting value for the in-plane conductivity between two overlaps, we compute the electrical conductance of a trapezoid as depicted in Fig. 4. The height of the trapezoid corresponds to the connection line between the node positions, i.e. the centers of the overlaps. The base edges are the intersections of the overlaps and the perpendiculars to the connection line at the two nodes. With the base edges − being the two terminals, the conductance Gin of such a trapezoidal resistor is given as
The trapezoidal description considers several important geometrical aspects: The distance between the overlaps, the part of the flake where the majority of the current is expected to flow, the area of the overlaps and their orientation towards each other. The conductances are normalized for every individual flake such that the total current within a given flake does not exceed the theoretical maximum −
Gin = Gin
tgr Δb ⎛ ⎛ b1 ⎞ ⎞ ⎜log ⎟ l in ⎝ ⎝ b2 ⎠ ⎠ ⎜
⎟
(2)
with the in-plane electrical conductivity σin , the thickness of a graphene sheet tgr , the distance between the two nodes, i.e. the height of the trapezoid l in , the lengths of the baselines b1 and b2 (b1 > b2 ) and their difference Δb = b1 − b2 > 0 . For b1 = b2 = b, Eq. (2) becomes −
Gin = σin·
(4)
Gin is the normalized conductance between two overlaps, Aflake signifies the area of the flake, and the Ai correspond to the areas of the individual trapezoids in the given flake. A normalization is necessary because for a large number of overlaps in a single flake, the sum of the areas of the trapezoids ∑i Ai can be much greater than Aflake . Without normalization to Aflake , the conductances could become unphysically large and allow more current flow than what is physically achievable. For a small number of overlaps on the other hand, the normalization to Aflake area guarantees that current is not limited to the region between two overlaps. During normalization, each trapezoid is equally weighted, irrespective of the amount of current actually flowing through it. This introduces a source of error that needs to be taken into account. An alternative approach that considers current flow requires a self-consistency loop where each connection is re-weighted after the electrical network has been set up and solved. However, we find that this procedure leads to substantially longer computation times and still unsatisfactory precision, which is why we opt for a different method: We introduce an empirical parameter k to compensate for the approximations in our model. We initially determine this parameter by generating test structures and comparing multiple methods. The computation of k only has to be performed once and will be explained comprehensively in the following section. Before that, we demonstrate how the electrical network is constructed and how it is solved: We set up the conductance matrix (Gij ) of the electrical network according to nodal analysis. On the main diagonal, for every node i , we write the sum of all conductances of the edges that are connected to i
−1
−
Gin = σin·
Aflake ∑i Ai
tgr b l in
Gii =
(3)
∑j Gij,
i≠j
(5)
The off-diagonal elements of the conductance matrix are the negative conductances of the edges between the nodes i and j : −Gij ∀ i ≠ j . If nodes i and j are not connected, Gij is zero. To compute the total electrical conductivity of the GCM, we simulate a current through the structure. We define a border region at the edges of the structure where we connect all flakes to current sources (Fig. 5). At the points where current is injected, the current density is kept constant. This corresponds to homogeneous current injection at the side faces of the GCM. Furthermore, current is conserved, i.e. the same absolute amount of current Itot that enters one side leaves the structure on the other side. Apart from the border regions, the external current at every node is zero, in accordance with Kirchhoff’s current law. For a total of n nodes, we obtain the matrix form of the current-
Fig. 4. Trapezoid used for the calculation of the in-plane conductance. 366
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
The thickness h is given as the total number of 2D layers times the average layer thickness. The width w corresponds to the vertical axis of the sample in Fig. 5. 2.3. Adjusting the empirical parameter k As shown in our previous publication [14], the total conductivity σtot of an idealized GCM is ultimately limited by the packing density p of the GCM and the in-plane conductivity σin of the individual graphene flakes according to
σtot ≤ p ·σin ≡ σmax
The empirical parameter k has to be chosen such that this fundamental upper limit is respected. However, since the limit value p ·σin was derived for an idealized GCM, we need to investigate σmax for structures that are more realistic. We generate geometries as described above with varying packing densities in the relevant range for GCMs p > 0.70 . Subsequently, we compute the limit value of the total connetwork . From the same geometries, ductivity within our network model σmax we randomly select sections where we perform finite element analysis FEA (FEA) and compute the limit value of the total conductivity σmax . For FEA, we use the commercial software COMSOL Multiphysics® [20]. FEA is computationally far more expensive than our network approach, which is why we have to limit FEA to exemplary sections of the complete geometries. By considering multiple sections per structure (Fig. 6) however, we minimize the effects of statistical fluctuation. The results are compared to the limit value of the full structure computed with our network network approach σmax . FEA of various structures investigated with The conductivity limits σmax FEA are shown in Fig. 7a. As expected from theoretical prediction, the packing density and the in-plane conductivity are decisive factors for the limit of the total electrical conductivity. Apart from that, the homogeneity of the spatial distribution of graphene flakes plays a major role. Most notably, our FEA shows that for a GCM with homogeneously distributed graphene flakes, the limit value of the electrical conductivity corresponds to the limit value of the idealized GCM. This holds true for arbitrarily shaped and oriented flakes. Correspondingly, network we obtain the largest values σmax within our network model for homogeneous structures. Our FEA shows that inhomogeneities such as holes and clusters within a layer substantially reduce the total electrical conductivity, even if the overall packing density and the in-plane conductivity of the graphene flakes remain the same. From a physical point of view, these findings are coherent. Inhomogeneities lead to constrictions in the current path and therefore longer effective pathways, which in turn cause a smaller electrical conductivity. As a prime example of an inhomogeneous spatial flake distribution, we model a checkerboard structure (Inset Fig. 7). Note that this geometry is regular but it is not homogeneous in terms of spatial distribution. In a homogeneous
Fig. 5. Top view of a GCM. Current enters or leaves the structure through the yellow flakes at the sides. All the flakes that intersect with the red lines are considered for voltage measurement.
voltage-equation
⎛ G11 −G12 ⋯ −G1n ⎞ ⎛ Φ1 ⎞ ⎛ I1 ⎞ ⎜ −G21 G22 ⋯ −G2n ⎟·⎜ Φ2 ⎟ = ⎜ I2 ⎟ ⋯ ⋯ ⋯ ⎟ ⎜ ⋯ ⎟ ⎜⋯⎟ ⎜⎜ ⋯ ⎟⎜ ⎟ ⎜ ⎟ ⎝ −Gn1 −Gn2 ⋯ Gnn ⎠ ⎝ Φn ⎠ ⎝ In ⎠
(6)
with the conductances Gij , the electric potential Φi at node i and the injected current values Ii at node i . The conductance matrix is singular and we are free to choose a reference node where the potential is set to zero. We subsequently solve the equation for the potential vector so that we obtain the electric potential at every overlap in the structure. To compute the global voltage drop, we define two probe surfaces parallel to the side faces of the GCM where current is injected (Fig. 5). They are located outside the border region and separated by a distance d that is at least five times as large as the largest graphene flakes. The probe voltage is computed by averaging all the potentials at the overlaps that intersect with the respective probe surface. The potential difference V between the two probes corresponds to the voltage drop measured in a four-point probes configuration. For GCM simulations with tens of thousands of graphene flakes and sufficiently large d , the statistical error of this approach is negligible. We chose this form of measurement, as we are interested in the pure material parameters and not in contact resistances or voltage probe effects. However, our code could be used to simulate surface voltage probes of arbitrary sizes as well. With the empirical parameter k , the probe distance d , the thickness h and the width w of the GCM, the total electrical conductivity σtot can be determined according to
σtot = k ·
Itot hw · V d
(8)
(7)
Fig. 6. The central structure is generated by our geometry algorithm and evaluated with our network model. Randomly chosen sections are investigated with FEA. The color code indicates the homogeneous voltage drop across the sections. 367
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Fig. 7. (a) Limit of the total electrical conductivity versus overall packing density. The symbols mark different degrees of spatial homogeneity. The green line is the theoretical maximum. The red line is the linear fit for a checkerboard structure, shown as the inset. (b) Empirical parameter k versus packing density.
structure, the probability to find a flake at a certain position is constant across the whole 3D object. In a checkerboard structure on the other hand, the regular arrangement of square flakes leads to periodic constrictions across all layers. This reduces the overall electrical conductivity significantly. From our FEA results, we determine the empirical parameter k : As reference value, we choose the limit of the electrical conductivity of extremely homogeneous structures. For those FEA cases, σmax is equal to the theoretical maximum p ·σin . We choose k such that the same holds true for the network model and we obtain FEA network σmax = k ·σmax = p ·σin . The value of k is dependent on the overall packing density but does not vary with flake size or other microscopic parameters within the boundaries described in the geometry section. It ranges from 2.25 at p = 0.70 to 2.12 at p = 0.99. The full dependence is shown in Fig. 7b. Once the empirical parameter is adjusted, computations such as the ones for our previously published results can be performed. No further comparisons between methods are required and the high computational speed of our network simulation can be exploited. Fig. 8. Total conductivity of a GCM as a function of out-of-plane conductivity. Three characteristic regimes are indicated.
3. Discussion 3.1. Convergence and robustness of the method
plane conductivity. For this exemplary study, we investigated a system of 0.8 packing density, 25 layers and 39,000 flakes in total. Three regimes can be identified: For small σout , the overall conductivity is suppressed (Regime I), for very large σout , saturation is achieved (Regime III). The largest effects take place in Regime II, the transitional regime. Over a comparatively small range of out-of-plane conductivity values, the total conductivity moves from negligibly small to the saturation value. In Eq. (8), we already introduced the saturation value as the upper limit for the total conductivity of a GCM σmax . It is defined as the product of the packing density p and the in-plane conductivity σin . However, this simple mathematical dependence only holds true for highly conductive out-of-plane connections – a state that is achieved by large flake overlaps and high out-of-plane conductivities. For smaller out-of-plane conductivities, the impact of both p and σin decreases. This effect is visualized in Fig. 9. Fig. 9a shows the total conductivity of a GCM with respect to the in-plane conductivity of the individual graphene flakes. The different curves correspond to different out-of-plane conductivities σout . In general, a larger in-plane conductivity is always beneficial for the total conductivity of the macro-material. But the smaller σout , the more the overall conductivity levels off, even for the largest values of σin. The physical interpretation is straightforward: The connections between individual flakes become less and less conductive. At some point, even infinitely large in-plane conductivities cannot compensate small out-of-plane conductivities. The four upper curves in Fig. 9a are located in the transition regime as indicated in Fig. 8. Here, a
As our simulation method relies on random geometric elements, a minimum number of flakes has to be included in the simulation to obtain statistically reliable results. Furthermore, to treat a simulated GCM as one material and apply Eq. (7), the dimensions of the modeled GCM have to exceed a certain threshold. In extensive parameter studies, we analyzed the necessary conditions for successful modeling with our network method. In brief, a minimum number of 20–30 layers is required for convergence, depending on the structural details. This corresponds to a thickness of 6–10 nm, which is easily fulfilled by real GCMs. The length and width of the simulated GCM should be at least 25 times as large as the average flake length. For statistical robustness, a total of 30.000 flakes across all layers is required. The most important results of the convergence tests are illustrated in more detail in the Supplementary Material. If extremely inhomogeneous materials are simulated, the system size should be 1.5 times as large. With these parameters, the statistical fluctuation of the resulting total electrical conductivity can be limited to less than 2%. 3.2. Application of the method In this section, we use our method to study the interplay between several material parameters in a GCM. We focus on a critical regime in the parameter space that we did not analyze in our previous work. In Fig. 8, the total conductivity of a GCM is plotted as a function of out-of368
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Fig. 9. Total conductivity of a GCM as a function of (a) in-plane conductivity and (b) packing density for various σout . The symbols represent simulated data points while the lines indicate the trend. For (a), p = 0.8 and for (b), σin = 100 MS/m .
small increment in out-of-plane conductivity has a strong impact on the total conductivity. Fig. 9b shows the total conductivity of a GCM as a function of packing density for different values of σout . Higher packing density always corresponds to more current-carrying material and thus a higher total conductivity. But a minimum out-of-plane conductivity is required in order to exploit the material to full extent. Analogous to the in-plane-conductivity, the influence of the packing density decreases for smaller σout . Mathematically, this can be deduced from the slope of the curves, which decreases continuously with decreasing σout . We conclude that when fabricating a high-performance GCM, the material parameters have to be adjusted carefully. As different influencing factors compete, they overlap or even dominate each other. Experimental results are thus difficult to interpret and our simulation helps to understand the underlying physical reasons. Further comprehensive parameter studies and a comparison to experimental data can be found in our previous publication [14].
In conclusion, we present the methodological details of a modeling approach for the electrical conductivity of 2D-based conductors. We generate randomized geometries where we consider a variety of microscopic material properties such as in-plane and out-of-plane conductivity, packing density, or flake size. These parameters can be adjusted to fit the desired materials and the experimental conditions. Starting from the geometric structure of a material, we create a 3D resistor network that allows us to compute the overall electrical conductivity. By varying the microscopic parameters, their influence on the macroscopic electrical conductivity can be analyzed. We demonstrate this by investigating the complex interplay of packing density, in-plane and out-of-plane conductivity in a GCM. The method is fast and efficient compared to other simulative approaches such that large systems with several ten thousand flakes can be investigated.
3.3. Conclusions and outlook
Leo Rizzi: Conceptualization, Methodology, Software, Validation, Investigation, Writing - original draft. Andreas Zienert: Methodology, Writing - review & editing. Jörg Schuster: Conceptualization, Methodology, Resources, Writing - review & editing. Martin Köhne: Conceptualization, Project administration. Stefan E. Schulz: Conceptualization, Resources, Supervision.
CRediT authorship contribution statement
The strength and speed of our modeling method for 2D-based conductors come from its simplicity. By using simple geometric objects and a variable parameter set, we can adjust the computation to various physical situations. For pure graphene-only GCMs, the modeling is performed as described above. The conductivities should be varied according to statistical distributions across the whole GCM. Multilayers and defects are considered by increasing the thickness and adjusting the conductivity of the respective flakes. In principle, any geometry can be used as a starting point. To include dopants, the polygonal objects have to be regarded as graphene platelets decorated with non-carbon atoms or molecules. In-plane and out-of-plane conductivities need to be modified as well as the thickness of the platelets. Again, statistical distributions should be applied to decide which flakes are decorated with dopants and which are not. In our previous publication, we used this procedure to compare experimental data to our simulations and we obtained excellent agreement [14]. Going beyond graphene-based materials, the method can be used to model any bulk material composed of quasi-2D elements or similar ultrathin components. Examples might include conductive polymer nanocomposites or composites based on transition metal dichalcogenides. In our study, we concentrated on GCMs since they appear to be the most promising materials for conductor applications as of right now. To the best of our knowledge, no other simulation method is able to model similar systems more efficiently. Compared to FEA, our approach requires less time and is less memory intensive at system sizes that are ten times as large or more. Due to that, our approach enables the usage of statistical methods to reduce the fluctuation of the results.
Acknowledgment The authors would like to thank W. Pieper and F. Teichert for fruitful discussions as well as E. Lorenz for technical support. Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. Declaration of interest None. Data availability The raw data required to reproduce these findings cannot be shared at this time due to technical limitations. The processed data required to reproduce these findings cannot be shared at this time due to time limitations. However, all the raw input parameters required to reproduce our results are presented in this article. The processed data is fully presented in the form of graphs and plots. 369
Computational Materials Science 161 (2019) 364–370
L. Rizzi, et al.
Appendix A. Supplementary material
Phys. Lett. 91 (2007) 163513. [10] S. Stankovich, D.A. Dikin, G.H.B. Dommett, K.M. Kohlhaas, E.J. Zimney, E.A. Stach, R.D. Piner, S.T. Nguyen, R.S. Ruoff, Graphene-based composite materials, Nature 442 (2006) 282. [11] S. Stankovich, D.A. Dikin, R.D. Piner, K.A. Kohlhaas, A. Kleinhammes, Y. Jia, Y. Wu, S.T. Nguyen, R.S. Ruoff, Synthesis of graphene-based nanosheets via chemical reduction of exfoliated graphite oxide, Carbon 45 (2007) 1558–1565. [12] Z. Xu, C. Gao, Graphene chiral liquid crystals and macroscopic assembled fibres, Nat. Commun. 2 (2011) 571. [13] Z. Xu, C. Gao, Graphene fiber: a new trend in carbon fibers, Mater. Today 18 (2015) 480–492. [14] L. Rizzi, A. Zienert, J. Schuster, M. Köhne, S.E. Schulz, Electrical conductivity modeling of graphene-based conductor materials, ACS Appl. Mater. Interfaces 10 (2018) 43088–43094. [15] L. Peng, Z. Xu, Z. Liu, Y. Guo, P. Li, C. Gao, Ultrahigh thermal conductive yet superflexible graphene films, Adv. Mater. 29 (2017) 1700589. [16] Z. Wang, B. Mao, Q. Wang, J. Yu, J. Dai, R. Song, Z. Pu, D. He, Z. Wu, S. Mu, Ultrahigh conductive copper/large flake size graphene heterostructure thin-film with remarkable electromagnetic interference shielding effectiveness, Small 14 (2018) 1704332. [17] S. Gillies et al., Shapely: manipulation and analysis of geometric objects, 2007. [18] C. O'Callaghan, C. Gomes da Rocha, H.G. Manning, J.J. Boland, M.S. Ferreira, Effective medium theory for the conductivity of disordered metallic nanowire networks, Phys. Chem. Chem. Phys. 18 (2016) 27564–27571. [19] F. Dalmas, R. Dendievel, L. Chazeau, J.-Y Cavaillé, C. Gauthier, Carbon nanotubefilled polymer composites. Numerical simulation of electrical conductivity in threedimensional entangled fibrous networks, Acta Mater. 54 (2006) 2923–2931. [20] COMSOL Multiphysics® v. 5.3, COMSOL AB, Stockholm, Sweden. www.comsol. com.
Supplementary data to this article can be found online at https:// doi.org/10.1016/j.commatsci.2019.02.022. References [1] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Electric field effect in atomically thin carbon films, Science 306 (2004) 666–669. [2] K.S. Novoselov, V.I. Fal′ko, L. Colombo, P.R. Gellert, M.G. Schwab, K. Kim, A roadmap for graphene, Nature 490 (2012) 192. [3] P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M.C. Asensio, A. Resta, B. Ealet, G. Le Lay, Silicene: compelling experimental evidence for graphene-like two-dimensional silicon, Phys. Rev. Lett. 108 (2012) 155501. [4] P. Miró, M. Audiffred, T. Heine, An atlas of two-dimensional materials, Chem. Soc. Rev. 43 (2014) 6537–6554. [5] Q.H. Wang, K. Kalantar-Zadeh, A. Kis, J.N. Coleman, M.S. Strano, Electronics and optoelectronics of two-dimensional transition metal dichalcogenides, Nat. Nanotechnol. 7 (2012) 699. [6] L. Song, L. Ci, H. Lu, P.B. Sorokin, C. Jin, J. Ni, A.G. Kvashnin, D.G. Kvashnin, J. Lou, B.I. Yakobson, P.M. Ajayan, Large scale growth and characterization of atomic hexagonal boron nitride layers, Nano Lett. 10 (2010) 3209–3215. [7] A.A. Balandin, Thermal properties of graphene and nanostructured carbon materials, Nat. Mater. 10 (2011) 569. [8] J.-H. Chen, C. Jang, S. Xiao, M. Ishigami, M.S. Fuhrer, Intrinsic and extrinsic performance limits of graphene devices on SiO2, Nat. Nanotechnol. 3 (2008) 206. [9] J. Moser, A. Barreiro, A. Bachtold, Current-induced cleaning of graphene, Appl.
370