Modelling Carbon-Based Nanomaterials (CNMs) and Derived Composites and Devices

1. Introduction

Recent advances in analytical and computational modelling enable the prediction and understanding of material properties and responses on scales ranging from the electronic/atomistic, through the microstructure or transitional, and up to the continuum (Figure 1). Multiscale material modelling is based on a systematic reduction of the degrees of freedom at different length scales in which a material can be described. Connections between these scales are established either by parametrization, by grouping (i.e., coarse graining) or homogenization procedures. These different scales in the description of materials have historically been associated with different disciplines, from physics to engineering, including chemistry and materials science (Figure 1). Each scale can be studied by means of different types of models (electronic, atomistic, mesoscopic and continuum), which employ distinct basic entities and usually work in different ranges of timescales and number of entities (Figure 2). However, the scale and type of model do not have a one-to-one correspondence, since there are areas of overlapping regarding the models that can be employed for a specific analysis.

This report presents a review of different modelling strategies and methodologies used in the field of carbon-based nanoparticles (graphene, carbon nanotubes—CNTs, nanofibers, etc.), carbon-based nanomaterials (CNMs) and devices (such as sensors) that can be obtained from these materials. Within this framework, the variety of models, methods and applications is extremely wide, and getting a full picture in a condensed report is hard-to-reach. This review, then, aims to cover some of the key modelling approaches usually employed at different length and time scales, ending with a discussion on a specific type of sensor, derived from these carbon-based materials, as an example of device-level modelling. This leads to the four main sections that constitute the present document.

At the lower scale, the first type of models focuses on electronic calculations, theoretical studies of single-particle and collective excitations. These models are relevant for photonics, plasmonic and nano-electronics applications. At the next level, molecular dynamics (MD) models enable the study of material properties/responses that are mainly dominated by interactions between atoms or molecules, defined through force fields. These models help in predicting mechanical, thermal, surface properties, and other physical properties such as density, viscosity, or glass transition temperature; they are used transversally in the design of materials and nanocomposites for a wide range of applications. At the continuum level, modelling is typically employed to predict the macroscopic-homogenized properties of materials by considering the constituents as a continuous medium rather than as discrete particles. Compared with the latter, MD models provide detailed atomic-level information but are limited to small systems due to computational constraints. Continuum models, on the other hand, can handle much larger systems but lack the detailed description of atomic or molecular interactions. These models can also be used to predict different types of properties (again mechanical, thermal, electrical, acoustic, and other physical properties). Certain kinds of products and devices need to be modelled as systems instead, which might require a combination of different type of properties and physics. Each system has its own requirements, which can lead to the integration of different sub-models. At this scale, as an application example relevant to sensors, this review focuses on the modelling and simulation of graphene field effect transistors (GFET) devices and circuits.

This review has been performed within the framework of COST Action EsSENce (High-performance Carbon-based composites with Smart properties for Advanced Sensing Applications—CA19118), funded by the European Cooperation in Science and Technology (COST). The EsSENce COST Action aimed to construct an international scientific & technological innovation hub focusing on advanced composite materials reinforced with CNMs for sensing applications. Particularly, working group WG5 deals with the modelling of materials and sensing mechanisms, identified as a key area to support the development of CNMs, advanced high-performance composites, and their implementation into products and devices. Thus, the objective is to exemplify and promote i) the use of “standard” or widely shared modelling protocols, or ii) the development/application of new modelling approaches. In addition, a second goal is to identify some of the needs in this field, considering both experimental characterization (needed to feed or validate the models) and the modelling perspective.

2. Theoretical Modelling of Interactions of Charged Particles with Graphene-Based Nanomaterials and Their Composites

Graphene is the best known and the most explored two-dimensional (2D) material. It is the world’s thinnest and strongest material, with the highest electrical and thermal conductivity known. Due to its outstanding physical, chemical, electrical and optical properties, graphene shows great potential in many fields, including sensors.

In nanoscale devices graphene typically appears in stacks separated by insulating layers of finite thickness [2,3,4], which usually support strong Fuchs-Kliewer or optical surface phonon modes [5]. Those phonon modes are active in the terahertz (THz) to mid-infrared (mid-IR) frequency range and can dampen the Dirac plasmon in doped graphene which operates in the same frequency range [6], or can hybridize with it [7]. As a prototype of layered nanostructures involving doped graphene sheets, a sandwich-like composite which consists of two layers of graphene separated by an insulating layer was studied and it was found that the structure supports a variety of interesting plasmon-phonon hybrid modes in the THz to mid-IR frequency range [8]. An insulating layer of sapphire (aluminum oxide, Al₂O₃) was chosen because it is often used as a dielectric spacer in experiments [9,10,11,12,13,14]. The effective dielectric function of the system was obtained by using a local dielectric function for the bulk Al₂O₃ and by using two approaches within the random phase approximation for graphene’s electronic response: an ab initio method based on the time-dependent density functional theory calculations and a method based on the massless Dirac fermion approximation for graphene π bands.

