Computational characterization of the structural and mechanical properties of nanoporous titania

Abstract

Nanoporous titania is one of the most commonly used biomaterials with good biocompatibility and mechanical strength. Understanding to the influence of pore structures on their performances is crucial for the design and preparation of titania-based materials. Two kinds of structural models for nanoporous titania were constructed and used to investigate the effect of pore size and/or porosity on their mechanical properties by using molecular dynamic simulations with the Matsui–Akaogi potentials. The porous structures were relaxed and their elastic constants were computed and used to evaluated their bulk, shear and Young's moduli. Overlap effect in small pores, pore size and porosity have considerable influence on computed elastic moduli. Compared to bulk rutile TiO₂, reduced mechanical moduli were predicted. Simulations on uniaxial tensile tests revealed an anisotropic stress–strain relationship and a brittle-to-ductile transition for structures with large porosities. Fracture failure was predicted for all the studied porous structures. The maximum stress decreases with pore size and porosity, while the corresponding strain decreases with pore size, but increases with porosity.

---

1. Introduction

Porous materials can be classified into three types by their pore sizes: microporous (<2 nm), mesoporous (2–50 nm) and macroporous materials (>50 nm). Nanoporous materials are defined as those porous materials with pore diameters of 1–100 nm.¹ Nanoporous materials have shown great potential in catalysis, hydrogen storage, drug delivery and biological medicine, etc.^2–5 because of their high surface-to-volume ratios, moderate and tunable pore sizes, pore shapes, porosity and framework compositions. Being used as biomedical biomaterials for tissue engineering and drug delivery, some porous materials with good biocompatibility and bioreductivity have been characterized, such as porous titanium and porous calcium phosphates. In fact, implant three-dimensional substitute materials with porous structures are widely used in the treatment of bone trauma in clinics because of their good osteoreductivity that is essential for bone regeneration.^6–9 The pore structures provide ideal anchors for biomolecules/cell adhesion, proliferation, differentiation and formation of new organs. On the other hand, the porous materials by careful design may possess mechanic strengths that are sufficient to support damaged organs or tissues.^10,11 Typically, most biological implant materials for bone repair have a Young's modulus higher than those of hard tissues. Stress occlusion can impede bone remodeling and healing, resulting in increased bone porosity. One of the major problems associated with implants in orthopedic surgery is the mismatch in Young's modulus between natural bone and implants. One way to alleviate this problem is to reduce the Young's modulus of the materials by introducing holes,^12–14 thereby minimizing damage to tissues near the implant and ultimately extending device lifetime. In addition, bone/implant fixation is achieved by the mutual intersection between the bone and the porous implant matrix.¹⁵

Among various bone graft substitutes from natural and synthetic materials, titanium dioxide (TiO₂) is considered to be a highly suitable biocompatible material for bone anchoring implants due to its good biocompatibility, osteoconductivity and osteoinductivity.^16–19 In fact, TiO₂ exhibits little or no toxicity in both in vitro and in vivo experiments.²⁰ Extensive studies from their preparation, characterization, through applications have been conducted for various TiO₂-based materials. Fattakhova-Rohlfing et al.²¹ surveyed the synthesis of various porous TiO₂ nanostructures including membranes, spheres, fibers, etc. Sabetrasekh et al.²² reported that porous TiO₂ particles exhibit better activity in cell proliferation, but less toxicity than commercial bone grafts. Natix®, Straumann® BoneCeramic and Bio-Oss®. Haugen et al.²³ developed a high-interconnectivity and high-porosity TiO₂ that favors osseous integration, union and regeneration. Santos et al.²⁴ explored the mechanical properties of biomedical titanium dioxide films prepared by anodizing on commercial pure Ti (cp-Ti), including nanohardness, elastic modulus and brittleness. The modified films have relatively high nano-hardness and small elastic modulus in comparison with the cp-Ti. The reduction in modulus of elasticity was attributed to the porosity and inherent roughness of the oxide surface.

A number of computational studies have been conducted to explore the structures and properties of rutile titania and to help with the design of rutile-based materials. The computations focused mainly on the electronic, thermodynamic and mechanical properties of bulk or nano-sized TiO₂. For example, Mashreghi²⁵ calculated the thermal expansion coefficients of rutile titania nanoparticles under 300–1000 K by means of molecular dynamics (MD) simulations with the Buckingham potentials. Kim et al.²⁶ reproduced the lattice structures, thermal expansion coefficients and bulk modulus of rutile, brookite and anatase by using molecular dynamics simulations with the Morse interactions. Chung and Buessem²⁷ used Voigt–Reuss–Hill (VRH) approximation to calculate the elastic modulus of isotropic polycrystalline TiO₂ based on the anisotropy single crystal elastic constant. The results are in good agreement with the measured polycrystalline elastic modulus. Mahmood et al.²⁸ performed a first-principles calculation of rutile TiO₂ using the general gradient approximation (GGA) proposed by Perdew–Burke–Ernzerhof (PBE), and found that the single crystal elastic constant (C_ij), bulk modulus (B), Young's modulus (Y), shear modulus (G), linear bulk modulus (B_a) and the corresponding single crystal elastic constant (S_ij) alternately increase or decrease with pressure. However, most of those studies focused on non-porous crystalloid rutile TiO₂. In the research of mesoporous titania, Solveyra et al.²⁹ performed molecular dynamics simulations of water confined in mesoporous TiO₂-rutile pores at different water contents, obtaining water density and diffusion coefficient as a function of distance from the surface. Similarly, Velasco et al.³⁰ performed a detailed analysis of the hydrodynamics in the pore space of titanium dioxide by a combination of experimental nuclear magnetic resonance NMR data and MD simulation. Gautam et al.³¹ used MD simulation to study the structural and kinetic properties of propane restricted in mesoporous TiO₂.

