The bond between carbon nanotubes (CNTs) and cementitious materials is the key characteristic for predicting theﬂexural strength. However, CNT dispersion quality may change the bonding mechanism. A probabilistic approach can bene ﬁt the deterministic models to capture the uncertainties affecting these characteristics. This study proposes a probabilistic model using the Kelly ‐Tyson theory to predict the ﬂexural strength of CNT ‐cement nanocomposites. The proposed model considers the effects of experimental variables on CNT dis- persion quality, bonding mechanism and the ﬂexural strength. To this end, a Bayesian methodology is employed to calibrate the unknown model parameters and various sources of uncertainty using extensive testdata. The model is then used to identify the optimum ranges of variables to maximize the ﬂexural strength through computing the failure probability which is de ﬁned as the probability of not meeting certain strength requirements (herein, 50% increase compared with the control). The model suggests that CNT aspect ratioranges from 400 to 800 and concentration between 0.08 and 0.18 c‐wt% yields the highest ﬂexural strength. Finally, the effect of changes in experimental variables on the probability estimates is examined using sensitiv- ity and importance measures. The analysis reveals that the proposed model can capture the experimentallyobserved trends with reasonable accuracy. For example, the importance of age increases as CNT concentrationincreases. 1. Introduction Concrete cracks when subjected to excessive tensile forces. Over the past decades, various types of ﬁbers (steel or synthetic) have been used to compensate for the low tensile/ ﬂexural strength of concrete by mitigating the crack propagation [1–4]. The contributions of ﬁbers to delay the crack propagation rely on the shape and size of the ﬁbers, their surface textures and the interfacial bond strength between the ﬁbers and the matrix [5,6] . In this regard, incorporating smaller ﬁbers might help to delay the initiation of cracks [7]. Recently, carbon nan- otubes (CNTs) have been extensively used to increase the ﬂexural strength of cementitious composites (cement paste, cement mortar and concrete) [8,9] . The nanoscale size of CNTs may reduce the dis- tance between adjacent nanotubes, delaying the crack propagation. Besides, CNTs have tensile strength of around 30 GPa with elastic modulus in the range of 1 TPa[10,11] . Because of the extraordinary physical and mechanical properties of CNTs, many studies have reported that the ﬂexural strength of CNT reinforced cementitious composites has been increased signi ﬁcantly [8,9,12,13] . However,the authors ’statistical analysis reported inconsistent results in the ﬂex- ural strength of cementitious materials reinforced with CNTs (herein, CNT ‐cement composites) [7]. The discrepancy can be attributed to the CNT dispersion and bonding states within the cement matrix[14–16]. The bond strength between CNTs and the cement matrix is governed by multiple factors including CNT type, surface condition and dispersion quality [17–19]. Through performing extensive statisti- cal analysis of experimental test data reported in the literature, the authors found that CNT dispersion quality signi ﬁcantly contributed to the mechanical properties [7]. The dispersion quality, itself, depends on multiple factors such as physical characteristics of CNTs, dispersion procedure and cement matrix composition. The improper combinations of these factors often lead to degradation of the interfa- cial bond and ﬂexural strengths. These phenomena are explained in the authors' previous works [7,17,20] . To address these phenomena, an analytical model can be a viable option to estimate the bond strength between CNTs and the matrix considering the contributions of multiple experimental variables for maximizing the ﬂexural strength of CNT ‐cement composites. https://doi.org/10.1016/j.compstruct.2020.112748 Received 27 April 2020; Revised 23 July 2020; Accepted 26 July 2020 Available online 30 July 20200263-8223/© 2020 Elsevier Ltd. All rights reserved.⇑Corresponding author. E-mail addresses: mahyar.ramezani@louisville.edu (M. Ramezani), young.kim@louisville.edu (Y.H. Kim), z.sun@louisville.edu (Z. Sun).Composite Structures 253 (2020) 112748 Contents lists available at ScienceDirect Composite Structures journal homepage: www.elsevi er.com/locate/compstruct Recently, various analytical models considering CNT physical char- acteristics (e.g., length, diameter and volume fraction) have been used to predict the tensile strength of polymers, metals and ceramics [21,22] . For example, the interfacial bond strength of CNTs in a metal matrix was increased when using longer CNTs, resulting in higher ten- sile strength [21]. The agglomeration of CNTs, however, reduced the frictional forces resulting in debonding from the metallic matrix [23]. Although the interfacial bond strength is very critical in deter- mining the tensile/ ﬂexural strength [24], the existing models (e.g., Halpin ‐Tsai model [25] and Kelly ‐Tyson model [26]) typically assume a homogenous dispersion and perfect bond between CNTs and the hosting matrix [27]. Therefore, few researchers have modi ﬁed the existing models by adding a multiplier of dispersion coef ﬁcient [28–30]. The different dispersion mechanisms of CNTs in cement matrix, however, limits the applicability of the existing models to CNT ‐cement composites. The dispersion coef ﬁcient should therefore be formulated by parameters such as the CNT physical characteristics, dispersion procedure and matrix composition. Recently, the literature data was analyzed to identify the most important variables and their interactions affecting the dispersion quality, bond strength and the ﬂexural strength [7]. Besides, the in ﬂuence of CNTs on the rate of cement hydration is another important variable affecting the ﬂexural strength that needs to be considered in the prediction model [7,17,31] . This study aims to propose a reliable model in a probabilistic man- ner to understand the effects of CNT dispersion, bond strength and cement hydration process on the ﬂexural strength of cementitious composites (herein, cement pastes and mortars). Note that some researchers have recently started to incorporate CNTs within concrete [32,33] . The proposed probabilistic model can be updated to include the in ﬂuence of CNTs on the ﬂexural strength of concrete (i.e., with coarse aggregate), if more data becomes available in the future. The following section presents the base model (Kelly ‐Tyson model [26]) and the required modi ﬁcations in order to improve its predictability for CNT ‐cement composites. Besides, because incorporating CNTs is correlated with a high degree of uncertainty, a probabilistic approach is adopted to consider various sources of uncertainty to affect the ﬂex- ural strength. Using existing experimental test data obtained from the literature, several probabilistic models are developed and the best model is selected. Then, the performance of the proposed model is evaluated using the reliability analysis. The proposed model in this study can explain the reported experimental trends in the literature including the contradictory test results and provide valuable informa- tion for future researchers and designers to identify the best combina- tion of the experimental variables to achieve certain ﬂexural performance. 2. Modi ﬁed Kelly-Tyson model The Kelly ‐Tyson model [26] which is the extension of the well ‐ known rule of mixtures has been extensively used to predict the tensile strength of brittle and ductile materials containing CNTs such as poly- mers, metals and ceramics [22,34 –37]. Besides, the Kelly ‐Tyson model was used to predict the ﬂexural strength of concrete reinforced with short discrete ﬁbers [38,39] . For condition of equal elastic strain in theﬁber and matrix, the rule of mixtures is considered as the simplest way to model the tensile strength of composites reinforced with aligned continuous ﬁbers [40]: σc¼σfvfþσmð1/C0vfÞð 1Þ where σcis the ultimate tensile strength of the composite ( MPa orpsi), σfis the tensile strength of the ﬁber ( MPa orpsi),σmis the tensile strength of the matrix at ultimate strain of the composite ( MPa orpsi) and vfis the volume fraction of the ﬁber (unitless). To consider the con- tribution of the semi ‐aligned short discrete ﬁbers, Eq. (1)can be mod- iﬁed to include the ﬁber length ( χL) and orientation ( χo) factors. Toaccount for the randomness in ﬁber orientations, χocan be calculated asχo¼∑nancos4Φnwhere anis the fraction of ﬁbers having the orien- tation angle Φnwith respect to the reference axis [41]. When using aligned ﬁbers, the value of χois one and it is equal to 3/8 and 1/6 for randomly orientated ﬁbers in two and three dimensions, respec- tively [42]. Due to the extremely small sizes of CNTs, they are assumed to be randomly arranged in three dimensions (i.e., χo¼1=6)[28,29] . Kelly and Tyson [26] have developed the χLforﬁbers having lengths ( L) shorter or longer than the ﬁber critical length (Lc¼σfd=2τi;dis the ﬁber diameter and τiis the interfacial bond stress) to determine the mode of failure: ﬁber fracture ( L≥Lc)o r pull ‐out ( L<Lc). Therefore, the tensile strength of composites can be calculated as follows: σc¼χo∑L<Lc L¼0τiL d/C18/C19 þ∑L L¼Lc1/C0Lc 2L/C18/C19 σf"# vfþσmð1/C0vfÞð 2Þ where the ﬁrst and second terms in the bracket are contributions from ﬁbers with L<Lcand L≥Lc, respectively. Because CNT length in cementitious composites is generally smaller than the critical length (L<Lc), the contribution of CNTs with L≥Lcis neglected. Also, due to the very low volume fraction of CNTs used in cemen- titious materials ( vf¼0.00004 –0.00128; based on the database), the predicted ﬂexural strength of CNT ‐cement nanocomposites using Eq. (2)is almost the same as the ﬂexural strength of the matrix. The authors ’preliminary study has shown that vf2=3is appropriate to ﬁt the experimental data. The modi ﬁed form of the Kelly ‐Tyson equation to predict the ﬂexural strength of CNT ‐reinforced cementitious com- posites can therefore be expressed as follows: σc¼1=6τiðARÞvf2=3þσm1/C0vf2=3/C0/C1 ð3Þ where ARis CNT aspect ratio ( AR¼L=d). To predict the ﬂexural strength using Eq. (3),τimust be quanti ﬁed. The experimental investi- gations to measure τiis challenging or impractical due to the limitations associated with nanoscale sizes. Therefore, theoretical approaches might be used to determine the interface properties [43]. In this study, the Pukanszky model [44] is used to correlate the composite strength with the interface properties. Note that Pukanszky model have been widely used to characterize the composites containing various nanopar- ticles including nanoclays, layered silica nanoparticles and CNTs [22,45,46] : σc¼σmð1/C0vf2=3 1þ2:5vf2=3ÞeBvf2=3ð4Þ where Bis the empirical adhesion parameter that re ﬂects the capacity of stress transfer between CNTs and matrix. In this study, Bis calculated using the available experimental data from the literature including the authors ’work. According to Lazzeri and Phuong [47], the adhesion parameter Bcan be correlated with τiin the Kelly ‐Tyson model. There- fore, in this study, τiis attained by connections between Eqs. (3) and (4): τi¼6σmðB/C02:5Þ ARð5Þ To calculate τifor each individual data, the average value of Bfor the entire dataset ( Bavg¼55) is used as the representative value for adhesion between CNTs and cement matrix. The statistics of Band τi are shown in Table 1 . 2.1. Critical factors affecting CNT mechanisms in cementitious composites To understand the mechanisms of CNTs within cementitious mate- rials (e.g., bonding, dispersion and cement hydration rate mecha- nisms), three critical factors were proposed in the authors ’previous work [17] to capture the interactions between the important experi- mental variables to predict the mechanical properties including the ﬂexural strength.M. Ramezani et al. Composite Structures 253 (2020) 112748 2 1) dispersion factor ( ηD) to correlate CNT dispersion as a functions ofAR, CNT volume fraction in aqueous solution ( vf/C0w), superplasticizer ‐to‐CNTs ( SP/CNTs ) ratio and Ultrasonication Energy Indicator ( UEI) 2) matrix factor ( ηM) to include the in ﬂuence of water ‐to‐cement ratio ( w/c) and sand ‐to‐cement ratio ( s/c) 3) hydration age factor ( ηH) to consider the interactions between CNT concentration ( κ) and hydration age of cement ( t) A brief description of these critical factors is provided in Table 2 . More details are provided in the authors ’previous study [17]. The applicability of the proposed critical factors to predict the ﬂexural strength is validated in Section 5.1 . To correct the bias in Eq. (3)and to capture the mechanisms of CNTs in cementitious materials, the proposed critical factors ( ηD,ηH and ηM) are added to the modi ﬁed form of the Kelly ‐Tyson model as follows: σc¼ηDηHηM1=6τiðARÞvf2=3þσm1/C0vf2=3/C0/C1 /C2/C3 ð6Þ In addition, Eq. (6)is formulated in a probabilistic manner using a Bayesian methodology to consider different sources of uncertainty such as the statistical uncertainty in the estimation of the unknown model parameters, the model errors associated with the inaccurate form of the model and missing variables [48] that might in ﬂuence theﬂexural strength. 3. Formulation of the probabilistic model To formulate the modi ﬁed form of the deterministic Kelly ‐Tyson model (see Eq. (6)) in a probabilistic manner, additional correction terms (herein, ηD,ηHand ηM) and model error are added to the deter- ministic expression to account for the inherent bias and to capture the remaining variability in residuals, respectively [48]. The probabilistic model relates the quantity of interest (herein, the ﬂexural strength; σc) to a set of measurable variables (de ﬁned as x¼ðx1;x2;x3;/C1/C1/C1 Þ) such as AR,κ,UEI,SP/CNTs ratio, w/cratio, s/cratio, t, etc. The probabilistic capacity model can be written as follows: Cx;ΘðÞ ¼ CxðÞ þ γx;θðÞ þ @ɛ ð7Þwhere Cx;ΘðÞ is the capacity model of interest (herein, σc),xis a set of measurable variables, Θ¼θ;σðÞ ;θ¼ðθ1;θ2;θ3;/C1/C1/C1 Þdenotes a set of unknown model parameters used to ﬁt the data, CxðÞis the selected deterministic model, γx;θðÞ is the correction term (herein, ηD,ηHand ηM),@is a standard deviation of the error term taken from VarCx;ΘðÞ½/C138 ¼ @2,ɛis a random variable with zero mean and unit vari- ance and @ɛis the error term of the model. The ﬁnal form of the prob- abilistic model used to predict the ﬂexural strength of CNT ‐cement composites is as follows: σc¼ηDxD;θD ðÞ /C2 ηHxH;θH ðÞ /C2 ηMxM;θM ðÞ /C21=6τiARðÞ vf/C0m2=3þσm1/C0vf/C0m2=3/C0/C1 /C2/C3 þ@ɛ ð8Þ (8)where vf/C0mis the volume fraction of CNTs within the cement matrix and θ¼ðθD;θH;θMÞare unknown model parameters used to ﬁt the data: ηDxD;θD ðÞ ¼ 1/C0ARθ1/C2vf/C0wθ3 UEIθ2/C2SP=CNTs ðÞθ4/C16/C17hi , ηHxH;θH ðÞ ¼t 4þ0:85t/C16/C17κθ5and ηMxM;θM ðÞ ¼1þs=c ðÞθ6 w=cðÞθ7: 3.1. Bayesian parameter estimation To calibrate the unknown model parameters, a Bayesian updating rule [49] is used because it allows to update their statistics (i.e., means, variances and correlation coef ﬁcient) as new data becomes available in the future [50]: fΘðÞ ¼ λLðΘÞpðΘÞð 9Þ where Θis a vector of unknown model parameters, pðΘÞis prior distri- bution of Θbased on the knowledge about Θthat represents the infor- mation available before collecting the data Đof size n,LðΘÞis the likelihood function that captures the information on Θfrom data Đ which is proportional to the conditional probability PDjΘðÞ of observ- ingĐfor given values of Θ,λ¼RLðΘÞpðΘÞdðΘÞ/C2/C3/C01is a normalizing factor and fΘðÞis posterior distribution that shows the updated state of information about Θ. Once fΘðÞis known, the posterior mean vector (MΘ) and the covariance matrix ( ∑ ΘΘ) are obtained. Because there is no prior information about Θis available before collecting data Đ,a non ‐ informative prior distribution is assumed [48]. If new data becomes Table 2 Critical factors affecting the ﬂexural strength. Critical factorEquation Variable (unit) Mechanism ηD 1/C0AR UEI/C2vf/C0w ðSP=CNTs ÞhiAR¼CNT aspect ratio (unitless) vf/C0w¼CNT volume fraction in aqueous solution (unitless) UEI¼ultrasonication energy indicator (index, express in minute)SP/CNTs ¼superplasticizer-to-CNTs ratio (unitless)Ifv f/C0worAR¼0 (i.e., no CNTs added), ηD¼1 (i.e., perfect dispersion). Asvf/C0wand/or ARincreases, ηDgets smaller. AsUEIand/or SP/CNTs increases, ηDgradually reaches one. Beyond certain UEIand/or SP/CNTs ,ηDgets constant. ηM1þs=c ðÞ w=cðÞs=c¼sand-to-cemet ratio (unitless) w=c¼water-to-cement ratio (unitless)Asw/cratio decreases, the value of ηMincreases. Ass/cratio increases, the value of ηMincreases. ηH t 4þ0:85t/C16/C17κt¼age of specimen (days) κ¼CNT concentration ( c-wt%)If CNT is not added ( κ=0 c-wt%), ηH¼1. At low CNT concentrations ( κ≤0.1c-wt%), the in ﬂuence of tonηHis minimal.Asκincreases, η Hreaches unity ( ηH¼1) at later ages (i.e., the higher value oft). Note: ηD: dispersion factor; ηM: matrix factor; ηH: hydration age factor.Table 1 Statistics of the Pukanszky adhesion parameter ( B) and interfacial bond stress ( τi). Parameter Minimum Maximum Mean Median Standard deviation B(unitless) -7.827 255.9 54.92 44.19 43.85 τi(MPa) 0.63 9.28 2.98 1.97 2.27M. Ramezani et al. Composite Structures 253 (2020) 112748 3 available in the future, Eq. (10)can still be used to update the estimates of the unknown model parameters by using the current fΘðÞas a new previous distribution. 4. Database A total of 105 data is collected from several literatures including the authors ’work [8,12,17,20,51 –56]to calibrate the unknown model parameters in Eq. (9)using a Bayesian methodology. To evaluate the robustness of the proposed probabilistic model, the database is split into two sets of data (training set and test set) using a hold ‐out method [57]. This process minimizes the issues of over ﬁtting and being overly optimistic on model error [58]. The training set is used to estimate the unknown model parameters and the test set is used to evaluate the per- formance of the developed model. The training set is selected using a random subsampling method without replacement [57]. Besides, two different random selections of AandBare provided. The two different random selections is to ensure that more data get a chance to be the test set, resulting in an unbiased evaluation of the model performance. In both random selections, 80% of the data (training set; 84 data points) is used to develop the model (datasets I and III for random selections Aand B, respectively) and the remaining 20% of the data (test set; 21 data points) is used to validate the model (datasets IIand IV for random selections Aand B, respectively). Table 3 shows the minimum, maximum, mean, median and stan- dard deviation of different variables used to develop the probabilistic model (random selection A: dataset I). When the ranges of data in Table 3 are considered in the application of the suggested model, the model can predict the ﬂexural strength within a certain con ﬁdence level. In this study, the average values of L,dand ARare used. For example, if Land dare reported to be between 10 and 30 μmand 20–40nm, respectively, then average length, diameter, and aspect ratio of 20 μm,3 0 nmand 667 (unitless) are used. 5. Probabilistic model assessment Tables 4 and 5 show the posterior statistics of Θestimated using dataset I (for random selection A) and dataset III (for random selection B), respectively. The mean values of the unknown model parameters for datasets I and III differs from one another by less than 20% (except forθ 5Aand θ5Bwhich exhibit 22.1% difference). This indicates the robustness of the proposed model. Besides, the accuracy of the pro- posed probabilistic model is evaluated as shown in Fig. 1 . The solid cir- cles are the datasets I and III used to develop the model and the open triangles are the datasets II and IV used to validate the model. The pre- dicted ﬂexural strength for most of the data shows good accuracy within the solid lines (i.e., within ± standard deviation of the error term in the proposed model; @A= 15% and @B= 18%). The robust- ness of the proposed probabilistic model can be validated using data-sets II and IV (open triangles) that are not used for the estimation of the unknown model parameters. Note that the three data points inside the oval in Fig. 1 (a) and (b) are related to the ﬂexural strength of cement mortar reinforced with mechanical functionalized CNTs at dif- ferent ages (3, 7 and 28 days). The mechanical functionalized CNTs provided superior bond properties [9], resulted in under ‐prediction of the ﬂexural strength. 5.1. Applicability of the proposed critical factors In this section, the corrected critical factors (i.e., ηD(xD,θD),ηH(xH, θH) and ηM(xM,θM); see Eq. (8)) are further analyzed in order to val- idate their applicability to predict the ﬂexural strength without bias. Note that because the mean values of the calibrated unknown model parameters for random selections Aand Bare close to one another (seeTables 4 and 5 ); the statistics of the unknown model parameters for random selection Ais used for the analysis throughout this research. The proposed probabilistic model after adding the mean val- ues of the unknown model parameters using random selection Ais as follows: σc¼1/C0AR0:39/C2vf/C0w0:48 UEI0:28/C2SP=CNTs ðÞ0:13 !"# /C2t 4þ0:85t/C18/C191:06κ /C21þs=c ðÞ0:17 w=cðÞ0:13/C21=6τiARðÞ vf/C0m2=3þσm1/C0vf/C0m2=3/C0/C1 /C2/C3 þ0:15ɛ ð10Þ The corrected critical factors are analyzed by comparing the esti- mated with the measured values. The estimated value is obtained by plugging in the value of the experimental variables in corresponding relations. For example, the estimated value of ηM(xM,θM) is directly calculated using ηMxM;θM ðÞ ¼1þs=c ðÞ0:17 w=cðÞ0:13hi . The measured value of each critical factor is back calculated using Eq. (10) and the experimental test results of σcand σm, assuming that ηD(xD,θD),ηH(xH,θH) and ηM(xM,θM) are not correlated. Figs. 2 through 4 are plotted to validate the applicability of ηD(xD, θD),ηH(xH,θH) and ηM(xM,θM), respectively. Each critical factor is plotted as a function of a main variable ( x‐axis) and its interaction with another independent variable at two different ranges: R1 for the exper- imental variable below its mean value and R2 for the experimental variable above the mean value (see Table 6 ). For example, Fig. 2 (a) shows the in ﬂuence of changes in ARand its interaction with vf/C0wat two different ranges of R1(vf/C0w≤0:0012) and R2(vf/C0w>0:0012) onηD(xD,θD). The solid lines in each plot represent the estimated val- ues of ηs. To plot the solid lines, the mean value of other variables is used except for the main variable. Also, the dashed lines represent the standard deviation of the model error (± @A). The open circles rep- resent the measured values of ηs. Table 3 Range of Variables to Develop Probabilistic Model (random selection A: dataset I). Variables Minimum Maximum Mean Median Standard deviation κ(c-wt%) 0.0125 0.5 0.1033 0.09 0.0987 L(μm) 1.5 55 21.7 20 13.45 d(nm) 8 32.5 25.29 27.5 7.58 AR(unitless) 157.89 2500 902.83 800 622.28 SP/CNTs 1.32 24 6.4 4 5.26 w/c 0.3 0.6 0.426 0.45 0.0748 s/c 0 3 1.527 2.25 1.235 Age(days) 3 28 14.643 7 11.18 σm(MPa) 1.94 9.3 5.36 4.88 1.91 UEt(J/ml) 25 2800 1703.65 2612.5 1208.2 UEm(J/ml/ min) 3.95 28 9.21 8 5.998 UA(%) 30 100 55.63 50 15.63M. Ramezani et al. Composite Structures 253 (2020) 112748 4 Figs. 2 through 4 exhibit that most of the measured ηs are within ±@Aand they follow the similar trends as the estimated ηs. The good correlations between the measured and estimated values validate the applicability of the proposed critical factors to capture the mechanisms of CNTs affecting the ﬂexural strength. As an example, comparing R1 andR2i nFig. 3 , it is obvious that the in ﬂuence of ton the value of ηH (xH,θH) becomes more pronounced at higher κ(see R2i nFig. 3 ). 5.2. Model selection criteria For model evaluation, it is assumed that there is a single correct model (or at least best model) that tradeoffs between goodness of ﬁt and parsimony (i.e., preserves the model accuracy while using less variables). In this study, the two most commonly used penalized like- lihood information criteria, Akaike information criterion ( AIC) andBayesian information criterion ( BIC), are used to select the best model. The AICand BICcan be calculated as follows: AIC¼/C02logLΘðÞ þ 2Np ð11Þ BIC¼/C02logLΘðÞ þ NplogN s ð12Þ where logLΘðÞis a measure of model ﬁt,Npis the number of unknown model parameters included in ΘandNsis the number of samples. The best model is selected from a set of candidate models containing differ- ent subsets of the critical factors as listed in Table 7 . The best model is the one that neither under ‐ﬁts nor over ‐ﬁts (i.e., the model with the smallest approximated values of AICandBIC). As seen in Table 7 , Both AICand BICestimated values are smallest for model VIII compared with other competing models, indicating that including all critical fac- tors is best to accurately predict the ﬂexural strength.Table 4 Posterior Statistics of the Unknown Model Parameters (random selection A: dataset I). Unknown parameter Mean Standard deviation Correlation matrix θ1A θ2A θ3A θ4A θ5A θ6A θ7A @A θ1A 0.3925 0.0280 1 −0.5390 0.9981 0.9495 0.8606 −0.5683 −0.9771 0.6282 θ2A 0.2798 0.0086 1 −0.5713 −0.6965 −0.6583 −0.0396 0.6082 −0.4281 θ3A 0.4806 0.0299 1 0.9483 0.8541 −0.5671 −0.9786 0.6275 θ4A 0.1317 0.0369 1 0.9458 −0.3876 −0.9627 0.6753 θ5A 1.0589 0.0575 1 −0.2129 −0.8844 0.5981 θ6A 0.1730 0.0046 Sym. 1 0.4226 −0.5294 θ7A 0.1269 0.0146 1 −0.6007 ∂A 0.1491 0.0008 1 Table 5 Posterior Statistics of the Unknown Model Parameters (random selection B: dataset III). Unknown parameter Mean Standard deviation Correlation matrix θ1B θ2B θ3B θ4B θ5B θ6B θ7B @B θ1B 0.3950 0.0876 1 0.5287 0.9739 0.4178 −0.4650 −0.4585 −0.7044 0.0411 θ2B 0.3350 0.0745 1 0.4056 −0.0169 −0.1806 −0.2618 −0.5600 −0.0174 θ3B 0.4751 0.0926 1 0.3769 −0.4299 −0.4915 −0.7011 0.0524 θ4B 0.1055 0.0341 1 −0.2875 −0.1602 −0.4694 0.1143 θ5B 0.8846 0.2263 1 0.1879 0.3906 −0.0287 θ6B 0.1697 0.0255 Sym. 1 0.1756 −0.0573 θ7B 0.1379 0.0555 1 −0.0903 ∂B 0.1768 0.0137 1 Fig. 1. Predicted value of ﬂexural strength versus experimental results (a) random selection A(b) random selection B(datasets I & III are used for model development and datasets II & IV are used for model validation).M. Ramezani et al. Composite Structures 253 (2020) 112748 5 In addition, ΔAICscore which is de ﬁned as the difference between AICvalues of the best model ( AIC min) and other competing models (AIC i) are calculated: ΔAIC¼AIC i/C0AIC min ð13ÞAccording to Burnham and Anderson [59], the models with ΔAIC≤2 have substantial evidence of better performance. The models with 4 ≤ΔAIC≤7 have less evidence and those models having ΔAIC> 10 have no evidence [59]. It can be seen that all competing Fig. 2. Estimated vs. measured ηD(xD,θD) as a function of: (a) ARand its interaction with vf/C0w(b)UEIand its interaction with SP/CNTs ratio. Fig. 3. Estimated vs. measured ηH(xH,θH) as a function of tand its interaction with κ.M. Ramezani et al. Composite Structures 253 (2020) 112748 6 models have ΔAICvalue of more than 10. This further con ﬁrm the superior performance of model VIII (Eq. (9)) compared with other can- didate models. 6. Probability of not meeting design speci ﬁcations After calibrating the unknown model parameters using dataset I of random selection A(seeTable 4 ), the probabilistic model (Eq. (10))i s used to identify the optimum ranges of each variable to maximize the ﬂexural strength. This could be achieved by computing the probability of not meeting a speci ﬁedﬂexural strength requirement (probability of failure). Every structure has a capacity ( C) and is subjected to some sort of demand ( D) where both of them depend on random variables. There- fore, each of Cand Dhas a probability distribution ( fDDðÞand fCðCÞ) which could be combined to form a joint probability density function(f CDðC;DÞ; see Fig. 5 (a)). If the joint probability density function is known, the probability of failure ( Pf) could be easily determined by Pf¼R1 /C01RC≤D /C01fCDC;DðÞ dCdD . However, because of the large number of variables, the probability of failure cannot be determined in a math- ematical way and therefore might be estimated using simulation ‐based analyses. Following the conventional notation in structural reliabilitytheory [60], a limit state function gð/C1Þis introduced such that the event [gð/C1Þ≤0] denotes the conditional probability of Cnot meeting D(i.e., C/C0D<0; failure occurs). Conversely, if C/C0D>0, the structure is safe. In this study, because there is no design criteria available in the lit- erature, the minimum ﬂexural strength requirement of CNT ‐cement composites is set to the deterministic value of 50% increase compared with the control specimen (i.e., D= 1.5 σm). Using the proposed prob- abilistic model in Eq. (9), the limit state function is formulated: gðασm;x;ΘÞ¼σcx;ΘðÞ /C0 ασm ð14Þ where αis the minimum ﬂexural strength requirement (i.e., α= 1.5). Therefore, the probability of not meeting the demand is written as follows: PΘðÞ ¼ P½gðασm;x;ΘÞ≤0jΘ/C138ð 15Þ where P/C1j/C1ðÞrepresents the conditional probability of gασm;x;Θ ðÞ ≤0 for given values of Θ. Fig. 5 (b) visualizes the probability of failure using the isoprobabil- ity lines which connect pairs of values of different variables corre- sponding to the same probability of failure ( p½gðασm;x;ΘÞ≤0/C138). The lower values of probability of failure ( Pf) indicates that the likelihood ofC/C0D>0 is higher (i.e., the ﬂexural strength meets the design criteria). 6.1. Relationship between variables involved in ultrasonication procedure Various dispersion techniques have been employed to uniformly disperse CNTs within the cement matrix; amongst which an ultrason- ication procedure is the most common practice [7,61] . Fig. 6 (a), (b) and (c) shows the contour plots of isoprobability lines for the ultrasonication energy indicator ( UEI) versus κ,ARand SP/CNTs ratio, respectively, using the mean value of other variables.Table 6 Ranges used for each Critical Factor. Critical factor (y-axis)Main variable (x-axis)Range (based on mean values of data) R1 R2 ηDðxD;θDÞ AR vf/C0w≤0.0012 vf/C0w> 0.0012 ηDðxD;θDÞ UEI SP/CNTs ≤6.4 SP/CNTs > 6.4 ηHðxH;θHÞ t κ≤0.1033 κ> 0.1033 ηMðxM;θMÞ s/c ratio w/c ≤0.426 w/c> 0.426 Table 7 Model Selection Criteria (using random selection A: dataset I). Model Np Information criteria R-squared # Critical factors AIC BIC ΔAIC I (Eq. (3)) Deterministic model 0 −18.73 −18.73 25.02 0.79 II ηD 4 −8.014 −8.317 31.40 0.84 III ηM 2 −11.15 −11.30 28.26 0.86 IV ηH 1 −17.58 −17.66 21.83 0.87 V ηD×ηM 6 −6.434 −6.889 32.98 0.86 VI ηD×ηH 5 −11.80 −12.18 27.61 0.86 VII ηM×ηH 3 −16.43 −16.66 22.98 0.87 VIII (Eq. (8)) ηD×ηM×ηH 7 −39.41 −39.94 0 0.95 Fig. 4. Estimated vs. measured ηM(xM,θM) as a function of s/cratio and its interaction with w/cratio.M. Ramezani et al. Composite Structures 253 (2020) 112748 7 Fig. 6 (a) shows that when using κ≤0.05 c‐wt% (based on weight percent of cement powder), UEIhas negligible impact on decreasing Pf. However, in case of a higher κ, the impact of UEIonPfbecomes more pronounced. For example, in case of κ= 0.1 c‐wt%, increasing UEIfrom 5 to 25 results in 25% decrease in Pffrom 0.85 to 0.64 (see crosses on the dashed line in Fig. 6 (a)). In addition, Pfis lowest when incorporating UEIof above 25 and κbetween 0.10 and 0.18 c‐wt% (see the hatched area in Fig. 6 (a)). Fig. 6 (b) shows as ARincreases, the contribution of UEItoPf becomes more pronounced. For example, in case of AR= 200, increas- ingUEIfrom 10 to 30 does not contribute to attaining lower Pf(see solid circles on the same probability of failure line ( Pf¼0.96) in Fig. 6 (b)). However, in case of AR= 500, Pfdecreases from 0.92 to 0.77 by increasing in UEIfrom 5 to 25 (see crosses on the dashed line in Fig. 6 (b)). Also, the lowest Pfis achieved when incorporating UEIof above 25 and ARof around 500 (see the hatched area in Fig. 6 (b)). Generally, as UEIincrease, the circles which represent the iso- probability line of failure becomes smaller. This suggests that UEIis reaching to its threshold value. Fig. 6 (c) shows that the impact of UEIon the Pfincreases as SP/ CNTs ratio increases. Generally, as SP/CNTs ratio increases, lower UEIis needed to attain the same Pf. For example, incorporating UEI and SP/CNTs ratio of 30 and 6, respectively, exhibits the same Pfas utilizing UEIandSP/CNTs ratio of 25 and 7 or 20 and 12, respectively, does (see the solid squares on the same probability of failure line ( Pf¼0.63) in Fig. 6 (c)). Besides, in case of a constant UEI, the probability of failure is not changed by increasing in SP/CNTs ratio beyond a certain limit. For example, in case of UEI= 10, Pfreaches a constant value of around 0.7 beyond SP/CNTs ratio of approximately 7 (see solid circles on the dashed line in Fig. 6 (c)) and this threshold reaches sooner in lower UEIvalues. This might be explained by the better dispersibility of CNTs using higher UEI. Thereafter, higher amount of SPcould be taken by CNTs to wet their surface [62]. Conversely, in the case of SP/CNTs ratio of below 3, higher amount of UEI does not con- tribute to lowering the Pf(see crosses on the same Pfline of 0.96 in Fig. 6 (c)). 6.2. Relationship between w/c and s/c ratios Fig. 7 shows the contour plot of the probability of failure lines as a function of s/cand w/cratios. The probability of failure decreases as s=cratio increases. For example, in case of a constant w/c= 0.45, increasing s/cratio from 0 (i.e., cement paste) to 3 results in decreas- ing the Pffrom 0.93 to 0.84 (see the solid squares on the dash ‐dotted line in Fig. 7 ). Conversely, Pfincreases as w/cratio increases. For example, in case of a constant s/cratio of 2, Pfincreases from 0.79 to 0.95 as the w/cratio increases from 0.35 to 0.60 (see the solid cir- cles on the dashed line in Fig. 7 ). Also, there are several options to attain the same Pfusing certain combinations of w/cand s/cratios. For example, Pf= 0.80 can be attained when using w/cands/cratios Fig. 5. Visualization of probability of failure: (a) joint probability density function (b) isoprobability lines. Fig. 6. Contour plots of the probability of failure as functions of (a) UEIand CNT concentration, (b) UEIand CNT aspect ratio (c) UEIand SP/CNTs ratio.M. Ramezani et al. Composite Structures 253 (2020) 112748 8 of 0.40 and 3 or 0.35 and 2, respectively (see crosses on the same Pf line of 0.8 in Fig. 7 ). 6.3. Relationship between κand t Fig. 8 shows the in ﬂuences of κand tonPf. In case of κ≤0.05 c‐wt%, the in ﬂuence of tonPfis negligible. For example, when κ= 0.05 c‐wt% (see dashed line in Fig. 8 ),Pfis not affected by the variation in tfrom 10 to 28 days (see the solid circles on the same probability of failure line of Pf¼0.78 in Fig. 8 ). The contribution of ttoPfincreases as κincreases. The same Pfcan be achieved using cer- tain combinations of κand t. For example, when κ= 0.1 c‐wt%,Pf reaches 0.66 at 14 days, while 25 days is needed to attain the same Pfatκ= 0.4 c‐wt% (see crosses on the same probability of failure line inFig. 8 ).6.4. In ﬂuence of optimum ranges of variables on probability of failure Fig. 9 shows contour plots of the probability of failure as a function ofκandAR.Fig. 9 (a) and (b) show the isoprobability lines for ﬂexural strength requirement of at least 50% increase compared with the con- trol (i.e., D= 1.5 σm) when using mean and optimum values of other variables, respectively. The optimum values of variables ( UEI= 30, SP/CNTs =6 , w/c= 0.4 and s/c= 3) is obtained using the optimal PfinFig. 6 (c) and 7. Generally, Pfdecreases as κorARincreases up to their optimum values beyond which Pfincreases. Also, In case of κ≤0.05 c‐wt%, the impact of ARonPfis minimal. Besides, in case of AR≤500, the impact of κonPfis negligible. Fig. 9 (a) shows that Pfis lowest ( Pf¼0.84) when using ARranging from 400 to 800 and κbetween 0.08 and 0.18 c‐wt% (see hatched areas in Fig. 9 (a)). Comparing Fig. 9 (a) and (b), it is seen that incor- porating optimum values of other variables signi ﬁcantly contribute to decreasing Pfby 64% from 0.84 to 0.30. Besides, the same Pfcan be attained when using CNTs having either lower κwith higher ARor higher κwith lower AR. For example, the probability of failure is 0.3 in case of using AR= 800 with κ= 0.1 c‐wt%o r AR= 500 with κ= 0.2 c‐wt% (see crosses on the same Pfline of 0.3 in Fig. 9 (b)). The similar trend was also observed in the experimental studies [12,19] . 7. Sensitivity and importance analyses This section presents the sensitivity and importance analyses for different variables in the limit state function (Eq. (14)) to assess their effect on the probability estimates, helping researchers to minimize Pf by designing more ef ﬁciently. 7.1. Sensitivity analysis The sensitivity analysis is used to determine which variables are most sensitive to the change of Pf[63]. To compute the sensitivity measures of different variables, the vector δis used: δ¼σ∇$Θfβ ð16Þ where σis diagonal matrix with diagonal elements given by the stan- dard deviation of each random variable in x,βis a reliability index from ﬁrst order reliability method (FORM) analysis, and ∇$Θfβis a gradient vector of βwith respect to Θfcomputed at the mean point. The sensi- tivity measures can be employed for optimal design and to gain insights into which variable(s) to be changed in order to increase the ﬂexural strength. The sensitivity measure of each variable is computed with respect to three main variables: κ,ARand UEI. These variables are believed to have the highest impact on the ﬂexural strength of CNT ‐cement nanocomposites. Fig. 10 shows the sensitivity measure of each variable as a function of κ(see Fig. 10 (a)), AR(see Fig. 10 (b)) and UEI (see Fig. 10 (c)) while incorporating mean values of other variables. The higher change in the sensitivity measure of each variable ( y‐ axis) with respect to the changes of x‐axis (i.e., main variable) indi- cates the higher sensitivity of Pfto the change in that variable (for D= 1.5 σm). Fig. 10 (a) shows that up to κ= 0.2 c‐wt%, the sensitivity measure of most variables changes signi ﬁcantly. However, the sensitivity mea- sure remains almost constant beyond κ= 0.2 c‐wt%. In case of κ≤0.2 c‐wt%,w/cratio, ARand s/cratio are the most sensitive variables to the changes in Pf, respectively, while SP/CNTs ratio and tare the least sensitive variables. In case of κ> 0.2 c‐wt%,Pfis most sensitive to the variation in tfollowed by AR. When computing the sensitivity mea- sures as a function of AR(seeFig. 10 (b)), κis the most sensitive vari- able below ARof around 500 beyond which the sensitivity of κto Fig. 7. Contour plot of probability of failure as a function of w/cand s/c ratios. Fig. 8. Contour plot of probability of failure as a function of CNT concen- tration and hydration age.M. Ramezani et al. Composite Structures 253 (2020) 112748 9 change the Pfgradually decreases. In contrast, w/cratio is the least sensitive variable to the change in Pf. When computing the sensitivity measure with respect to the change in UEI(seeFig. 10 (c)), ARandSP/ CNTs ratio are the most sensitive variables to the changes in Pf. In addition, the negative or positive value of the sensitivity mea- sure for each variable indicates that Pfdecreases or increases, respec- tively, as the amount of that variable increases. For example, the positive values of sensitivity measure for ARand/or w/cratio indicates that the probability of failure increases as the value of ARand/orw/c ratio increases. On the other hand, Fig. 10 (a) shows that Pfdecreases (i.e., the chance of attaining at least 50% increase in the ﬂexural strength increases) as the value of t,UEI,s/cratio and SP/CNTs ratio increases. Fig. 10 (c) shows that as UEIincreases, the sensitivity mea- sure of ARgets closer to zero, indicating that higher ARis bene ﬁcial to increase the ﬂexural strength in case of a high UEI(i.e., proper CNT dispersion). 7.2. Importance analysis In the limit state function (Eq. (14)), some variables have larger impacts on the variance of the limit state function and some variables have smaller impacts, and therefore are less important. The vector of importance measures ( γ) is computed as follows [64]: γT¼αTJu/C3z/C3SD0 kαTJu/C3z/C3SD0kð17Þwhere αis row vector of the negative normalized gradient of the limit state function evaluated at the design point in standard normal space, z is vector of random variables ( z¼ðxp;ΘÞÞ,Ju/C3z/C3is Jacobian of the prob- ability transformation from the original space (i.e., z) to the standard normal space (i.e., u) with respect to the coordinates of the most likely failure point (i.e., design point, z/C3) and SD0is the standard deviation diagonal matrix of equivalent normal variables z0deﬁned by the lin- earized inverse transformation z0¼z/C3þJz/C3u/C3ðu/C0u/C3Þat the design point. The elements in the matrix SD0are the square root of the corre- sponding diagonal elements of the covariance matrix ∑0¼Jz/C3u/C3JT z/C3u/C3of the variables in the z0: The importance measure is used to rank the importance of different variables on the variability of the limit state function with respect to κ (seeFig. 11 (a)), AR(seeFig. 11 (b)) and UEI(seeFig. 11 (c)) when using mean values of other variables. In importance analysis, the higher absolute value of the importance measure vector ( γ) for each variable indicates the higher importance of that variable. Fig. 11 (a) shows that w/cratio is the most important variable at κ≤0.05 c‐wt%. The importance of w/cratio gradually decreases as κincreases and beyond κ= 0.3 c‐wt%, the w/cratio is the least impor- tant variable to the variation in the limit state function (the smallest absolute value of importance measure ( y‐axis) amongst other vari- ables). As κincreases, the importance of AR,tand UEIalso increases, while the importance of SP/CNTs ratio remains almost constant. For Fig. 9. Contour plots of probability of failure as a function of CNT concentration and aspect ratio using (a) mean of other variables (b) optimum values of othe r variables. Fig. 10. Sensitivity measures for ﬂexural strength as functions of (a) κ(b)AR(c)UEI.M. Ramezani et al. Composite Structures 253 (2020) 112748 10 example, tis the least important variable at κ< 0.05 c‐wt%. However, atκ= 0.5 c‐wt%,tis the second most important variable after AR. Fig. 11 (b) shows that as ARincreases, the importance of κ decreases. Below ARof around 1000, κis the most important variable to the change in Pf. However, beyond ARof 2500, κis the least impor- tant variable to the change of the ﬂexural strength. Fig. 11 (c) shows that when UEI<5 , ARandSP/CNTs ratio are the most important variables to the change in the ﬂexural strength, respec- tively. However, the importance of SP/CNTs ratio decreases by increasing in UEI and it becomes the least important variable at UEI > 10. In addition, as UEIincreases, the importance of κand w/c ratio increases. When UEI≥15,κis the most important variable to the variance of the limit state function. 8. Conclusions This paper proposes a bias ‐free probabilistic model based on the available deterministic Kelly ‐Tyson model to predict the ﬂexural strength of CNT ‐cement nanocomposites (pastes and mortars). The proposed model considers the in ﬂuences of CNT physical prop- erties, dispersion state, matrix composition and cement hydration using a Bayesian approach. The information obtained from the proba- bilistic model provides better understanding of the mechanisms ofCNTs improving the ﬂexural strength, leading to more ef ﬁcient designs. The following conclusions can be made: 1. The proposed probabilistic model is able to predict the ﬂexural strength within 15% error of the experimental test data. In addi- tion, the proposed model can be updated if new data becomes available in the future, leading to more accurate predictions. 2. Up to certain limits, the higher CNT aspect ratio ( AR) and concen- tration ( κ) result in superior improvement in the ﬂexural strength. Generally, ARranging from 400 to 800 and κranging between 0.08 and 0.18 c‐wt% (based on weight percent of cement powder) results in the highest improvement of the ﬂexural strength. 3. There are several options to achieve comparable improvements in theﬂexural strength with varied combinations of the in ﬂuential variables. For example, when using mean of other variables, utiliz- ingw/cand s/cratios of 0.4 and 3 or 0.35 and 2, respectively, results in similar improvements of the ﬂexural strength. 4. The AR,w/cratio and s/cratio are the most sensitive variables to the change in the ﬂexural strength with respect to κ. Also, AR and SP/CNTs ratio are the two most important variables with respect to UEI.A s κincreases, the importance of ARand age of specimen increases, while the importance of w/cand s/cratios decreases. In addition, as UEIincrease, the importance of SP/CNTs ratio to attain the superior ﬂexural strength decreases.CRediT authorship contribution statement Mahyar Ramezani: Conceptualization, Methodology, Data cura- tion, Software, Formal analysis, Writing ‐original draft. Young Hoon Kim: Conceptualization, Methodology, Writing ‐review & editing, Project administration, Resources, Funding acquisition, Supervision. Zhihui Sun: Writing ‐review & editing, Resources, Supervision. Declaration of Competing Interest The authors declare that they have no known competing ﬁnancial interests or personal relationships that could have appeared to in ﬂu- ence the work reported in this paper. Acknowledgments The authors wish to express gratitude and sincere appreciation of the partial ﬁnancial support by the Intramural Research Incentive Grant and Civil and Environmental Engineering Department at the Univer- sity of Louisville. Data availability statement The data that support the ﬁndings of this study are available on request from the corresponding author. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.compstruct.2020.112748 .