The most efficient way to probe the plasmon-phonon hybridization between graphene and the nearby insulator(s) is by means of an externally moving charged particle. Such hybridization has an effect on the energy loss of an incident particle [15], as well as on the resulting wake effect in the induced electrostatic potential [16]. The wake potential (the total potential in the plane of the upper graphene sheet) induced by an external charged particle moving parallel to the graphene-Al₂O₃-graphene composite system (Figure 3) was investigated in Ref. [17]. The distance of the charged particle from the top graphene, the thickness of the Al₂O₃ layer, and the doping density (i.e., Fermi energy) of graphene were fixed at their respective typical values. The effects of variations of all those parameters were studied in Ref. [18].

The plasmon-phonon hybridization also has an impact on the stopping force (the dissipative force which opposes the particle’s motion) and the image force (the perpendicularly oriented conservative force which bends the particle’s trajectory towards the upper graphene sheet). A thorough analysis of the stopping and image forces on a charged particle moving parallel to the graphene-Al₂O₃-graphene composite was performed in Ref. [19], covering broad ranges of the particle speeds and distances, as well as the doping densities of the two graphene sheets. It should be noted that the stopping force is the negative of the usual stopping power, whereas the image force is related to the familiar image potential.

Note that the wake effect and the stopping and image forces were previously investigated theoretically in free-standing and supported graphene [15,16,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39].

In the last few years van der Waals (vdW) heterostructures based on graphene and hexagonal boron nitride (hBN) layers with different stacking modes have attracted a great deal of interest because of their potential applications [40,41,42,43,44,45,46,47,48]. Graphene/hBN vdW heterostructures were studied very recently in Refs. [49,50].

Electron energy loss spectroscopy (EELS) is a commonly used experimental technique for investigating electronic and plasmonic properties of 2D materials and vdW heterostructures [51,52,53]. Theoretical modelling of the experimental EELS data for free-standing (single and multilayer) graphene sheets obtained by (scanning) transmission electron microscope is presented in Refs. [54,55,56,57], whereas the theoretical modelling of the experimental EELS data for monolayer graphene supported by different substrates is given in Refs. [58,59,60,61].

Acoustic plasmon (AP) in graphene or in graphene-dielectric-metal structure has been studied very intensively in the last few years [62,63]. In Ref. [64] the authors focused on the AP in graphene doped by alkali metals and demonstrated that two isoelectronic systems, KC₈ and CsC₈, support substantially different plasmonic spectra: the KC₈ supports a sharp Dirac plasmon (DP) and a well-defined AP, while the CsC₈ supports a broad DP and does not support an AP at all. These findings could be very useful in the area of chemical or biological sensing [65,66].

3. Molecular Dynamics Applied to CNM Properties Prediction

Molecular Dynamics (MD) is a computational method to simulate the behavior of atoms and molecules over time. In an MD simulation, the trajectories of atoms (considered as rigid spheres) and molecules are determined by numerically solving Newton’s equations of motion. Forces between particles are calculated by exploiting calibrated interaction potentials (force fields), and from these forces the accelerations, velocities, and subsequent positions of the atoms are determined [67]. For materials science, MD provides detailed insights into phenomena that would otherwise be difficult to observe directly, such as diffusion, phase transitions, interfacial effects, and mechanical behavior at the nanoscale.

One of the key aspects of MD simulations lies in the force fields, which are mathematical models used to calculate interactions between atoms due to the presence of covalent bonds (bonded interactions) and electrostatic forces (Coulomb/Van der Waals, non-bonded interactions) [68]. In MD simulations of CNMs, some specific force fields have gained prominence [69]:

• The Adaptive Intermolecular Reactive Bond Order (AIREBO) potential is tailored for carbon systems and describes long-range van der Waals interactions and torsional effects. It is versatile for modelling both sp² and sp³ hybridized carbon structures [70]. AIREBO might not perform well for systems with significant charge transfer or in the case of interactions with elements outside its parameterization.
• Tersoff potential considers both the distance between atoms (bond lengths) and their relative orientation (bond angles) to provide a detailed representation of the complex interactions that occur in carbon-based materials [71]. This potential may not be ideal for modelling weak interactions, and it might require recalibration for systems different from its original parameterization.
• ReaxFF is a reactive force field capable of simulating bond formation and breaking during MD simulations. This dynamic nature is achieved by not predefining specific bond types but allowing the system to evolve based on atomic positions and interactions. Due to its reactive nature, ReaxFF can be computationally demanding. It also requires careful system-specific parameterization to ensure reliable results, e.g., in the case of condensed carbon phases [72].
• Machine Learning (ML) interatomic potentials differ from traditional ones, as they do not depend on fixed mathematical formulas. Instead, they learn representations of the potential energy surface of the system through trainings based on lower-scale simulations. Several implementations for certain carbon forms with near DFT-level accuracy have been reported in the literature, e.g., Gaussian Approximation Potential (GAP) [73], hybrid neural network potential [74], GAP-20 potential for various crystalline phase of carbon and amorphous carbon [75]. Furthermore, MACE—a transferable force field for organic molecules created using ML trained on first-principles reference data—was recently implemented [76]. Despite the good accuracy of current ML-based force fields in predicting the properties of carbon allotropes, various challenges still exist, especially regarding the description of mechanical properties and the curation of reliable training datasets.

MD is well-suited for investigating various properties of CNMs and composites made thereof [77]. The mechanical properties of these materials, for example, can be determined through MD simulations that allow the stresses and strains experienced by the system to be evaluated [78]. In these tests, a strain is systematically applied to the system, and the resulting stress responses of the material are recorded, thus providing information on its elastic constants, tensile strength, and potential fracture points [79]. For instance, simulations of cellulose nanocrystal-graphene composites revealed enhanced mechanical properties due to covalent bonding and van der Waals interactions [80]. Similarly, MD analyses of single-walled carbon nanotubes (SWCNT) have demonstrated Young’s moduli in agreement with experimental values, showcasing their exceptional mechanical strength and stiffness [81]. However, these simulations face some challenges, as the simulated strain rates should be much higher than those typically found in experimental setups to provide meaningful results in a feasible computational time. Another technique for analyzing the mechanical properties of a material is nanoindentation, which involves the simulation of a virtual indenter pressing on the surface and allows the hardness and localized stress response to be derived [82]. For example, nanoindentation MD simulations on polymer nanocomposites highlighted the effect of nanoparticle interactions and temperature on mechanical reinforcement [83] The results, though, can be influenced by the chosen shape of the indenter and the interaction potentials utilized. Pull-out tests simulations [84] and density profile analysis [85], on the other hand, are valuable aids in characterizing the interface region between a CNM and the surrounding environment (e.g., the polymer in the case of composites), which is a key interaction region that strongly determines the properties of the composite material. These analyses provide useful insights into the mechanisms of load transfer through materials and interfacial adhesion behavior.

The thermal properties of CNMs can be explored using various protocols [86,87]. In the Non-equilibrium Molecular Dynamics (NEMD) technique, for instance, a temperature gradient is established within the simulation domain, thus enabling the calculation of thermal conductivity [88]. NEMD was employed to calculate the thermal conductivity of multi-walled CNTs with different geometrical features, such as diameter, length, chirality, and number of walls [89]. Still, this method entails an issue, as the artificial imposition of a gradient might not realistically replicate actual experimental scenarios. The Equilibrium Molecular Dynamics (EMD) method offers another approach, relying on the analysis of heat current fluctuations within a system at equilibrium. A notable example involves the use of EMD simulations to determine the thermal conductivity of graphene nanoribbons [90]. By employing the Green-Kubo method, researchers studied the effects of ribbon width, edge roughness, and hydrogen termination on thermal conductivity. The results showed that smooth edges yield the highest conductivity, while edge roughness significantly reduces it. However, EMD often requires extended simulation times. Another critical thermal property that can be computed by MD simulations is the thermal boundary resistance (TBR, also known as Kapitza resistance, originating from phonon scattering in the presence of defects or interfaces) [91]. As carbon nanomaterials are often embedded into other materials (e.g., polymeric matrices, fluids), understanding the efficiency of heat transfer across these interfaces becomes vital [92]. Both NEMD and EMD protocols offer quantitative insights into this property. Furthermore, theoretical models such as the Acoustic Mismatch Model and the Diffuse Mismatch Model, informed by inputs from MD simulations, provide complementary perspectives on TBR [93].

