\documentclass[11pt,reqno]{amsart} % Use if you need to include figures \usepackage{graphicx} \usepackage{mathtime} % Comment this out if you don't have this package \textheight=51pc \textwidth=38pc \title[Short title of up to 45 letters]{Full title of paper} \author{Name of first author} \address{Affiliation of first author\\ Full postal address\\ Anytown, XX 12345\\ Country} \email{Email of first author} \urladdr{(Optional) home page} \author{Name of second author} \address{Affiliation of second author\\ Full postal address\\ Anytown, XX 12345\\ Country} \email{Email of second author} \keywords{Optional here; must be supplied on submission form} \thanks{This is an acknowledgement to a funding agency. Acknowledgments to individuals are better in a section by themselves (see end of paper); this way an entry appears in the online index.} % Author definitions come here \def\thesubsection{\thesection\Alph{subsection}} \numberwithin{equation}{section} \def\theequation{\thesection\hbox{-}\arabic{equation}} \newtheorem{example}{Example}[section] % examples numbered within sections \newtheorem{theorem}[example]{Theorem} % theorems numbered with examples \begin{document} % must come before the abstract \begin{abstract} A brief abstract of about 150 words or less must be included. The abstract should be self-contained and not make any reference to the bibliography. \end{abstract} \maketitle % must come after the abstract \section{Introduction}\label{intro} Also acceptable is an untitled introduction (omit section header) or an introduction without section number (use \verb+\section*{Introduction}+). You should use symbolic labels for your cross references. For example, at the beginning of this section we wrote \verb+\label{intro}+. Later, when we want to refer to this section, we type $$\verb+Section~\ref{intro}+,$$ which appears as Section~\ref{intro}. The same thing can be done with examples, theorems and equations. \section{Some examples}\label{examples} In this section we look at the preferred numbering of subsections, theorems, examples and so on. \begin{example}\rm \label{ex1} This is our first example. Inside this example we have the equation \begin{equation} \label{einstein} e=mc^2. \end{equation} \end{example} This equation is numbered because it is important and will be referred to later; this is also done with \verb+\label+ and \verb+\ref+. However, an equation that won't be referred to later doesn't need a number, even if it is displayed (within double dollar signs). Too many numbered equations clutter the page for no good reason. \begin{theorem}[Pythagoras] \label{pyth} In a right triangle, the square on the hypotenuse equals the sum of the squares on the perpendicular sides. \end{theorem} \begin{example}\rm This is our second example; it comes after Example~\ref{ex1} and Theorem~\ref{pyth}, so it gets a higher number than both. We prefer to have examples, theorems and other important statements all in a single sequence; this helps the reader by providing plenty of guideposts. (If Lemma 1 comes after Theorem 2 and so on, the only way to find results is to look through the whole paper.) \end{example} \subsection{About subsections} We also like to identify subsections with a letter to distinguish their numbers from those of examples and theorems. \section{Lists and enumerations} Lists can be either numbered (enumerations) or unnumbered. The default numbering for top-level, second-level and third-levels lists is this: \begin{enumerate} \item an entry \item another entry \begin{enumerate} \item a subentry \item another subentry \begin{enumerate} \item a subsubentry \item another subsubentry \end{enumerate} \end{enumerate} \end{enumerate} It there is good reason to change the numbering method, you can do so; consult any \LaTeX\ book. % ftp://ftp.ams.org/pub/tex/doc/amsmath/short-math-guide.pdf \section{Figures} \end{document} Pure titanium and titanium-based alloys exhibit excellent mechanical and biological properties, thus make titanium-based foams potential materials for load-bearing sandwich cores and orthopaedic implants \cite{Banhart, Spoerke, Wen2002a, Wen2002b}. In particular, as a potential implant material, stiffness of the porous material drops with the square of relative density to be comparable to bone stiffness and the open porosity allows complete bone ingrowth \cite{Gibson, Spoerke}. These properties make porous titanium a promising solution to the inherited problems of the monolithic metallic implants such as the ``stress shielding'' effect \cite{Chang, Dunand, Li2004, Spoerke, Wen2002a, Wen2002b}. However, the porous microstructure of the foam leads to stress and strain concentrations near pores under load-bearing conditions, resulting in reduction of strength and ductility. Finite element (FE) simulation on a microstructural level is therefore needed to understand and predict the macroscopic and microscopic responses to better target and optimize the \end{enumerate} \end{document} Pure titanium and titanium-based alloys exhibit excellent mechanical and biological properties, thus make titanium-based foams potential materials for load-bearing sandwich cores and orthopaedic implants \cite{Banhart, Spoerke, Wen2002a, Wen2002b}. In particular, as a potential implant material, stiffness of the porous material drops with the square of relative density to be comparable to bone stiffness and the open porosity allows complete bone ingrowth \cite{Gibson, Spoerke}. These properties make porous titanium a promising solution to the inherited problems of the monolithic metallic implants such as the ``stress shielding'' effect \cite{Chang, Dunand, Li2004, Spoerke, Wen2002a, Wen2002b}. However, the porous microstructure of the foam leads to stress and strain concentrations near pores under load-bearing conditions, resulting in reduction of strength and ductility. Finite element (FE) simulation on a microstructural level is therefore needed to understand and predict the macroscopic and microscopic responses to better target and optimize the application of these porous materials. The titanium foam considered here was processed by a solid-state foaming technique in which individual high-pressure argon pores expand at elevated temperature and coalesce to form large pores \cite{Dunand, Murray}. As porosity is less than 25{\%}, pores are mostly rounded, generally equiaxed, and un-merged, as illustrated in Figure~\ref{fig1}(a). As the material is foamed to higher porosity, small pores coexist with large pores which have a complex, tortuous shape as shown in Figure~\ref{fig1}(b) \cite{Murray}. The microstructures are locally heterogeneous, although sufficiently large experimental samples behave homogeneously. Regarding such a heterogeneous material, a question of fundamental importance is raised, i.e., whether a FE model simulating a fragment of the heterogeneous microstructure is large enough to be a representative volume element (RVE) to describe the responses of the titanium foam. The determination of the size of RVE relies on basic definitions. A few definitions of the minimum RVE have been developed for the study of heterogeneous materials. One definition was proposed by \cite{Hill}, according to which the RVE was described as ``\textit{a sample that (a) is structurally entirely typical of the whole mixture on average, and (b) contains a sufficient number of inclusions for the apparent overall moduli to be effectively independent of the surface values of traction and displacement, so long as these values are macroscopically uniform}'' \cite{Hill}. In this definition, (a) requires the RVE to statistically include all possible microstructural configurations and (b) demands the effective properties obtained from the RVE to be independent of the uniform displacement and traction boundary conditions (BCs). Regarding the (b) part in the definition, \cite{Huet} analytically proved that the effective elastic modulus obtained from a RVE is bounded by the average apparent responses of finite size domains under uniform displacement boundary condition (UDBC) and uniform traction boundary condition (UTBC) \cite{Huet}. This conclusion was extended to a nonlinear elastic heterogeneous material by \cite{Hazanov} and to elastoplastic materials under proportional loading by \cite{Jiang2001a, Hazanov, Jiang2001b}. Following this line, homogenization convergence was shown by many researchers that as the domain size increases to the minimum RVE, the two bounds converge to the effective property \cite{Hollister, Huet, Jiang2001a, Ostoja-Starzewski, Pecullan, Terada}. However, the convergence of the two bounds was found to be extremely slow for heterogeneous materials with soft inclusions in a hard matrix, and the higher the contrast between the moduli of matrix and inclusions, the slower the convergence \cite{Bouyge, Jiang2001a, Ostoja-Starzewski, Pecullan}. Given that zero modulus pores provide extreme moduli contrast to the titanium matrix for the titanium foam, the minimum RVE according to Hill's definition becomes prohibitively large to be used for the FE analyses, especially for three dimensional (3D) FE analyses due to the limitation of computational power. Another more pragmatic definition of RVE was introduced by \cite{Drugan}: ``\textit{the smallest material volume element of the composite for which the usual spatially constant `overall modulus' macroscopic constitutive representation is a sufficient accurate model to represent mean constitutive response}.'' \cite{Drugan}. Base on this RVE concept, they analytically proved that it is possible to obtain relatively accurate estimations of elastic constants with small RVEs for a microstructure with non-overlapping spherical inclusions. In particular, the effective moduli obtained by a finite size domain at a length of only two sphere diameters can be close to those by an infinite length domain within a few percent of error even in case of void inclusions. This conclusion is contradictory to the results obtained according to Hill's definition. The BC applied on RVE is the key factor responsible for the difference. While the RVE should be independent of BCs according to Hill's definition, the right BC is needed for the RVE defined by Drugan and Willis. As pointed out by \cite{Jiang2001a}, the derivation of Drugan and Willis for the minimum RVE implied the RVE should be associated with periodic BC and the microstructure is actually periodic \cite{Jiang2001a}. Jiang et al.\ verified the RVE theory of Drugan and Willis for unidirectional fiber--matrix composites with the modulus contrast of inclusion to matrix in a range of 0.001 to 1000 \cite{Jiang2001a}. Gusev numerically verified that the average of elastic constants of many small RVEs is close to the result for large RVEs for 3D microstructures with randomly distributed non-overlapping spherical inclusions under periodic BC \cite{Gusev}. Their results also indicate that smaller analyzed windows lead to more scattered results. Following the concept of the RVE defined by Drugan and Willis but accounting for the fluctuations of the apparent properties predicted by finite size domains, \cite{Kanit} recently proposed a more quantitative definition \cite{Kanit}. Based on their definition, the concept of one single minimum RVE should be abandoned and a sufficient number of small domain RVEs should be used to obtain the average linear properties. While UTBC and UDBC on small domains result in large, oppositely biased errors for the effective property, they found a periodic BC gives smaller error for the same window size. However, it should be noted that periodic BC requires the continuity of the inclusions on opposite boundaries to ensure the periodicity of the microstructure. Because such unnatural periodicity is seldom observed in real heterogeneous materials, periodic BCs are not appropriate for FE models developed by cutting out fragments of actual microstructures or simulated microstructures based on actual microstructures. Other BCs might be considered as alternatives which maintain the philosophy of the RVE concept of Kanit et al.\ emphasizing the statistic average sense of the responses. \cite{HazanovandHuet} proved analytically that the elasticity tensor predicted by a model smaller than Hill's RVE definition submitted to mixed BCs falls between the predictions associated with UTBC and UDBC \cite{HazanovandHuet}. They concluded that it is possible to get relatively accurate results by using small domains under mixed BCs. Jiang et al.\ numerically verified this conclusion on unidirectional fiber-matrix composites \cite{Jiang2001a}. However, what kind of mixed BCs to be associated with RVEs to obtain accurate estimates for porous microstructure is still an open question. No satisfactory solution for practical RVEs which can be solved with reasonable computational power, to the best of our knowledge, exists yet for porous microstructure, due to the extreme property contrast between inclusion and matrix. The objective of this paper is to find reliable and practical RVEs for porous titanium combining the homogenization convergence concept and statistical averages from solving the boundary value problems. As discussed above, solutions of finite size domains under mixed BCs approach effective properties much faster than UDBC and UTBC. However, since results from finite size domains under a single mixed BC approach the effective properties from one side, it is difficult to determine the point where the finite size domains are large enough. In addition, such a one side bias can only be eliminated by increasing the size of analyzed domains. Therefore, two mixed BCs simulating a uniaxial loading condition providing opposite bias are designed and imposed on the porous models in this paper in order to determine the convergence of RVEs and obtain relatively accurate results with small RVEs by averaging the results associated with the two BCs. By this approach, we trade-off the large RVEs with two mixed BCs and more small domains. In other words, the prohibitively large RVE according to the Hill's definition are divided into to reasonable small ones associated with special BCs. While previous definitions of RVE focused only on macroscopic linear properties, our current study also aims at verifying the convergence of both the macroscopic and the microscopic elastoplastic responses. The results of this paper elucidate a method for numerical prediction of the global and local response with small size models for porous materials. It should be noted that a porosity of $40\sim 50${\%} is generally considered to be an ideal range for orthopaedic applications of porous titanium providing reduced stiffness (for reduced stress-shielding) as well as sites for bone ingrowth while maintaining mechanical durability simultaneously. However, the microstructure at such high porosity is extremely complex (see, e.g., Figure~\ref{fig1}(b)), complicating the current RVE study. Since we focus on the methodology exploration in the current study, low porosity titanium foam is used for the study and the results and implications for a 3D model with a higher porosity of 48{\%} are discussed. The RVE for higher porosity titanium foams will be considered in our future work. In next sections, 2D and 3D FE models with different sizes are created based on the simulated microstructure of the experimental material using the methodology presented in our previous work \cite{Shenetal}. BCs for the FE simulations are described. Then the uniaxial stress-strain response is simulated based on the FE models associated with the two mixed BCs. Both macroscopic and microscopic responses are demonstrated for the convergence study. \section{Finite element modeling}\label{sec:2} \subsection{Simulated microstructure}\label{subsec:2.1} A simulated microstructure was generated based on the actual microstructure of the titanium foam at 14.7{\%} porosity using the methodology presented in our previous work \cite{Shenetal}. 3D pore size and location distribution information was first derived from 2D sections of a sample at 14.7{\%} porosity as shown in Figure~\ref{fig1}(a). The location distribution of pores is random. The pore size distribution follows Weibull distribution of which the density function is denoted by $f(x)$ as follows \cite{Tobias}: \begin{equation} \label{eq1} f(x) = \frac{m}{x}\Bigl(\frac{x}{c}\Bigr)^me^{ - (x / c)^m}, \end{equation} where the parameters, $m$ and $c$, are shape and scale parameter, respectively. $m = 2.47$ and $c = 64.68$ were obtained for the pore sections observed in 2D sectioning planes providing average diameter of $63.3\,\mu$m; $m = 2.29$ and $c = 60.52$ were derived by stereology study providing 3D average pore diameter of $53.6\,\mu$m \cite{Shenetal}. A 3D simulated microstructure was developed based on this data which retains the essential geometry features of the random microstructure. The microstructure of the foam at higher porosity was achieved by a geometric simulation of pore growth and movement during the foaming process. In the simulated 3D multi-pore microstructure (3D-MP-70271) at porosity 14.7{\%}, 70,271 spherical pores were randomly located in a $(4000\,\mu\textrm{m})^{3}$ cube. No pore impingement was permitted. Sections parallel to the $x_{1} x_{2}$ plane passed through the 3D-MP-70271 model were cut every $10\,\mu$m along $x_{3}$ axis to simulate the serial sectioning process which is a reliable technique for constructing microstructures \cite{Li1998}. The distribution of pore area fractions in window sizes of $1000\,\mu$m, $2000\,\mu$m, $3000\,\mu$m, and $4000\,\mu$m is shown in Figure~\ref{fig2}. It can be seen that as the size of the window increases, the range and variation of the area fraction decreases and the area fractions converge to the overall porosity of the microstructure. The simulated microstructure is comprised of two distinct phases, pores and titanium matrix. To avoid computational difficulties, it was assumed that pores are linear elastic with a very low modulus of $10^{ - 7}$\,GPa and Poisson's ratio of 0.3. The titanium matrix having elastic modulus of 110 GPa, Poisson's ratio of 0.33, and yield strength of 275\,MPa is representative of CP Ti-40 \cite{ASM}. The matrix yield surface follows the Von~Mises yield criterion with isotropic hardening. All material properties and simulations are for room temperature. \subsection{2D and 3D finite element models}\label{subsec:2.2} 2D and 3D finite element models were created by cutting out the fragments of the simulated microstructure. To study the convergence of responses of various size windows, four groups of 2D sections with sizes (side length) from $1000\,\mu$m to $4000\,\mu$m having area fraction matching the overall porosity of 14.7{\%} were selected for constructing 2D models. The four groups of 2D finite element models with side lengths of $1000\,\mu$m, $2000\,\mu$m, $3000\,\mu$m, and $4000\,\mu$m have ratios of model length to the average pore size $(L/d)$ of 15.8, 31.6, 47.4, and 63.2, respectively. Each model group contains eight models and is named based on their size. For example, 2D-1000 refers to the model group with side length of $1000\,\mu$m. Examples of these 2D model groups are shown in Figure~\ref{fig3}. It can be seen that some pores having a centroid near an edge were truncated. On average, there are 59 pores and 10,930 elements in the 2D-1000 model; 235 pores and 39,010 elements in the 2D-2000 model; 542 pores and 88,385 elements in the 2D-3000 model; 943 pores and 159,768 elements in the 2D-4000 model. Two 3D model groups with side lengths of $170\,\mu$m and $340\,\mu$m having porosity of 14.7{\%} were selected from different locations of the 3D-MP-70271 microstructure. Each group contains four 3D models. The ratios of model size to the average pore size $(L/d)$ are 3.17 and 6.34, respectively. Models were named 3D-170 and 3D-340 after the model size. Examples of the 3D-170 and 3D-340 models are shown in Figure~\ref{fig4}. Some pores having a centroid near a face were truncated. On average, the 3D-170 model contain 12 pores and 15,439 elements; the 3D-340 model 64 pores and 125,536 elements. One 3D model with higher porosity of 48{\%} in a $340\,\mu$m cube was created by cutting a block of further foamed microstructure. As shown in Figure~\ref{fig5}, pores are connected to form large pores which have a complex, tortuous shape. The model contains 130,552 elements. All 2D models were meshed with eight-node biquadratic plane strain elements and all 3D models were meshed with ten-node modified tetrahedral elements with hourglass control to prevent volumetric lock during plastic deformation \cite{Hibbitt}. All finite element analyses were performed using ABAQUS software. Mesh convergence was verified based on overall and local stress values. The overall strain and stress are calculated by using volume averages, \begin{align} \label{eq2} \bar {\varepsilon }_{ij} &= \frac{1}{V}\int_V {\varepsilon _{ij} } dV = \frac{1}{\sum_{m = 1}^N {V^{(m)}} }\sum_{m = 1}^N {\varepsilon _{ij}^{( m )} V^{( m )}} \\ \label{eq3} \bar {\sigma }_{ij} &= \frac{1}{V}\int_V {\sigma _{ij} } dV = \frac{1}{\sum_{m = 1}^N {V^{(m)}} }\sum_{m = 1}^N {\sigma _{ij}^{( m )} } V^{( m )} \end{align} and the standard deviation of the microscopic stress distribution is weighted by the volume of the element: \begin{equation} \label{eq4} SD = \sqrt {\frac{\sum_{m = 1}^N {V^{(m)}(\sigma _{ij}^{(m)} - \overline \sigma _{ij} } )^2}{\sum_{m = 1}^N {V^{(m)}} }}, \end{equation} where $V^{(m)}$ is the volume of element $m$, $N$ is the total number of elements, $\sigma _{ij} $ is the Cauchy stress tensor, and the total strain tensor is decomposed into elastic and plastic components. It should be noted that the stress and strain tensors, $\sigma _{ij} $ and $\varepsilon _{ij} $, are obtained at the centroid of each element. \subsection{Boundary conditions}\label{subsec:2.3} The experimental procedure to obtain a macroscopic stress-strain response for porous materials is typically a uniaxial compression test, which has been performed for porous titanium \cite{Davis, Shenetal}. Our numerical study thus focuses on the study of the responses of the titanium foam under uniaxial proportional compressive loading. As discussed in section~\ref{sec:1}, it is possible to obtain accurate estimates of the material properties with relative small RVEs under mixed BCs \cite{HazanovandHuet}. Two types of mixed BCs were designed in the numerical simulation as described in Table~\ref{tab1} in which $x_{2}$ direction is the loading direction. For boundary condition 1 (BC1), uniform displacements are imposed on the faces perpendicular to loading direction without friction; other faces parallel to loading direction are traction free. This boundary condition is used to simulate the experimental setup in the mechanical test. For boundary condition 2 (BC2), the same conditions are given for the faces perpendicular to loading direction; however the faces parallel to loading direction remain straight and parallel during deformation to simulate an interior domain to be compatible with the surrounding material. Therefore, BC2 is periodic mechanically, but not microstructurally as the pores on opposing edges are not continuous. The boundaries for 2D models are the same as those applied to the left, right, top, and bottom faces of the 3D models to impose uniaxial compressive loading. Note that 2D models are restricted to uniaxial and biaxial load cases. \section{Results and discussion}\label{sec:3} \subsection{Macroscopic response}\label{subsec:3.1} The uniaxial stress-strain responses predicted by the sixteen individual 2D simulations (eight models deformed along $x_{1}$ and $x_{2}$ directions) of the 2D-1000, 2D-2000, 2D-3000, and 2D-4000 models under the two mixed BCs are shown in Figure~\ref{fig6}. These collected results reveal that as 2D window size increases from $1000\,\mu$m to $4000\,\mu$m, the dispersion of the predicted stress-strain curves for each BC decreases. The range of the two groups of stress-strain curves associated with BC1 and BC2 also become smaller. The average responses of the sixteen individual 2D simulations under the two BCs converge to a common value as the window size increases as shown in Figure~\ref{fig7}. The convergence trend of the two BCs suggests BC1 is slightly underrestrictive which is similar to UTBC, and BC2 is slightly overrestrictive which is similar to UDBC. The asymptotic relationship of the average elastic modulus and overall stress at 1{\%} strain corresponding to the length scale of the window can be seen in Figure~\ref{fig8}. The further overall averages of the total thirty-two responses (eight models deformed along $x_{1}$ and $x_{2}$ directions under BC1 and BC2 ) for each model group are shown in Figure~\ref{fig9} and the curves are nearly indistinguishable. These observations in the figures are confirmed in Table~\ref{tab2} which shows the 2D model predictions for elastic modulus and overall stress at 1{\%} uniaxial strain. The table shows that under BC1 and BC2 the standard deviation of elastic modulus decreases from 2.34{\%} and 2.12{\%} for the 2D-1000 models to 0.45{\%} and 0.44{\%} for the 2D-4000 models and for the overall stress from 3.69{\%} and 3.0{\%} to 1.55{\%} and 1.44{\%}, respectively. For all the responses for the two BCs (thirty-two responses), the standard deviation decreases from 2.29{\%} to 0.48{\%} for modulus and 4.23{\%} to 1.7{\%} for overall stress. As the scatter of the data decreases, the overall averages (BC1{\_}2avgs) of each model group under the two BCs remain almost constant. For example, only 0.71{\%} difference exists between the averages of the overall stresses of the 2D-1000 and 2D-4000 models. Based on the observations above, we believe that the averages over the individual simulations under the two BCs should be close to the ``exact'' solution. This hypothesis originates from the homogenization convergence concept that the estimates under different BCs converge to the effective response as the calculated domain approach to a Hill's RVE for any locally heterogeneous but globally homogeneous microstructures \cite{HazanovandHuet, Hollister, Huet, Jiang2001a, Ostoja-Starzewski, Pecullan, Sab, Terada}. Therefore, we conclude that a number of small models can obtain convergent results equivalent to larger models. At the same time, since the averaged stress-strain curves associated with BC1 and BC2 approach the ``exact'' solution which lies in between, the error of the predicted result of each BC to the effective one can be estimated. For example, since 1.7{\%} difference exists between the average stresses at 1{\%} strain for the 2D-4000 models under the two BCs as shown in Figure~\ref{fig8}, the error of the result from each BC to the ``exact'' solution should be less than 1.7{\%}, and the error of the average of BC1 and BC2 (BC1{\_}2avg in Table~\ref{tab2}) less than 0.85{\%}. This method suggests that by selecting a certain number of models associated with two selected slightly underrestrictive and overrestrictive boundary conditions providing opposite bias of the result, the convergent results and degree of accuracy can be estimated. The fast homogenization achieved here can be credited to the fact that all the selected models have fixed porosity and are subjected to the two mixed BCs. These conditions are in contrast to the work of \cite{Kanit}. In their work on the determination of the RVE size, small models have large variance of volume fractions which leads to large variance of the prediction results. They showed that bias exists for the average of the predictions of a number of small domains associated with one periodic BC and the bias decreases to zero as the size of analyzed domains reach a certain level. This one sided approach can be eliminated only by increasing domain size. The study can also be generalized to 3D FE analysis. The uniaxial stress-strain responses predicted by the twelve individual 3D simulations (four models deformed along the three perpendicular directions) of the 3D-170 and 3D-340 models under the two mixed BCs are shown in Figure~\ref{fig10}. Similar to the 2D simulations, the scatter of the curves decreases as the size of the 3D models increases. The average stress-strain curves over the individual simulations associated with BC1 and BC2 for 3D-340 models are closer than 3D-170 models as shown in Figure~\ref{fig11}. The overall averages of the stress-strain responses for each model group of both the BCs (BC1{\_}2avg) are plotted in Figure~\ref{fig12} and the two curves almost superpose. The results for elastic modulus and overall average stress at 1{\%} uniaxial strain for the 3D models are shown in Table~\ref{tab3}. For all the response for the two BCs, i.e., twenty-four responses, the standard deviation for modulus decreases from 3.05{\%} for 3D-170 models to 0.97{\%} for 3D-340 models, and 5.95{\%} to 2.47{\%} for overall average stress. The two mixed BCs offer biased predictions as in the 2D models. Therefore, the discussions and conclusions based on the 2D FE analyses above still hold true for the 3D FE study. A detailed comparison between 2D and 3D FE simulations was conducted in our companion work \cite{ShenandBrinson}. It was found that the macroscopic responses predicted by the 3D models are in reasonable agreement with the experimental and theoretical results. The macroscopic plastic responses predicted by 2D models are lower than the 3D predictions while the elastic responses are close. 2D models overpredict the probability of high Von~Mises stress and equivalent plastic strain and therefore overestimate the failure probability for porous materials. It should be noted that the current BCs for the model groups are only appropriate for the simulation of titanium foam at low porosities (less than 25{\%}) under uniaxial loading conditions. However, this methodology elucidates a method to find relatively small RVEs for heterogeneous materials, especially for the difficult case of high contrast properties between the phases. As the material is foamed to high porosity, pores are connected to form larger pores complicating the microstructure as shown in Figure~\ref{fig5}. The uniaxial stress-strain responses predicted by the 3D model with higher porosity of 48{\%} (the 3D model deformed along $x_{1, }x_{2}$ and $x_{3}$ directions) under the two mixed BCs are shown in Figure~\ref{fig13}. The results indicate a much larger scatter than the 3D-340 models which have the same model size but lower porosity. The convergent response prediction for titanium foam at high porosity is therefore quite challenging and will be the focus of our future work. It can be expected that more analyzed domains are needed to obtain convergent result while using the boundary conditions providing biased error. \subsection{Microscopic response}\label{subsec:3.2} The microscopic field variable distributions are very important for failure analysis because failure is a local event rather than a volume averaged event. The Von~Mises stress distributions in the matrices of all the 2D models deformed along $x_{1}$ and $x_{2}$ directions under 1{\%} macroscopic strain are plotted in Figure~\ref{fig14}. Since BC1 is less restrictive than BC2, the distributions of Von~Mises stress under BC1 are broader than those under BC2. Similar to the observation regarding the macroscopic responses, as the model size increases, the curves become less dispersed. Both the mean value and standard deviation of the Von~Mises stress distribution converge with increasing model size. For example, the maximum difference of the mean values of Von~Mises stress under BC2 is 8.5{\%} for 2D-1000 models and 4.2{\%} for 2D-4000 models. The average distribution curves of the sixteen individual simulations of 2D-1000, 2D-2000, 2D-3000, and 2D-4000 models associated with the two mixed BCs converge to a common mean value and a common standard deviation as shown in Figure~\ref{fig15}. The difference between the averaged mean values under BC1 and BC2 for 2D-1000 models is 5.9{\%} and for 2D-4000 models is 1.8{\%}. This common mean value and common standard deviation of the Von~Mises stress distribution can be clearly seen in Figure~\ref{fig16} in which the further overall averages of the total thirty-two distribution curves of all the 2D models are plotted. Although not shown here, the equivalent plastic strain distribution follows the similar convergence trend. To obtain a clear view of the influence of the BCs on the microscopic distributions for models with different sizes, the equivalent plastic strain distributions predicted by one of the 2D-1000 and 2D-4000 models under 1{\%} far field strain are plotted in Figure~\ref{fig17}. In the 2D-1000 model, the plastic strain distributions are influenced by BCs such that differences exist even in the middle of the analyzed domain. In the 2D-4000 models, differences in the plastic strain distributions away from the boundaries become negligible. However, the deficiency of the individual small 2D models can be compensated by averaging results of more models associated with the two BCs. The accumulated frequencies of the Von~Mises stress of the 2D models are plotted in Figure~\ref{fig18}. These values were obtained by calculating the percentage of all the matrix elements exceeding a certain value in each model group. As is evident in Figure~\ref{fig18}, the results of the Von~Mises stress distribution of the 2D models under BC1 and BC2 also converge in the same manner as the macroscopic stress-strain responses. With this frequency averaged over all the 2D simulations under the two BCs, relatively accurate results can be obtained by the 2D-1000 model which is very close to the prediction of 2D-4000 models as shown in Figure~\ref{fig19}. Therefore the convergence discussion on the macroscopic responses should be still valid for the microscopic variable statistic distributions. To generalize the microscopic study to 3D FE analysis, the corresponding results of Von~Mises stress predicted by the twelve individual 3D simulations (four models deformed along the three perpendicular directions) of the 3D-170 and 3D-340 models under the two mixed BCs are shown in Figure~\ref{fig20}. Similar convergence trend can be observed as the size of the 3D models increases. The average responses over the individual simulations associated with the two BCs are shown in Figure~\ref{fig21}. The difference between the mean Von~Mises stress at 1{\%} strain for the 3D-170 models under BC 1 and BC 2 is 8.5{\%} and for 3D-340 models 3.6{\%}. The overall averages of the total responses for each 3D model group are shown in Figure~\ref{fig22}. Although the curves are not as close as 2D models, the convergence trend is similar to the 2D study. The accumulated frequencies of the Von~Mises stress of the 3D models are plotted in Figure~\ref{fig23}. The frequency averaged over all the 3D simulation under the two BCs is shown in Figure~\ref{fig24}. Similar convergence trend can also been observed. It can be expected that if the 3D models increase to $640\,\mu$m, the individual stress and strain curves would become more convergent. However, due to the computational limitation, it is not practical to analyze the model with $640\,\mu$m side length. Since the two BCs are bounds for the properties, the ``exact'' effective response should lie in between. In other words, since the difference between the mean Von~Mises stress at 1{\%} strain for the 3D-340 models under BC1 and BC2 is 3.6{\%}, the average should be within 3.6{\%} of the ``exact'' solution. It can be found that for microscopic response, the 3D models predict higher mean Von~Mises stress than 2D models but relatively uniform distributions with smaller standard deviation. For example, the distribution of the Von~Mises stress is $357.7\pm 46.7$\,MPa for 3D-340 models and $298\pm 95.5$\,MPa for 2D-4000 models. An extensive comparison between the stress and strain distribution has been studied in our companion work \cite{ShenandBrinson}. \section{Conclusions}\label{sec:4} An approach to determine RVEs of porous titanium has been presented in this paper. The method adopts the RVE concept of Kanit et al.\ according to which the RVEs can be small domains as long as the average of the small domains provides an unbiased result \cite{Kanit}. Since the estimates associated uniform traction and displacement BCs provide the lower and higher bounds for the effective properties providing the largest bias, two mixed boundaries conditions are therefore designed to obtain results with a smaller bias. Four groups of 2D models and two groups of 3D models with various sizes and fixed porosity have been constructed based on a simulated microstructure of an experimental material. As the length scale of the model increase, the individual responses of models become less dispersed. At the same time, the averages of the estimates associated with the two mixed BCs show opposite bias and converge to a common value. The error can be estimated and the relatively convergent result can be obtained from small models by averaging the responses associated with the two mixed BCs. This method developed here can be used to simulate microstructures of real materials. Although only a special case for titanium foam at low porosity under uniaxial loading has been studied, this method elucidates a way to study other heterogeneous materials and higher volume fraction of pores/inclusions. By choosing appropriate boundary conditions providing biased error, convergent result can be achieved by averaging the results of a number of small analyzed domains while the number of the domains is dependent upon the nature of the microstructure, the designed boundary, and the domain size. \section*{Acknowledgements} The authors acknowledge the financial support of the National Science Foundation through grant number DMR-0108342 and thank Professor Dunand's group for providing titanium foam and images. \bibliographystyle{plain} \bibliography{main} \newpage \begin{figure} \centerline{\includegraphics{..//main/fig-1a}} \centerline{\includegraphics{..//main/fig-1b}} \caption{Optical micrograph of metallographic cross-section for titanium foam with (a)~14.7\% and (b)~50\% porosity} \label{fig1} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-2a.eps}} \centerline{\includegraphics{..//main/fig-2b.eps}} \end{figure} \begin{figure} \centerline{\includegraphics{..//main/fig-2c.eps}} \caption{(a)~The cubic domain of the 3D-MP-70271 model with side length of $4000\,\mu$m and a series of window sizes of $1000\,\mu$m, $2000\,\mu$m, $3000\,\mu$m, and $4000\,\mu$m, (b)~Pore area fraction as a function of $x_3$ position for the series of windows, (c)~Distributions of pore area fraction of each window series} \label{fig2} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-3}} \caption{Geometry of (a) 2D-1000, (b) 2D-2000, (c) 2D-3000, and (d) 2D-4000 model} \label{fig3} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-4a}} \centerline{\includegraphics{..//main/fig-4b}} \caption{Geometry of (a) 3D-170 and (b) 3D-340 model} \label{fig4} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-5}} \caption{Geometry of the 3D model with porosity of 48\% in a $340\,\mu$m cube} \label{fig5} \end{figure} \newpage \ \centerline{\includegraphics{..//main/fig-6a}} \newpage \ \centerline{\includegraphics{..//main/fig-6b}} \newpage \ \centerline{\includegraphics{..//main/fig-6c}} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-6d}} \caption{Macroscopic stress-strain responses predicted by the sixteen individual simulations of (a)~2D-1000, (b)~2D-2000, (c)~2D-3000, and (d)~2D-4000 models associated with the two mixed boundary conditions} \label{fig6} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-7}} \caption{Macroscopic stress-strain responses averaged over the sixteen individual simulations of 2D-1000, 2D-2000, 2D-3000, and 2D-4000 models associated with the two mixed boundary conditions} \label{fig7} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-8a}} \centerline{\includegraphics{..//main/fig-8b}} \caption{(a)~Apparent modulus and (b)~overall average stress at 1\% far-field strain associated with the two mixed boundary conditions as a function of the ratio of the window size to pore diameter} \label{fig8} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-9}} \caption{Macroscopic stress-strain responses averaged over the thirty-two individual simulations of 2D-1000, 2D-2000, 2D-3000, and 2D-4000 models} \label{fig9} \end{figure} \newpage \ \centerline{\includegraphics{..//main/fig-10a}} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-10b}} \caption{Macroscopic stress-strain responses predicted by the twelve individual simulations of (a)~3D-170 and (b)~3D-340 models associated with the two mixed boundary conditions} \label{fig10} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-11}} \caption{Macroscopic stress-strain responses averaged over the twelve individual simulations of 3D-170 and 3D-340 models associated with the two mixed boundary conditions} \label{fig11} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-12}} \caption{Macroscopic stress-strain responses averaged over the twenty-four individual simulations of 3D-170 and 3D-340 models} \label{fig12} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-13}} \caption{Macroscopic stress-strain responses predicted by the 3D model of 48\% porosity associated with the two mixed boundary conditions} \label{fig13} \end{figure} \newpage \ \centerline{\includegraphics{..//main/fig-14a}} \newpage \ \centerline{\includegraphics{..//main/fig-14b}} \newpage \ \centerline{\includegraphics{..//main/fig-14c}} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-14d}} \caption{Von Mises stress distributions predicted by the sixteen individual simulations of (a)~2D-1000, (b)~2D-2000, (c)~2D-3000, and (d)~2D-4000 models associated with the two mixed boundary conditions} \label{fig14} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-15}} \caption{Von Mises stress distributions averaged over the sixteen individual simulations of 2D-1000, 2D-2000, 2D-3000, and 2D-4000 models associated with the two mixed boundary conditions} \label{fig15} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-16}} \caption{Von Mises stress distributions averaged over the thirty-two individual simulations of 2D-1000, 2D-2000, 2D-3000, and 2D-4000 models} \label{fig16} \end{figure} \begin{figure} \centerline{\includegraphics{..//main/fig-17}} \caption{Equivalent plastic strain distribution predicted by (a)~one of 2D-1000 model under BC1, (b)~the 2D-1000 model under BC2, (c)~one of 2D-4000 model under BC1, and (d)~the 2D-4000 model under BC2} \label{fig17} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-18}} \caption{Accumulative frequency of Von Mises stress exceeding a certain value predicted by 2D-1000, 2D-2000, 2D-3000, and 2D-4000 models associated with the two mixed boundary conditions} \label{fig18} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-19}} \caption{Accumulative frequency of Von Mises stress exceeding a certain value, averaged over the two mixed boundary conditions predicted by 2D-1000, 2D-2000, 2D-3000, and 2D-4000 models} \label{fig19} \end{figure} \newpage \ \centerline{\includegraphics{..//main/fig-20a}} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-20b}} \caption{Von Mises stress distribution predicted by the twelve individual simulations of (a)~3D-170 and (d)~3D-340 models associated with the two mixed boundary conditions} \label{fig20} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-21}} \caption{Von Mises stress distributions averaged over the twelve individual simulations of 3D-170 and 3D-340 models associated with the two mixed boundary conditions} \label{fig21} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-22}} \caption{Von Mises stress distributions averaged over the twenty-four individual simulations of 3D-170 and 3D-340 models} \label{fig22} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-23}} \caption{Accumulative frequency of Von Mises stress exceeding a certain value predicted by 3D-170 and 3D-340 models associated with the two mixed boundary conditions} \label{fig23} \end{figure} \newpage \ \begin{figure} \centerline{\includegraphics{..//main/fig-24}} \caption{Accumulative frequency of Von Mises stress exceeding a certain value, averaged over the two mixed boundary conditions predicted by 3D-170 and 3D-340 models} \label{fig24} \end{figure} \newpage \ \centerline{\includegraphics{last}} \begin{table} \begin{tabular} {|p{80pt}|p{125pt}|p{125pt}|} \hline \textbf{Model Face}& \textbf{Boundary Condition 1}& \textbf{Boundary Condition 2} \\ \hline Top ($x_2 = a$)& $u_{2} = -0.01a$;\quad $t_{1}=t_{3} = 0$& $u_{2} = -0.01a$;\quad $t_{1}=t_{3} = 0$ \\ \hline Bottom ($x_2 = 0)$& $u_{2} = 0$;\quad $t_{1}=t_{3} = 0$& $u_{2} = 0$;\quad $t_{1}=t_{3} = 0$ \\ \hline Left ($x_1 = 0$)& $t_{1}=t_{2}=t_{3} = 0$& $u_{1} = -\bar {u_1 }$;\quad $t_{2}=t_{3} = 0$ \\ \hline Right ($x_1 = a$)& $t_{1}=t_{2}=t_{3} = 0$& $u_{1}=\bar{u_1 } $;\quad $t_{2}=t_{3} = 0$ \\ \hline Front ($x_3 = a$)& $t_{1}=t_{2}=t_{3} = 0$& $u_{3}=\bar{u}_3 $;\quad $t_{1}=t_{2} = 0$ \\ \hline Back ($x_3 = 0$)& $t_{1}=t_{2}=t_{3} = 0$& $u_{3} = -\bar{u}_3 $;\quad $t_{1}=t_{2} = 0$ \\ \hline \end{tabular} \vspace{6pt} \caption{Boundary conditions ($\bar{u}_1 $ and $\bar{u}_3$ are displacements of the four side faces upon loading on the top)} \label{tab1} \end{table} \begin{table} \begin{tabular} {|p{37pt}|p{42pt}|p{40pt}|p{31pt}|p{51pt}|p{35pt}|} \hline \multicolumn{2}{|p{79pt}|}{\raisebox{-1.50ex}[0cm][0cm]{}}& \multicolumn{2}{|p{83pt}|}{$E$\,(GPa)} & \multicolumn{2}{|p{87pt}|}{$\sigma $\,(MPa)} \\ \cline{3-6} \multicolumn{2}{|p{139pt}|}{} & AVG.& SD{\%}& AVG.& SD{\%} \\ \hline \raisebox{-3.00ex}[0cm][0cm]{2D-1000}& BC1& 78.06& 2.34& 241.37& 3.69 \\ \cline{2-6} & BC2& 79.08& 2.12& 254.30& 3.00 \\ \cline{2-6} & \textbf{\textit{BC1{\_}2avg}}& \textbf{\textit{78.57}}& \textbf{\textit{2.29}}& \textbf{\textit{247.84}}& \textbf{\textit{4.23}} \\ \hline \raisebox{-3.00ex}[0cm][0cm]{2D-2000}& BC1& 77.19& 0.97& 242.59& 2.61 \\ \cline{2-6} & BC2& 77.92& 0.91& 250.1& 2.10 \\ \cline{2-6} & \textbf{\textit{BC1{\_}2avg}}& \textbf{\textit{77.56}}& \textbf{\textit{1.04}}& \textbf{\textit{246.27}}& \textbf{\textit{2.79}} \\ \hline \raisebox{-3.00ex}[0cm][0cm]{2D-3000}& BC1& 77.67& 0.54& 243.72& 1.92 \\ \cline{2-6} & BC2& 78.03& 0.51& 249.45& 1.29 \\ \cline{2-6} & \textbf{\textit{BC1{\_}2avg}}& \textbf{\textit{77.85}}& \textbf{\textit{0.57}}& \textbf{\textit{246.59}}& \textbf{\textit{1.99}} \\ \hline \raisebox{-3.00ex}[0cm][0cm]{2D-4000}& BC1& 77.69& 0.45& 244.04& 1.55 \\ \cline{2-6} & BC2& 78& 0.44& 248.15& 1.44 \\ \cline{2-6} & \textbf{\textit{BC1{\_}2avg}}& \textbf{\textit{77.85}}& \textbf{\textit{0.48}}& \textbf{\textit{246.09}}& \textbf{\textit{1.7}} \\ \hline \end{tabular} \vspace{6pt} \caption{2D model predictions for elastic modulus and overall stress in loading direction at 1{\%} uniaxial strain} \label{tab2} \end{table} \begin{table} \begin{tabular} {|p{37pt}|p{32pt}|p{45pt}|p{41pt}|p{51pt}|p{35pt}|} \hline \multicolumn{2}{|p{79pt}|}{\raisebox{-1.50ex}[0cm][0cm]{}}& \multicolumn{2}{|p{77pt}|}{$E$\,(GPa)} & \multicolumn{2}{|p{77pt}|}{$\sigma $\,(MPa)} \\ \cline{3-6} \multicolumn{2}{|p{129pt}|}{} & AVG.& SD{\%}& AVG.& SD{\%} \\ \hline \raisebox{-3.00ex}[0cm][0cm]{3D-170}& BC1& 80.36& 3.26& 266.73& 5.93 \\ \cline{2-6} & BC2& 81.93& 2.62& 287.91& 2.97 \\ \cline{2-6} & \textbf{\textit{BC1{\_}2avg}}& \textbf{\textit{81.15}}& \textbf{\textit{3.05}}& \textbf{\textit{277.32}}& \textbf{\textit{5.95}} \\ \hline \raisebox{-3.00ex}[0cm][0cm]{3D-340}& BC1& 80.61& 0.91& 275.63& 1.89 \\ \cline{2-6} & BC2& 81.42& 0.78& 286.66& 0.96 \\ \cline{2-6} & \textbf{\textit{BC1{\_}2avg}}& \textbf{\textit{81.01}}& \textbf{\textit{0.97}}& \textbf{\textit{281.14}}& \textbf{\textit{2.47}} \\ \hline \end{tabular} \vspace{6pt} \caption{3D model predictions for elastic modulus and overall stress in loading direction at 1{\%} uniaxial strain} \label{tab3} \end{table}