Compared to tremendous studies on crystallized titanium dioxides, fewer researches have been conducted to their porous samples. Precise design and control over the pore structures is crucial for modulating the properties of porous TiO₂. Computational studies provide an effective way to establish the correlation between pore structures and performances. In this work we designed various computational models for porous titania with different pore sizes and porosities, and studied their structural and mechanical properties. Our purpose is to figure out the influence of pore structures on the properties of porous titania, which would be helpful for the design and preparation of titania-based materials.

2. Methodology

Titania has three phases in nature, rutile, anatase, and brookite.³² Rutile TiO₂ is thermodynamically stable in its bulk form with a large crystallite size at normal pressure and temperature up to its melting point 1830 K.³³ Rutile belongs to P4₂/mnm space group in which a titanium atom sits at the center of the unit cell and is surrounded by an octahedron formed by six oxygen atoms. The centered Ti atom has six coordinates, while each O atom is coordinated by three Ti atoms. Starting with the rutile TiO₂ structure, two series of computational models for porous titania were designed and constructed by excavating one cylinder along the z-axis direction. In Model I, all of the rutile TiO₂ porous structures were established at a fixed porosity of 8.1% and the cylinders are 1.3, 2.8, 3.4 and 5.1 nm in diameter, labeled as IA, IB, IC and ID, respectively. This model was used to evaluate the effect of pore size on the structures and properties of porous titania. In Model II, cylinders with a fixed diameter of 2.8 nm were excavated in the supercell. The porosity (from 8.1% to 22.6%) was controlled by changing the number of unit cells in the supercell. This model was used to evaluate the effect of porosity on the structures and properties of porous titania, labeled as IIA, IIB, IIC, IID and IIE, respectively. In both models, the Ti/O ratio is kept to 1 : 2 (ref. 29) and the wall thickness between adjacent pores are larger than 1 nm, as shown in Fig. 1. The three-dimensional structures of Model I and Model II were also given in Fig. S1.† The structural parameters of the constructed models are presented in Table S1.†

Fig. 1 Pore structures with different diameters and porosity in titania.

Arising from the extensive research interest, several kinds of force fields have been proposed for titanium dioxide.^34–46 The Matsui–Akaogi (MA) force fields⁴⁶ are commonly used to describe the TiO₂ systems. Other potentials used for the TiO₂ systems include modified-MA^42,43 and MS-Q.^44,45 Using molecular dynamics simulations with the MA force fields, Matsui and Akaogi⁴⁶ calculated the structure and physical properties of four kinds of TiO₂. Naicker et al.³⁶ and Fuertes et al.³⁸ simulated the surface energy and optical properties of TiO₂ nanoparticles using MD simulations with the MA potential. The MA potential was also used to simulate the surface structure of TiO₂ (ref. 35) and the deformation of anatase-type nanocrystals under uniaxial stretching and compression at room temperature.³⁷ In addition, the MS-Q potential^44,45 also produced accurate descriptions for the structures and properties of TiO₂ crystal.^39–41 The MA potential is defined as

where U_ij(r), C_ij and r_ij are the interaction energy, van der Waals (vdW) coefficient and distance between atoms i and j, respectively. ρ and q are ionic-pair dependent length parameter and atomic charge, respectively. ε₀ is the dielectric constant in vacuum. Alderman and Skinner et al.⁴³ further proposed a modified-MA potential aswhere σ is an interaction-dependent length parameter. The MS-Q potential model is defined aswhere A_ij, r₀, and B_ij are the parameters of Morse potential. r₀ is the equilibrium interatomic separation. The above mentioned three potentials were employed in this work. The parameters are listed in Table S2–S4 in ESI.†