For thermodynamic properties, MD simulations lean on specialized techniques, like free energy calculations. Properties such as adhesion can be probed using advanced sampling methods, such as umbrella sampling and metadynamics [85], but these techniques often require care in choosing the force field and can be computationally burdensome. The wetting properties of CNMs are fundamental as well for optimizing their performance in suspensions and composites. Wettability, quantified by contact angle measurements, can indeed affect the dispersion of CNMs in various solvents, along with the ability of polymers to spread and adhere to CNMs thus influencing the mechanical and thermal properties of the composite [94]. At the atomistic level, two main approaches are used to measure the contact angle. The free energy perturbation (FEP) method involves calculating the free energy changes to determine the interaction parameters between a liquid and a surface [95]. This approach allows the interaction parameters to be calibrated and the work of adhesion and friction coefficients to be evaluated. In the second method, called the droplet method, a liquid droplet is placed on the tested surface and allowed to relax until equilibrium is reached [96]. The contact angle is then measured by analyzing the shape of the droplet along the three-phase contact line [97]. The latter approach is widely used to study the effects of surface roughness and interfacial properties. For instance, MD simulations revealed that difunctional epoxy and cyanate ester resins exhibit high wettability on CNT surfaces, while polyether ether ketone resins show poor wetting properties [98]. Graphene oxide, due to its hydrophilic functional groups, generally offers better wettability than pristine graphene, improving dispersion and interfacial bonding in polymer matrices [99]. Additionally, carbon nanoparticle coatings synthesized through controlled flame deposition can be tailored from hydrophilic to superhydrophobic states by adjusting synthesis conditions, which optimize surface interactions and enhances composite performance [100].

Other type of properties that can also be studied through MD models are the dielectric ones (dielectric constants, relaxion) using specific approaches as the dipole moment fluctuation method [101].

While offering mechanistic insights into the properties and behavior of CNMs, MD simulations do not lack challenges. For instance, MD approaches typically operate within specific temporal (fs to µs) and spatial (nm to µm) scales, and phenomena outside these scales might not be detected and analyzed effectively [102]. Moreover, the choice of the force field can significantly impact the results, as not all force fields describe the intricate interactions in nanomaterials with the same accuracy [103]. Furthermore, high-resolution simulations, especially those involving long timescales or large systems, can be computationally demanding, requiring substantial resources and time. The imposition of periodic boundary conditions, instead, can lead to artifacts in the results, especially if the size of the simulated system is not large enough compared to the phenomena under observation [104]. Finally, the rates used in simulations (e.g., thermal or mechanical ones), due to computational constraints, often exceed experimental rates, potentially leading to discrepancies between numerical and experimental data [105].

Looking ahead, the MD research landscape on CNMs presents promising prospects [106]. One of the most significant improvements in recent years has been the development of enhanced force fields. Ongoing research in this area aims to refine these force fields specifically for carbon nanomaterials, with the goal of increasing the accuracy in predicting their properties. Moreover, the growing popularity of multiscale modeling as a robust approach offers significant opportunities: by linking/coupling MD with other simulation methodologies, such as electronic, mesoscopic and/or continuum models [107], researchers aim to bridge the spatiotemporal gaps between scales. Yet, one of the most transformative shifts in the field of molecular dynamics is probably the integration of data-driven approaches within MD simulations. The support of machine learning and artificial intelligence is not merely augmenting computational efficiency, but it is reshaping the paradigms of simulations. These tools offer optimized parameter selections, predictive capabilities, and the prospect of devising new force fields. Lastly, shared efforts between experimentalists and computational researchers are fostering iterative refinements in simulation methodologies. This synergy is guiding MD studies closer to experimental observations, ensuring a more harmonized understanding of CNMs.

4. Continuum Models

At the next level within the multiscale materials modelling framework, continuum models are powerful tools for predicting the properties of nanocomposites. These models consider the microstructure of the nanocomposite, including the size, shape, and distribution of the nanoparticles, to predict the macroscopic-homogenized properties of the material. Continuum models can be used to predict different types of properties or responses: mechanical, thermal, electrical, acoustic, among other physical properties.