Molecular dynamics simulations were carried out using the LAMMPS code.⁴⁷ The time integral of the Newton's equation of motion was performed using the velocity Verlet algorithm, and the time step was set to 1.0 fs. The Nosé–Hoover thermostat^48,49 was used to control the temperature. The periodic boundary conditions were applied in the three directions. All the structures were firstly relaxed with the conjugate gradient (CG) algorithm,⁵⁰ and then the constant-pressure constant-temperature (NPT) ensemble was employed for simulations at given temperatures. Particle–particle particle-mesh (PPPM) method was applied to minimize errors in long range terms coulombic. The cut-off of electrostatic interaction and van der Waals interaction used in all the simulation was set to 10 Å. Each simulation was carried out for 1.2 ns. The variations in potential energy with simulation time of all the structures are shown in Fig. S2–S4 in ESI.† All the systems reach their equilibrium at about 200–500 ps. The simulations continue and the data of last 400 ps were used for analysis.

The elastic modulus of rutile TiO₂ are computed with the elastic constants, which are obtained using the energy-strain theory:

andwhere V₀ and E(V₀,0) represent the volume and energy of the unstrained structure. μ, δ and C_ij denote the stress tensor, strain tensor and elastic constants respectively. Six independent elastic constants, C₁₁, C₁₂, C₂₃, C₃₃, C₄₄ and C₆₆, are defined for the tetragonal bulk and mecroporous rutile TiO₂ structures.⁵¹ The maximum strain was set to 1% in the evaluation of elastic constants. Based on the constants, bulk modulus (K), shear modulus (G) and Young's modulus (E) and Poisson ratio (η) were then computed with the Voigt–Reuss–Hill method.^52–54 The calculations of elastic constants were carried out for the equilibrated structures at 300 K. The details of K, G, E, η calculations are given in the ESI.† The volume thermal expansion coefficient (β) is defined aswhere V₀ and T₀ represent the initial reference volume and temperature. V_T represents the volume at temperature T. The β variations of titanium dioxide were explored at temperatures between 300 K and 1000 K with a temperature separation of 100 K.

3. Results and discussion

First, we examined the reliability of the potentials by examining the structures and properties of bulk rutile TiO₂. A 9 × 8 × 18 supercell with 7776 atoms was established. The simulated lattice parameters are given in Table 1, together with the measured and calculated results in previous studies.^55,56 The three potentials produce comparable results, and the MA and MS-Q results match well with the measurements.⁵⁶

Table 1 Lattice parameters of bulk rutile TiO₂

	a (Å)	b (Å)	c (Å)	α (°)	β (°)	γ (°)
a Measured in ref. 56.b First-principle calculations in ref. 55.
MA	40.515	36.013	54.332	90.00	90.00	90.00
Modified MA	43.722	38.864	54.844	90.00	90.00	90.00
MS-Q	40.155	35.694	52.772	90.00	90.00	90.00
Rao et al.^a	41.346	36.752	53.262	90.00	90.00	90.00
Camargo et al.^b	41.652	37.024	53.748	90.00	90.00	90.00

---

The computed six elastic constants of rutile TiO₂ are given in Table 2. Previous calculations^33,57 produced similar results with experiments.⁵⁸ Our simulations with the MA and the MS-Q potentials are in good agreement with other studies. The modified-MA results deviate considerably from the measurements, in particular to the shear moduli. The MA potential was thus selected to simulate the elastic properties of porous rutile TiO₂.

Table 2 Elastic constants (in GPa) and moduli (in GPa) of rutile TiO₂

	C11	C33	C12	C23	C44	C66	K	G
Calc.^33	292	471	192	147	114	236	233	113
Calc.^57	267	483	165	152	122	212	217	128
Expt.^58	268	484	175	147	124	190	210	113
MA	313	430	220	139	121	220	228	115
Modified-MA	283	392	309	185	143	155	255	1
MS-Q	331	535	186	158	129	211	241	132

---

Fig. 2 depicts the β values of bulk rutile TiO₂ calculated by using the three force fields as a function of temperature. In general, the β values increase with temperature, indicating that the volume of rutile TiO₂ positively expands when temperature increases. Touloukian et al.⁶⁰ and Saxena et al.⁵⁹ observed similar β values at high temperature of 700–1000 K, but distinct values at low temperature of 300–700 K. The three potentials predict a similar trend for the β values, which are however different in magnitude. The values by the modified-MA are the largest at the same temperature, followed by the MA and the MS-Q results. Our calculations with the MA and the MS-Q results are close to the measurements of Touloukian et al.⁶⁰ The MS-Q potential was thus selected to calculate the volume thermal expansion coefficients of the porous structures.

Fig. 2 Volume thermal expansion coefficient (β) of rutile TiO₂ at different temperatures. Experimental values are given for comparison.

Fig. S5 and S6,† which have similar shapes with Fig. 1, in ESI,† show the pore structures with different pore sizes and porosities in rutile TiO₂ after relaxation with the MA force field. Compared to Fig. 1, the atoms on the surface of the pores are rearranged to some extent. Both the surface oxygen and titanium atoms basically retain their interatomic connections with only small displacements around their original locations. The locations of inner atoms far away from the holes are almost unchanged. Table 3 lists the structural parameters of the pores in Model I and II. Compared to the unrelaxed structures (Table S1†), in both models the total volumes of the cubic boxes decrease after relaxation under the isobaric–isothermal ensemble. Accordingly, the pore sizes become smaller, but the changes are less than 4%. Atoms have a tendency to move toward the vacuum holes, resulting to the reduction of pore size and the distortion of pore edges. In both models, moreover, the surface areas of the pores and porosity also have a little decrease. Table 3 also presents the specific surface area (γ) of all the porous models. γ decreases with pore size in Model I, and increases with porosity in Model II. Fig. 3 shows the volume thermal expansion coefficient (β) of the porous rutile TiO₂ structures. In both models, β increases with temperature, which is in agreement with many other materials. At some temperatures, β increases with pore size and porosity, but, in general, their changes are rather small. Moreover, the β values are comparable with that of bulk rutile TiO₂.

Table 3 Structural properties of rutile TiO₂ pore structures after relaxation using MA force field

Model	Pore size (nm)	Surface area (Å^2)	Total volume (Å^3)	Porosity (%)	Specific surface area (×10^−2 Å^−1)
Model I
IA	1.2	1384	39 053	8.1	3.5
IB	2.7	2351	149 432	7.9	1.6
IC	3.3	3125	247 644	8.3	1.2
ID	5.0	4449	523 809	8.0	0.8
Model II
IIA	2.7	2351	149 432	7.9	1.6
IIB	2.7	2390	95 796	12.5	2.4
IIC	2.7	2377	76 179	15.7	3.1
IID	2.7	2396	59 139	20.4	4.1
IIE	2.7	2389	53 688	22.4	4.4

---

Fig. 3 Volume thermal expansion coefficients (β) of the porous rutile TiO₂ structures with different pore sizes (up) and porosities (bottom).

Six independent elastic constants, C₁₁, C₁₂, C₂₃, C₃₃, C₄₄ and C₆₆ were computed for the porous structures with the MA force fields, which are given in Table 4. The mechanical stability of these structures was examined using the following conditions:^61,62 C₁₁ > 0, C₃₃ > 0, C₄₄ > 0, C₆₆ > 0, C₁₁–C₁₂ > 0, C₁₁ + C₃₃ − 2C₂₃ > 0, and 2C₁₁ + C₃₃ +2C₁₂ + 4C₂₃ > 0. These conditions are satisfied by the computed results of all the structures. The results of bulk rutile TiO₂ are also given for comparison. For Model I, the elastic constants of R = 1.3 nm (IA) have rather large changes compared to the corresponding bulk ones. Because of the serious overlap effect of surface atoms in the small pore, IA exhibits different features from the other models. The overlap effect reduces in the models from IB, IC to ID. The computed elastic constants are comparable for IB and IC, and change to some extent for ID. Moreover, changes in C₁₁, C₁₂ and C₆₆ are more remarkable than in C₂₃, C₃₃ and C₄₄. For Model II, the elastic constants decrease with increasing porosity. Similarly, the decrease in C₁₁, C₁₂ and C₆₆ are more remarkable than in C₃₃, C₂₃ and C₄₄. Therefore, both pore size and porosity have influence on the magnitude of elastic constants, especially in C₁₁, C₁₂ and C₆₆.

Table 4 Elastic constants (GPa), elastic moduli (GPa) and Poisson's ratio of porous rutile TiO₂ structures

Model	C11	C33	C12	C23	C44	C66	K	G	E	η
Bulk	313	430	220	139	121	220	228	115	295	0.28
Model I
IA	204	382	132	87	96	222	153	87	220	0.26
IB	313	388	220	120	102	220	215	108	277	0.28
IC	310	384	213	117	102	213	211	107	275	0.28
ID	261	388	189	123	104	93	197	84	220	0.31
Model II
IIA	313	388	220	120	102	220	215	108	277	0.28
IIB	287	369	202	118	100	207	202	102	262	0.28
IIC	223	354	160	103	95	184	170	92	235	0.26
IID	167	334	106	78	87	140	129	78	196	0.24
IIE	131	324	79	68	83	116	107	72	176	0.22

---

Bulk modulus (K) and shear modulus (G) are used to measure the rigidity of a material resisting to the volume deformation and shape deformation under external forces. Young's modulus (E), a proportion of stress and strain, is customarily used as a measure of rigidity of a crystalline solid. Fig. 4 and 5 display the variation of these moduli with pore diameter and porosity, respectively. The moduli of bulk rutile TiO₂ are also given in the figures. At the same pore diameter or porosity, an order of E > K > G is always noted. Those three moduli are smaller than the corresponding bulk values. Except IA, the three moduli decrease with pore size and porosity, which again illustrates the overlap effect in small nanopores. The reduced moduli in porous rutile TiO₂ are in agreement with the observations for the rigidity reduction of other porous materials. For example, Li et al.⁶³ verified that the elastic moduli of porous nickel decrease as the pore size increases. The reduction was attributed to the increased load and moment on the pore walls. It is interesting to note that the decrease of E and G with porosity (Fig. 5) is approximately between linear and quadratic functions. Kováčik et al.^64,65 proposed two models of permeability, which are given as