There are basically three types of continuum models: analytical, semi-analytical and numerical. The most common analytic models used for nanocomposites are based on the early work of Eshelby on ellipsoidal inclusions on an infinite elastic matrix, such as the relevant Mori-Tanaka’s mean field approach, which introducing as a hypothesis that each particle sees as far field the average matrix strain, allowed the application to non-diluted concentrations, calculating the properties of the composite by a physical combination of the properties of the nanoparticles and the matrix, weighted by their respective volume fractions. Other type of analytical models are the self-consistent models (SCM), where each nanoparticle in the nanocomposite is embedded in a matrix that has the same properties as the nanocomposite itself (this leads to a system of coupled equations that must be solved numerically); In general, these models can be solved relatively quickly compared with other approaches, making them useful for initial approximations and quick calculations. Despite this, as it is reported later, more and more advanced models are being developed attempting to capture complex geometric and material aspects. There are hundreds of works in the literature applying these models to determine standard fiber-based composites properties first, and later applied also to nanocomposites, for example considering the matrix-particle interphase as a third material phase constituent. In this sense, interesting references comparing/reviewing all these methods are the article of Y. Wang and Z. Huang [108], the chapter of Yehia A. Bahei-El-Din about Averaging Models of Fibrous Composites in [109], and the recent publication of A. Elmasry, W. Azoti, S. El-Safty and A. Elmarakbi that presents a review comparing different models for effective properties calculation of nano- and micro-composites [110]. Works focused on mechanical properties include: i) nanocomposites’ stiffness prediction (i.e., response in the elastic range), as the research done by M.M. Shokrieh, M. Esmkhani, Z. Shokrieh and Z. Zhao for an epoxy resin modified with graphene nanoplatelets [111], the work presented by A. Chiminelli and M. Laspalas for an epoxy resin modified with MWCNTs [112], the work of A. Singh and D. Kumar studying the influence of the functionalization of graphene nanoplatelets in the elastic properties of a modified polyethylene [113] or a more recent work of D. Shin, I. Jeon and S. Yang for graphene modified PET [114]; ii) non-linear behavior and strength predictions of CNMs, as the approach proposed by J. Nafar Dastgerdi, G. Marquis and M. Salimi introducing interfacial damage/debonding processes in CNTs reinforced polymers [115] or the model developed by W. Azoti and A. Elmarakbi applied to a graphene platelets GPL-reinforced polymer PA6 composite (Generalized Mori-Tanaka) [116]. Other works/studies can be found about the utilization of this type of analytic approaches to develop more advanced models introducing viscoelasticity [117] or creep effects [118], among others. Overall, good correlations with experimental results are obtained for the elastics properties especially at low nanoparticles concentrations, while at higher concentration more significant deviations are usually observed. This is generally explained due to interactions between the nanoparticles at high contents that produce non-homogeneous dispersions and agglomerations. For predictions in the non-linear regime and in terms of strength, it is generally seen that the results are strongly dependent on the particle/polymer interface representation and on the aspect ratios of the particles. This also observed in the modelling of other type of properties of nanocomposites (thermal, electrical, acoustic…). For example, in the work published by J. Shao et al. [119] for predicting dielectric properties of polymers modified with nanoparticles, the Knott model was modified with the Tanaka formula for hybrid particles to take into consideration the interphase region. Obviously, the issues mentioned at high concentrations are also present in these cases. Apart from this, these models are valuable to evaluate different functionalization of carbon-based nanoparticles.

In addition to the models mentioned above, some semi-empirical analytical models are also well-known for composites, starting with the classical rule of mixture, evolving to Chamis’ model and the Halpin–Tsai’s model (as a simplified version of the SCM), and leading to more advanced models introducing (again) elastoplasticity (approaches to yield stress and linearization methods [108]). The Halpin–Tsai approach for aligned reinforcement has been employed widely for the analysis of graphene-based nanocomposites [120]. In this sense, works as the ones developed by Weon et al. [121], Chong et al. [122] and Zarasvand and Golestanian [123] can be highlighted. In the last, the nonlinear tensile stress-strain behavior of randomly-distributed graphene nanocomposites have been obtained, presenting a good correlation with experimental results in the whole range of the stress-strain curves. In a work published by M. Yang et al. [124], the Halpin-Tsai model was adapted to quantitatively characterize the effect of temperature on the yield strength of nanofiller-reinforced polymer-matrix nanocomposites. Compared with the classical ones, the model showed better agreement with the available experimental data from sub-zero temperature to the full glass transition region. Progressing with elastoplasticity, several Non-linear Mean Field Methods have been developed from the elastic ones, the Tangent and Secant approaches being the most known. Some reference publications implementing these methods are the ones of Doghri [125,126]. An advantage of these models is the straightforward implementation. They also allow to study more complex loading cases than other models, including cyclic loads.

Semi-analytical methods are based on global constitutive equations that are evaluated from the local scale using analytic/explicit relations that link the microscopic and the macroscopic properties. Analytical relations are usually dependent on mean field procedures. The best known is the Transformation fields analysis (TFA). It connects analytical and computational approaches by a computational evaluation of the localization operators. The main concept is replacing the plastic strain field with piecewise uniform fields to reduce the number of macroscopic internal variables. Thereby, a set of reduced constitutive relations for the heterogeneous material can be established, leading to a computational time reduction compared to full numerical models [110]. Reference works of this type of model are the ones published by Dvorak [127,128]. A more recent work published by I.Al. Khattab and M. Sinapius present an interesting implementation of the TFA model as a user routine integrated into RVE-micro scale model [129]. The results reveal that the TFA is a proper method for solving inelastic deformation and other incremental problems in heterogeneous media with many interacting inhomogeneities on a nanoscale level.