E = E₀(1 − P/P_c)^f_E (7)

G = G₀(1 − P/P_c)^f_G (8)

where E and G are the effective Young's modulus and shear modulus of porous material, respectively with the porosity of P. E₀ and G₀ are the Young's moduli and shear moduli of bulk material. P_c is the porosity at which E or G becomes zero. f_E and f_G are the characteristic exponents of the porous samples. In general, the characteristic indices f_E and f_G are 1.10–1.70 and 1.00–1.40, complying with the computed trend shown in Fig. 5.

Fig. 4 Elastic moduli of porous rutile with different pore sizes. The corresponding values of bulk rutile are presented as dotted line.

Fig. 5 Elastic moduli of porous rutile with different porosities. The corresponding values of bulk rutile are presented as dotted line.

The ductility/brittleness variations have been observed in other porous materials. Meille et al.⁶⁶ reported that porous alumina, which has a fixed pore size, becomes brittler when its porosity increases. The variation was attributed to the thinner walls in the structures with larger porosities, which are more difficult to resist against shear loads. Poisson's ratio (η) represents the ratio of transverse strain to axial strain when a material is deformed along the direction of load. It is an elastic constant reflecting transverse deformation of a material. The material with a η indicates that the lateral deformation amount is larger than the longitudinal deformation before the plastic deformation occurs. On the contrary, the lateral deformation is smaller than the longitudinal deformation amount. The value of Poisson's ratio (η) is often used to estimate the ductility/brittleness with a critical point of 0.25 or lower for a brittle material and 0.33 or higher for a ductile material.^67,68 Bulk rutile TiO₂ is brittle, but little is known about the brittleness changes when the pores are introduced. The η values are between 0.26 and 0.31 for the structures in Model I and between 0.28 and 0.22 for the structures in Model II. However, one notes that the η values were computed by assuming the structures are isotropic, but the porous rutile TiO₂ are anisotropic. Therefore, the estimation with isotropic η values is not feasible for predicting the ductility/brittleness of the porous structures.

Uniaxial tensile measures the mechanical responses of a material subjected to a uniaxial load in an external environment. Molecular dynamics simulations have been conducted by many authors^69–72 to investigate the elastic limit, tensile strength, yield point, yield strength and other tensile properties of crystals, polymers and ceramics. The stress–strain behavior of porous rutile TiO₂ was simulated in this work by performing uniaxial tensile tests on the model structures. The simulations were carried out on the x, y and z directions of rutile TiO₂, which are perpendicular to the (001), (100) and (010) facets, respectively. In each stress–strain relationship, only one direction is stretched, while pressure in the other two vertical directions was set to zero. The strain rate of stretching was set to 1.0 × 10¹⁰ s⁻¹ and temperature was kept at 298 K. The stress–strain curves were obtained for all the structures in Model I and Model II, as shown in Fig. 6. The structures were stretched from 0 to 20% under an external load, and the corresponding stresses are σ_xx, σ_yy and σ_zz in the three directions, respectively. The fluctuations in the two directions perpendicular to the tensile direction are rather small and are not presented.

Fig. 6 Stress–strain curves of porous rutile with different pore sizes (top) and different porosities (bottom).

For the structures in Model I, their stress–strain curves follow the same behaviors: elastic at small strains, and then failure at great strains. However, the curves vary with pore size and are different from those of bulk rutile TiO₂. The porous structures have smaller elastic moduli than the bulk and their failure points are smaller than those of the bulk in the three directions. In the x and y directions, the structure with a large pore has a slightly smaller elastic modulus than the structure with a small pore. The maximum stresses σ_xx of Model IA–ID in the x direction are 17.4, 16.5, 15.4 and 12.7 GPa, the strains corresponding to the maximum stresses are 15.2, 13.6, 12.5 and 10.4%, respectively. In the y direction, the maximum stresses are 17.5, 16.2, 15.7 and 12.6 GPa and their corresponding strains are 16.0, 13.5, 12.9 and 10.3%, respectively. The maximum stresses σ_xx and σ_yy, as well as their corresponding strains decrease with pore size, indicating that the structures with larger pore sizes at the same porosity have lower tensile strengths. Therefore, the structure with a large pore tends to be fractured at a small load in the x and y directions. In the z direction, the maximum σ_zz values are larger those in the other two directions and have few changes with pore size. Their corresponding strains also have rather small changes, 12.7%, 12.5%, 12.7% and 12.6%, respectively. In the direction parallel to the cylindrical pore, therefore, the mechanical behaviors have smaller changes than the other two directions.

One notes that the curves for Model IA in the x and y directions fluctuate when the loads are close to the maximum stresses, which indicates a plastic deformation at a strong stress. Such plastic deformation was not noted in the other three structures in Mode I. IA has the smallest pore size and thinnest wall among the four structures when their porosity is fixed. The thin wall favors a plastic deformation for the rutile TiO₂. However, the plastic deformation changes into a brittle fracture when the strain increases further. Therefore, the stretching failure in all the structures in Model I is characterized by their brittle fracture, a sudden drop in stress at a given strain.

For the structures in Model II whose pore sizes are fixed and porosities are varying and larger than those of the structures in Model I, their stress–strain curves reflect elastic responses at small strains. However, the curves deviate from the elastic responses at larger strains. The stresses fluctuate with strain, which represents a plastic deformation at this stage. The deformation is more remarkable for the structures with large porosities, resulting from their thinner wall thickness. Brittle fractures occur when the stains increase further. Over the studied stain range, elastic and inelastic, the moduli of the porous structures are smaller than those of the bulk. The stress–strain relationship varies with porosity. In the x direction, the maximum stress reaches 16.73, 14.87, 13.9, 12.09 and 12.02 GPa at a strain of 13.68%, 14.2%, 14.78%, 14.58% and 15.54%, respectively. The maximum stress decreases with porosity. The structure with a larger porosity has a lower maximum stress, which, however, corresponds to a larger strain. This is in contrast to the variations found in Model I. The structures in Model II have larger porosities than those on Model I. When the external load is larger than 10%, the stress varies in an inelastic pattern with the strain. The deviation from elastic response is more remarkable for the structures with larger porosities. The five structures are in the order of IIA > IIB > IIC > IID > IIE. The brittle-to-ductile transition results from the rearrangement of atoms around the pores, leading to a small stress/stain over a certain range of strains. Therefore, IIE has the smallest maximum stress but the largest failure strain. Similar variations were noted in the stress–strain curves in the y direction, the maximum stress and the failure strain vary with porosity, but their maximum stresses are lower than those in the x direction. The brittle-to-ductile transition was also noted. In the z direction, the variations with porosity become less remarkable. The five structures have almost the same stress–strain relationship at small strains. Their maximum stresses are higher than those in the other two directions, and are slightly higher for the structures with larger porosities. The maximum stresses correspond to almost the same strain values. The brittle-to-ductile transition is not clear in this direction, but brittle fractures are noted for all the five structures.

4. Conclusion

Two structural models were constructed to mimic the structures of porous rutile TiO₂. In the first model, the four structures have different pore sizes of 1.3, 2.8, 3.4 and 5.1 nm in diameter, respectively, but have the same porosity. In the second model, the five structures have different porosities from 8.1%, 12.7%, 15.9%, 20.5% and 22.6%, respectively, but have the same pore size. The influence of pore size and porosity on the structures and mechanical properties of porous rutile TiO₂ were then investigated using MD simulations. The MA potential was selected to predict the structures and elastic properties, and the MS-Q potential was used for computing the volume expansion coefficients.

In the relaxed porous structures, the surface atoms basically retain their interatomic connections with only small displacements around their original locations. The locations of inner atoms far away from the holes are almost unchanged. Their volume thermal expansion coefficients are comparable with that of bulk rutile TiO₂, and increase with pore size and porosity at the same temperature, but the increments are rather small. The computed elastic constants vary with pore size and porosity, especially for C₁₁, C₁₂ and C₆₆. The structure with the smallest pore size of 1.3 nm exhibit exceptional variations in its elastic constants because of the overlap effect in its small pore. The elastic modulus (E, K and G) were evaluated based on the computed elastic constants. The three moduli are smaller than the corresponding bulk values, decreasing with pore size and porosity.

The uniaxial tensile tests were performed for all the structures in the two porous models with MD simulations. The stress–strain curves were obtained in three directions. The stress varies in an anisotropic way. In the z direction parallel to the cylindrical pore the stress behaves close to that of the bulk. In the other two directions, the maximum stress decreases with pore size and porosity. At small strains, elastic response was noted for all the structures. When the strain increases, inelastic response was noted for the structures with large porosities. The deviation from elastic response increases with porosity, and brittle-to-ductile transition was noted. When the strain increases further, brittle fracture where a sudden drop in stress occurs was noted for all the structures. The brittle failure occurs at smaller strains for the structures with larger pores in Model I, but at larger strains for the structures with large porosities, resulting from the brittle-to-ductile transitions. Our simulations show that the mechanical behaviors of porous structures depend closely with pore structures, including their size, density, etc. Pore shape and their connectivity may also have great influence on their mechanical properties, which will be investigated in our future work.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Key R&D Program of China (2017YFB0702600, 2017YFB0702601), the National Natural Science Foundation of China (21773159), Science Challenge Project (TZ2016001) and financial support from National High Technology Research and Development Program of China (No. 2015AA034202). Part of the calculations were carried out at High-Performance Computing Center of Sichuan University.

References