Finally, numerical continuum micromechanics models, such as the finite element method (FEM), are considered also powerful approaches for CNMs modelling. In numerical approaches, an acknowledged constitutive model (elasticity, viscoelasticity, elastoplasticity, viscoplasticity…) is assumed in a Representative Volume Element (RVE) to induce an explicit macroscopic model. On one hand, these models allow to introduce more accurate description of the materials nano-morphology [130]. This is one of the limitations of analytic models, that usually require an idealization of the nanoreinforcements’ shape, representing them as discs, cylinders or spheres. On the other hand, these models allow to introduce complex nonlinear multi-phase material behaviors. The main disadvantages of finite elements are computational costs, size dependency of the results and limitations in the development of RVEs with high inclusions volume concentrations and aspect ratios.

RVE based models apply to statistically homogeneous materials. They can be used as repeating unit cell (RUC) when the (nano)composite has a periodic microstructure. Sufficient number of randomly distributed fibers or particles to be contained in the RVE is needed so that the microstructure of a composite could be reflected precisely [108]. This is also linked with the size of the volume studied. Various works reflect the importance of defining a proper RVE size [131]. In this sense, it happens that a larger RVE size gave better prediction accuracy but resulted in lower computational efficiency. Finally, the boundary conditions are a critical aspect in RUC models, and they have to be defined carefully to represent properly the effect of the reinforcements/particles distribution patterns and the loads applied.

Some first reference publications that can be found about RVE and FEM to estimate properties of polymer nanocomposites are the works of Liu and Chen [132,133]. Particularly, they studied the elastic properties of CNT-reinforced polymers. Chwał and Muc [134,135] applied a similar approach with various boundary conditions to calculate the mechanical properties of SWCNT-polymer nanocomposite. More recent works present a combination of this type of FE continuum models with MD, in similar way to the multi-scale framework presented in this report. For example, recently Barakat et al. published an article where the distribution of mechanical properties of graphene-based polymer nanocomposites is computed using a micro-meso up to macro hierarchical computational approach employing non-equilibrium atomistic MD simulations and continuum FE models [136]. In the same line, the work published by Muhammad et al. can be highlighted [137], where not only mechanical but thermal properties are predicted for graphene-reinforced epoxies. For SWCNT-polymer nanocomposites, Malague et al. proposed a procedure to assess size effects using the atomistic simulations and equivalent continuum model with a large number of CNTs [138].

Elastic RVE/RUC models can be extended to nonlinear cases, providing that the nonlinear constitutive laws for the constituents are available. For example, based on a three-dimensional RVE model, Yuan and Lu [139] conducted a numerical investigation on the elastoplastic behavior of carbon nanotubes (CNTs) reinforced polymer composites. A. Zarasvand et al. [140] also conducted experimental, numerical, and micromechanical studies to determine the nonlinear behavior of CNT-reinforced polymer. When non-linear constitute behaviors are considered, FEM approaches become computationally consuming. To address this limitation, several numerical strategies with reduced computational effort have been developed. Examples are the Voronoi Cell Finite Element Method (VCFEM) [141,142], the Generalized Method of Cells (GMC) [143], and the Finite Volume Direct Averaging Micromechanics (FVDAM) [144,145,146].

Figure 5. Examples of continuum models to study CNMs. A) Mean field homogenization scheme [112], B) Combined FE and analytic Mori-Tanaka models for CNT reinforced polymer taking into account viscoelastic properties [117], C) Utilization of Halpin–Tsai approach for Graphene Nanoplatelet reinforced epoxy [123] and D) REV/FEM models analyses considering elastoplasticity for wavy and random CNTs distribution in nanocomposites [139].

Figure 5. Examples of continuum models to study CNMs. A) Mean field homogenization scheme [112], B) Combined FE and analytic Mori-Tanaka models for CNT reinforced polymer taking into account viscoelastic properties [117], C) Utilization of Halpin–Tsai approach for Graphene Nanoplatelet reinforced epoxy [123] and D) REV/FEM models analyses considering elastoplasticity for wavy and random CNTs distribution in nanocomposites [139].

5. CNM Devices—Graphene Field Effect Transistors

When we arrive to discuss carbon nanomaterial devices that can be considered mature enough to be used in sensing systems, one of the most used configurations is based upon graphene field effect transistors (G-FET). High frequency applications have been in a first time a major goal for this class of devices, and a series of interesting results reporting the fabrication of research level devices working in the GHz domain were published in the 2010’s years [147,148,149,150]. Anyway, the problems that were initially addressed for the graphene transistors to replace silicon MOSFET as the choice alternative for commercial high-performance logic and high frequency electronics [151], prevented until now their large-scale adoption for this class of applications.