1. G. Q. Lu and X. S. Zhao, Nanoporous Materials: Science and Engineering, Imperial College Press, London, 2003.
2. C. Mertens, D. Wiens, H. G. Steveling, A. Sander and K. Freier, Clin. Implant Dent. Relat. Res., 2014, 16, 365–373.
3. S. Park and Sujandi, Curr. Appl. Phys., 2008, 8, 664–668.
4. P. Pei, Z. F. Tian and Y. F. Zhu, Microporous Mesoporous Mater., 2018, 272, 24–30.
5. M. E. Davis, Nature, 2002, 417, 813–821.
6. D. Green, D. Walsh, S. Mann and R. O. C. Oreffo, Bone, 2002, 30, 810–815.
7. L. H. Nguyen, N. Annabi, M. Nikkhah, H. Bae, L. Binan, S. Park, Y. Kang, Y. Yang and A. Khademhosseini, Tissue Eng., Part B, 2012, 18, 363–382.
8. C. Szpalski, M. Wetterau, J. Barr and S. M. Warren, Tissue Eng., Part B, 2012, 18, 246–257.
9. K. Gulati, H. J. Moon, T. Li, P. T. Sudheesh Kumar and S. Ivanovski, Mater. Sci. Eng., C, 2018, 91, 624–630.
10. K. Rezwan, Q. Z. Chen, J. J. Blaker and A. R. Boccaccini, Biomaterials, 2006, 27, 3413–3431.
11. S. J. Hollister, Nat. Mater., 2005, 4, 518–524.
12. M. Takemoto, S. Fujibayashi, M. Neo, J. Suzuki, T. Kokubo and T. Nakamura, Biomaterials, 2005, 26, 6014–6023.
13. M. Niinomi, M. Nakai and J. Hieda, Acta Biomater., 2012, 8, 3888–3903.
14. C. E. Wen, Y. Yamada, K. Shimojima, Y. Chino, H. Hosokawa and M. Mabuchi, J. Mater. Res., 2002, 17, 2633–2639.
15. I. K. Triantafillopoulos and N. A. Papaioannou, The Effect of Pharmacological Agents on the Bone-Implant Interface, Springer London, London, 2014, pp. 221–237.
16. J. S. Carson and M. P. Bostrom, Injury, 2007, 38(suppl. 1), 533–537.
17. F. D. Beaman, L. W. Bancroft, J. J. Peterson and M. J. Kransdorf, Radiol. Clin. North Am., 2006, 44, 451–461.
18. S. T. Kao and D. D. Scott, Oral Maxillofac. Surg. Clin., 2007, 4, 513–521.
19. Y. Khan, M. J. Yaszemski, A. G. Mikos and C. T. Laurencin, J. Bone Jt. Surg., Am. Vol., 2008,(suppl. 1), 36–42.
20. G. Gertler, G. Fleminger and H. Rapaport, Langmuir, 2010, 26, 6457–6463.
21. D. Fattakhova-Rohlfing, A. Zaleska and T. Bein, Chem. Rev., 2014, 114, 9487–9558.
22. R. Sabetrasekh, H. Tiainen, S. P. Lyngstadaas, J. Reseland and H. Haugen, J. Biomater. Appl., 2010, 25, 559–580.
23. H. J. Haugen, M. Monjo, M. Rubert, A. Verket, S. P. Lyngstadaas, J. E. Ellingsen, H. J. Ronold and J. C. Wohlfahrt, Acta Biomater., 2013, 9, 5390–5399.
24. E. Santos, N. K. Kuromoto and G. A. Soares, Mater. Chem. Phys., 2007, 102, 92–97.
25. A. Mashreghi, Comput. Mater. Sci., 2012, 62, 60–64.
26. D. Kim, N. Enomoto, Z. Nakagawa and K. Kawamura, J. Am. Ceram. Soc., 1996, 79, 1095–1099.
27. D. H. Chung and W. R. Buessem, J. Appl. Phys., 1968, 39, 2777–2782.
28. T. Mahmood, C. Cao, W. S. Khan, Z. Usman, F. K. Butt and S. Hussain, Phys. B Condens. Matter, 2012, 407, 958–965.
29. E. González Solveyra, E. de la Llave, V. Molinero, G. J. A. A. Soler-Illia and D. A. Scherlis, J. Phys. Chem. C, 2013, 117, 3330–3342.
30. M. I. Velasco, M. B. Franzoni, E. A. Franceschini, E. Gonzalez Solveyra, D. Scherlis, R. H. Acosta and G. J. A. A. Soler-Illia, J. Phys. Chem. C, 2017, 121, 7533–7541.
31. S. Gautam and D. Cole, Chem. Phys., 2015, 458, 68–76.
32. D. Dambournet, I. Belharouak and K. Amine, Chem. Mater., 2010, 22, 1173–1179.
33. J. Zhu, J. X. Yu, Y. J. Wang, X. R. Chen and F. Q. Jing, Chin. Phys. B, 2008, 17, 2216–2221.
34. S. Ogata, H. Iyetomi, K. Tsuruta, F. Shimojo, R. K. Kalia, A. Nakano and P. Vashishta, J. Appl. Phys., 1999, 86, 3036–3041.
35. R. S. Kavathekar, P. Dev, N. J. English and J. M. D. MacElroy, Mol. Phys., 2011, 109, 1649–1656.
36. P. K. Naicker, P. T. Cummings, H. Z. Zhang and J. F. Banfield, J. Phys. Chem. B, 2005, 109, 15243–15249.
37. X. Zhang, H. Gao and X. Li, Extreme Mech. Lett., 2018, 22, 131–137.
38. V. C. Fuertes, C. F. Negre, M. B. Oviedo, F. P. Bonafe, F. Y. Oliva and C. G. Sanchez, J. Phys.: Condens. Matter, 2013, 25, 115304.
39. B. S. Thomas, N. A. Marks and B. D. Begg, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 144122.
40. A. Hallil, R. Tétot, F. Berthier, I. Braems and J. Creuze, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 165406.
41. X. J. Han, L. Bergqvist, P. H. Dederichs, H. Müller-Krumbhaar, J. K. Christie, S. Scandolo and P. Tangney, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 134108.
42. L. B. Skinner, C. J. Benmore, J. K. Weber, S. Tumber, L. Lazareva, J. Neuefeind, L. Santodonato, J. Du and J. B. Parise, J. Phys. Chem. B, 2012, 116, 13439–13447.
43. O. L. G. Alderman, L. B. Skinner, C. J. Benmore, A. Tamalonis and J. K. R. Weber, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 094204.
44. V. Swamy and J. D. Gale, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 5406–5412.
45. A. K. Rappé and W. A. Goddard, J. Phys. Chem., 1991, 95, 3358–3363.
46. M. Matsui and M. Akaogi, Mol. Simul., 1991, 6, 239–244.
47. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
48. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
49. W. G. Hoover, Phys. Rev. A, 1985, 31, 1695–1697.
50. R. Fletcher and C. M. Reeves, Comput. J., 1964, 7, 149–154.
51. W. J. Alton and A. J. Barlow, J. Appl. Phys., 1967, 38, 3817–3820.
52. W. Voigt, Lehrbuch der kristallphysik, Taubner, Leipzig, 1928.
53. R. Hill, Proc. Phys. Soc., 1952, 65, 349–354.
54. A. Reuss, Z. Angew. Math. Mech., 1929, 9, 49–81.
55. A. C. Camargo, J. A. Igualada, A. Beltrán, R. Llusar, E. Longo and J. Andrés, Chem. Phys., 1996, 212, 381–391.
56. K. V. K. Rao, S. V. N. Naidu and L. Iyengar, J. Am. Ceram. Soc., 1970, 53, 124–126.
57. X. G. Ma, P. Liang, L. Miao, S. W. Bie, C. K. Zhang, L. Xu and J. J. Jiang, Phys. Status Solidi B, 2009, 246, 2132–2139.
58. D. G. Isaak, J. D. Carnes, O. L. Anderson and H. Cynn, Phys. Chem. Miner., 1998, 26, 31–43.
59. S. K. Saxena, N. Chatterjee, Y. Fei and G. Shen, Thermodynamic Data on Oxides and Silicates, Springer, 1993.
60. Y. S. Touloukian, R. K. Kirby and R. E. Taylor, Thermal Expansion: Nonmetallic Solids, in Thermophysical Properties of Matter, Plenum, New York, 1970.
61. M. Grimsditch, A. Polian, V. Brazhkin and D. Balitskii, J. Appl. Phys., 1998, 83, 3018–3020.
62. G. V. Sin'ko and N. A. Smirnov, J. Phys.: Condens. Matter, 2002, 14, 6989–7005.
63. Z. Li, L. Yang, Y. Li, Y. Yang, C. Zhou, Y. Ding, J. Zhao and Y. Li, Mater. Des., 2013, 45, 52–55.
64. J. Kováčik, J. Mater. Sci. Lett., 1999, 18, 1007–1010.
65. J. Kovácik, J. Mater. Sci. Lett., 2001, 20, 1953–1955.
66. S. Meille, M. Lombardi, J. Chevalier and L. Montanaro, J. Eur. Ceram. Soc., 2012, 32, 3959–3967.
67. V. V. Bannikov, I. R. Shein and A. L. Ivanovskii, Phys. Status Solidi RRL, 2007, 1, 89–91.
68. H. Z. Fu, D. H. Li, F. Peng, T. Gao and X. L. Cheng, Comput. Mater. Sci., 2008, 44, 774–778.
69. M. S. Lavine, N. Waheed and G. C. Rutledge, Polymer, 2003, 44, 1771–1779.
70. X. Ma, L. Zhang, H. Wu, X. Zhang and M. Yang, J. Mol. Model., 2017, 23, 156.
71. F. Yuan and L. Huang, J. Non-Cryst. Solids, 2012, 358, 3481–3487.
72. D. Hossain, M. A. Tschopp, D. K. Ward, J. L. Bouvard, P. Wang and M. F. Horstemeyer, Polymer, 2010, 51, 6071–6083.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9ra02298h

---