As a follow up, the intense research work realized during these years has produced many interesting results, experimental and theoretical as well, that have permitted a good insight into the physic mechanism of the GFET structures, paving the way to the definition of an application playground in the biosensing domain, avoiding the HF requirements of a G-FET based ASIC project. A nice summary of the present state of this approach can be found in the recently published review paper of S. Szunerits [152]. According to the authors not only graphene is proposed for FET biosensing applications, but also other 2D materials like metal dichalcogenides, hexagonal boron nitride or black phosphorus have been investigated for FET gate-channel technology. Indeed, graphene maintains his position as a first-choice material commercial applications for low term potential.

Although it is not currently available on the market, it may be foreseen that the use of GFET circuitry will soon find commercial application in wearable flexible sensing systems [153], in Point-Of-Care systems [154,155], and even in skin bioelectronics for health monitoring [156] exploiting its mechanical/electrical properties together with biocompatibility characteristics. Advanced experimental results in this field have been recently presented for applications of major interest, like, for example, cancer monitoring [157], COVID19 screening [158] or glucose monitoring in diabetic patients [159].

While no clinical validation has yet been performed, the proposed technology for these important biomedical playgrounds has been fully demonstrated in laboratory environment. The next step still needed for allowing the fabrication of large-scale point-of-care biosensing systems is the definition of a reproducible fabrication process, integrated with standard silicon CMOS technology. Once some variability can be expected from the fabrication condition, modelling and simulation will play a major role in the G-FET circuitry design. Looking at the recent past, modelling and simulation of GFET device and circuits, in a variety of configurations has certainly given an important contribution to the actual advanced state of the art for design and fabrication.

Regarding the electrical properties, a model to calculate the DC characteristics of large-area graphene field-effect transistors was presented initially by Thiele and Schwierz [160]. This paper opened the door to the heuristic approach of correlating the experimental results with the specific characteristics of the graphene material, like the sheet carrier density and transport properties. More specifically, the Thiele&Schwierz model suc-cessfully describe the experimental DC behavior of the drain current (Id), as a func-tion of the Drain-Source Voltage (Vds), for different Gate Source voltage (Vgs), using a reduced set of parameters about the device geometry like the gate length (L) and width (W), the insulator thickness (tox), drain and source contact resistances (Rs, Rd), together with material properties like charge mobilities in the graphene channel, sheet carrier density, carrier saturation velocity.

As a general description, a G-FET has a structure that can be thought similar to a thin film transistor (TFT), where the 2D graphene layer is used for channel forming. The traditional approach for the electrical simulation of semiconductor FET devices is based on the drift-diffusion model, including the Poisson equation [161]. Following an approach similar to the one proposed by Thiele&Schwierz, yet including the standard drift-diffusion approach, a Verilog compact model suitable for circuit level design model based of the drift and diffusion scheme has been presented and made freely available by Landauer [162]. The main problem for the direct application of this model to the G-FET structure is the definition of the charge mobility value, that in a 2D material can be largely affected by the presence of impurities, lattice defects and substrate quality. A nice approach to this problem was recently described by Nastasi [163]. They describe extensively the physical and numerical model, including the model for the mobility and its dependence on the applied electric field, and the results obtained on a top gate configuration described by a 2D geometry.

The approach directed to modelling the behavior of the device when embedded in a more general circuit is of major importance when the target application is related to the sensing of some specific bio-element. From this point of view, the correlation between material properties and electrical behavior and output extraction is the main objective of any modelling task. The G-FET structure better suited for targeting the sensing of biomolecules (such as proteins and nucleic acids) in biological fluids, is the co-planar and liquid-immersed-gate configuration, which is effectively largely preferred in biological G-FET design. This configuration overcomes the standard top-gate or bottom-gate inherited by the flat panel semiconductor industry and presents a new challenge for modelling and simulation approach. The main results of this activity have the production of G-FET models suitable to be used in a circuit simulator, to support the fabrication of complete G-FET based circuit for biosensing applications. A dual gate G-FET model for circuit simulation, suitable to be used for biosensing use, was presented by Umoh. The main value of this model is the direct application in a SPICE implementation, in a traditional configuration settling on standard FET SPICE model parameters, used internally to calculate junction resistances and capacitances. Within this context, Jmai [164] presented a model freely distributed in its MATLAB and Verilog implementation, allowing the user to select an appropriate topology for a system-level design suitable to be used in real life applications. Once again, the charge mobilities are a main parameter that need to be assessed as a function of the fabrication condition. As already mentioned, considering the multiscale concept, these models at system level can be fed with results obtained from models at lower scales, in the specific case of GFETs for example to determine the dielectric properties that are required. In this way, the nanocomposite properties can be investigated and designed with MD or continuum models and these results can be introduced in the transistor model for device optimization.

As a matter of fact, from the biosensing point of view, it is possible to assume that the presence of the target bio-element on the graphene surface would be the main driving cause for the mobility fluctuation, and the correspondent variation of the drain current of the transistor. This observation focusses the modelling attention on a differential current extraction, targeting the information about the presence of the biomarkers, instead that the mobility value itself. The model for the dual gate configuration is directly extendible to the liquid gate configuration, and from a high-er-application-level approach, Fuente-Zapico has adapted the Synopsys Sentaurus Commercial TCAD suite for the simulation of a liquid gate graphene field effect transistor, GFET, used as antibody-based biosensor [165]. In this work the authors have successfully included the liquid electrolyte, the antibody functionalization of the graphene surface and the biomarkers trapping effects, reporting application driven parameters to the current-voltage characteristics, calculated in a standard commercial device-circuit simulator.

As a bottom line, it must be outlined that The Graphene Flagship’s 2D Experimental Pilot Line (2D-EPL), has been offering, starting in 2022, prototyping services to academics, SMEs and companies which can benefit from the progresses of graphene related materials integration with silicon [166]. The 2D-EPL provide in a multi-project wafer run (MPW) for G-FET circuit including a top/bottom contact with an optional local or global back and liquid gate. Directed to Bio/Gas/Chemical sensors this MPW include a stand-alone G-FET circuit as well as a fully CMOS integration.

Figure 5. G-FET in a typical dual gate configuration (A) and for biosensing application with the analyte in contact with the graphene layer (B). The application of a voltage at the liquid analyte permits a state-of-the-art dual liquid gate configuration.

Figure 5. G-FET in a typical dual gate configuration (A) and for biosensing application with the analyte in contact with the graphene layer (B). The application of a voltage at the liquid analyte permits a state-of-the-art dual liquid gate configuration.

6. Conclusions

A review of different modelling strategies and methodologies used within the framework of carbon-based nanoparticles (graphene, carbon nanotubes—CNTs, nanofibers, etc.), and carbon-based nanomaterials (CNMs) and devices (as sensors) that can be obtained from these materials has been presented. The range of modelling and simulation techniques generally employed covers from electronic to continuum models, passing through atomistic and meso-scale. For each case, their typical uses, some of their weaknesses and strengths, and the latest developments are presented, including application examples and key references. Sometimes, multi-scale workflows are required, where these different models and scales are sequentially linked as part of a single analysis. Modelling devices require multi-physics models, integrating different phenomena that can be present simultaneously. Then, system type models are used. In this sense, the modelling of G-FET for sensing applications is included in this review, as specific example of carbon-based devices.

As regards the theoretical study of interactions of charged particles with graphene-based nanomaterials, a mathematical framework for studying such interactions uses various response functions of such materials as an input, which may be available from analytical phenomenological models, ab initio calculations, or even experimental data from optical spectroscopies. This allows various research groups to quickly develop reliable models for analyzing plasmonic properties of layered nanostructures including graphene, which can be used to predict the properties for novel designs of such structures.

Regarding molecular dynamics models, they are typically used to investigate various properties of CNMs at the atomistic scale, thus helping to predict their physico-chemical behavior and optimize their performance for specific applications. MD simulations can provide valuable insights into the thermomechanical properties of car-bon-based structures (e.g., thermal conductivity or elastic constants) and interface interactions in composite materials, yielding increasingly accurate results due to ongoing research efforts and the integration of data-driven approaches.

Continuum models include mainly analytical, semi-analytical and numerical approaches. Analytical models can be solved relatively quickly and often provide a clear and straightforward understanding of the underlying physics and insights into the relationships between different variables and parameters. However, sometimes they are limited in terms of the geometrical complexity that they can address. Numerical models can partially deal with these issues, but are constrained by their computational cost, size dependency of the results and restrictions in the development of RVEs with high inclusions volume concentrations and aspect ratios. Overall, good reproducibility of the properties of the materials are obtained especially at low nanoparticles concentrations, even for non-linear and time dependent responses, while at higher concentration more significant deviations are usually observed. This is generally explained due to interactions between the nanoparticles at high contents that produce non-homogeneous dispersions and agglomerations. For mechanical properties predictions in the non-linear regime or in terms of strength, it is generally seen that the results are strongly dependent on the particle/polymer interface representation and on the aspect ratios of the particles. Similar influences are also observed for other types of properties, such as thermal ones.

Finally, electrical behavior of G-FET, their analysis and design can be supported through Spice and Verilog models made freely available by the authors. The theoretical formulation of such models is available in literature and the main simulation software houses are starting to include them in their commercial packages. We can expect, for the next future, the introduction of integrated models, describing mechanical, thermal electrical and optical models, addressing the sensing applications where this class of devices can offer breakthrough solutions. Modelling techniques are powerful tools to support materials research in the development of novel applications, particularly in the case of carbon-based materials and their derived composites and devices. It provides the key information for identifying new materials, tailoring materials and/or design materials for structures and systems.