Multi-scale method for simulating mechanical behaviors of multiphase composite materials

Abstract

A computer simulation analysis method suitable for describing the mechanical behavior of multiphase composites based on the real microstructure of materials relates to a multidisciplinary field such as computational material science, simulation and high throughput calculation. Through the first-principles calculation under nano scale, the molecular dynamics simulation under micro scale, and the thermodynamic calculation under mesoscopic scale, various physical parameters needed for the finite element simulation under macro scale can be obtained, including the elastic and plastic physical parameters of each phase in the composite at different temperature and different grain sizes. Focused ion beam experiment and image processing are adopted to obtain real material microstructure. Through the parameter coupling and parameter transfer among the calculated results of various scales, combining the microstructure of the material, stress-strain relationship, stress distribution and its evolution law, plastic deformation and other mechanical behaviors of the multiphase composites under complex stress and different temperature can be simulated.

Claims

1. A multi-scale simulation method for mechanical behavior of a multiphase composite material, comprising the following steps: (1) first-principles calculation under nano scale: conducting first-principles calculation under nano scale to respectively obtain elastic properties of metal phase and ceramic phase of the multiphase composite material at temperature 0 K and an energy corresponding to a crystal structure of different volumes; calculation scheme is as follows: 1.1 calculating Young's modulus, bulk modulus, shear modulus and Poisson's ratio of metal and ceramic phases of the multiphase composite material under temperature 0 K; firstly, structural relaxation of a single crystal structure of metal and ceramic phase materials is carried out respectively; elastic constant C as matrix element of a stiffness matrix is calculated by using the crystal structure obtained after the structural relaxation; $\begin{array}{cc}\left(\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right)=\left(\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{12}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{13}& c_{23}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{14}& c_{24}& c_{34}& c_{44}& c_{45}& c_{46}\\ c_{15}& c_{25}& c_{35}& c_{45}& c_{55}& c_{56}\\ c_{16}& c_{26}& c_{36}& c_{46}& c_{56}& c_{66}\end{array}\right)\left(\begin{array}{c}{\mathrm{.Math.}}_1\\ {\mathrm{.Math.}}_2\\ {\mathrm{.Math.}}_3\\ {\mathrm{.Math.}}_4\\ {\mathrm{.Math.}}_5\\ {\mathrm{.Math.}}_6\end{array}\right)& \left(1\right)\end{array}$ wherein σ.sub.i ε.sub.i c.sub.ij (i,j=1,2,3,4,5,6) represent stress, strain and elastic constants of single crystal material respectively; then the calculated elastic constants are converted into Young's modulus, bulk modulus, shear modulus and Poisson's ratio by using Voigt-Reuss-Hill approximation;
B=(B.sub.v+B.sub.r)/2  (2)
G=(G.sub.v+G.sub.r)/2  (3)
E=9B.Math.G/(3.Math.B+G)  (4)
υ=0.5.Math.(3B−2G)/(3B+G)  (5) wherein B represents bulks modulus, G represents shear modulus, E represents Young's modulus and V represents Poisson's ratio; Poisson's ratio is approximately unchanged with temperature; the bulks modulus, the shear modulus and the Young's modulus will change with the change of temperature; in formula (2) and (3), subscript v and r represent the approximate calculation methods of Voigt and Reuss respectively;
B.sub.v=(c.sub.11+c.sub.22+c.sub.33)/9+2(c.sub.12+c.sub.23+c.sub.13)/9  (6)
G.sub.v=(c.sub.11+c.sub.22+c.sub.33−c.sub.12−c.sub.13−c.sub.23)/15+(c.sub.44+c.sub.55+c.sub.66)/5  (7)
B.sub.r=1/(s.sub.11+s.sub.22+s.sub.33+2s.sub.12+2s.sub.13+2s.sub.23)  (8)
G.sub.r=15/(4s.sub.11+4s.sub.22+4s.sub.33−4s.sub.12−4s.sub.13−4s.sub.23+3s.sub.44+3s.sub.55+3s6)  (9) wherein S.sub.ij (i,j=1,2,3,4,5,6) represents matrix element of a compliance matrix S, which is compliance coefficient of the material; the compliance matrix S and a stiffness matrix C are inverse to each other; 1.2 calculating volume-energy relations of the metal and ceramic phases in the multiphase composite material at temperature 0 K; the crystal structure obtained after structural relaxation of metal and ceramic phases in step 1.1 is shrunk or expanded in a range of −10% to +10%, and energies corresponding to different cell volumes are calculated; (2) thermodynamic calculation of mesoscopic scale: conducting thermodynamic calculation of mesoscopic scale, based on quasi-harmonic Debye model, to obtain elastic properties and thermal expansion coefficients of metal phase and ceramic phase in the multiphase composite material at different temperatures respectively; calculation scheme is as follows: 2.1 calculating the Young's modulus of metal and ceramic phase in the multiphase composite material at different temperatures: non-equilibrium Gibbs free energy can be described as:
G*=E(V)+pV+A.sub.vib(θ,T)  (10) wherein, E(V) represents energies corresponding to the crystal structure with different volumes in step 1.2, P and V are external pressure and cell volume respectively, T represents temperature, A.sub.vib represents vibrational Helmholtz free energy, θ represents Debye temperature; vibrational Helmholtz free energy A.sub.vib can be expressed as: $\begin{array}{cc}{A}_{\mathrm{vib}}\left(\theta ,T\right)={\mathrm{nk}}_{B}T\left[\frac{9}{8}\frac{\theta }{T}+3\ln \left(1-e^{-\theta \text{/}T}\right)-D\left(\theta \text{/}T\right)\right]& \left(11\right)\end{array}$ wherein, D(θ/T) represents the Debye function, $D\left(\theta \text{/}T\right)=\frac{3}{{\left(\theta \text{/}T\right)}^3}{\int }_0^{\theta \text{/}T}\frac{x^3}{e^x-1}\mathrm{dx}\left(x=\theta \text{/}T\right),$ n represents the number of atoms in each cell, k.sub.n, represents Boltzmann constant, and the Debye temperature can be expressed as: $\begin{array}{cc}\theta ={\frac{\hslash }{k_B}\left[6{\pi }^2{V}^{1\text{/}2}n\right]}^{1\text{/}3}f\left(v\right)\sqrt{\frac{B_s}{M}}& \left(12\right)\end{array}$ wherein, M represents cell mass, v represents the Poisson's ratio, h represents reduced Planck constant, and f(v) represents the Poisson's ratio function $f\left(v\right)={\left\{{3\left[2{\left(\frac{2}{3}\frac{1+v}{1-2v}\right)}^{3\text{/}2}+{\left(\frac{1}{3}\frac{1+v}{1-v}\right)}^{3\text{/}2}\right]}^{-1}\right\}}^{1\text{/}3};$ B.sub.s represents adiabatic elastic modulus, which can be obtained by using volume-energy relationship data obtained in step 1.2, a calculation method is as follows: $\begin{array}{cc}{B}_s=V\left(\frac{d^{2}E\left(V\right)}{{\mathrm{dV}}^2}\right)& \left(13\right)\end{array}$ for a given (p,T), non-equilibrium Gibbs free energy minimizes the volume: $\begin{array}{cc}{\left(\frac{\partial G^{*}}{\partial V}\right)}_{p,T}=0& \left(14\right)\end{array}$ according to equation (14), thermodynamic state of system can be obtained, an universal equation of system state can be obtained, and the relationship between p, T and V can be determined; further, the bulk modulus of a certain temperature can be expressed as: $\begin{array}{cc}{B}_T=-{V\left(\frac{\partial p}{\partial V}\right)}_T& \left(15\right)\end{array}$ temperature effect can be introduced into the bulk modulus by calculating formula (15), and the bulk modulus B.sub.T at a specific temperature can be obtained; by formula (16):
E.sub.T=3B.sub.T(1−υ)  (16) the Young's modulus at different temperatures can be obtained; the Young's modulus E.sub.T at temperature 0 K calculated by formula (16) and the Young's modulus E at temperature 0 K calculated in step 1.1 has a difference ΔE=E.sub.T,0−E; this difference is used as the zero correction term of Young's modulus of other arbitrary temperatures calculated by formula (16), and the Young's modulus of any temperature is calculated accordingly E(T)=E.sub.T−ΔE; then the shear modulus G(T) at a particular temperature is calculated by $G\left(T\right)=\frac{E\left(T\right)}{2\left(1+v\right)};$ 2.2 calculating the thermal expansion coefficients of metal and ceramic phases at different temperatures: heat capacity at constant pressure Cv and Gruneisen parameters γ can be described as: $\begin{array}{cc}{C}_V=3{\mathrm{nk}}_B\left[4D\left(\theta \text{/}T\right)-\frac{3\theta \text{/}T}{e^{\theta \text{/}T}-1}\right]& \left(17\right)\\ \gamma =-\frac{d\ln \theta }{d\ln V}& \left(18\right)\end{array}$ based on formula (15), (17) and (18), the coefficient of thermal expansion α(T) can be obtained: $\begin{array}{cc}\alpha \left(T\right)=\frac{\gamma C_V}{B_{T}V}& \left(19\right)\end{array}$ therefore, the temperature effect can be introduced into the coefficient of thermal expansion to calculate the coefficient of thermal expansion at different temperatures; (3) molecular dynamics simulation under micro scale: molecular dynamics simulation under micro scale can obtain the plastic properties of metal phase in the multiphase composites at a specific temperature (that is, different temperatures can be taken) and a specific grain size, plasticity of hard and brittle ceramic phase is negligible compared with that of metal phase; firstly, a single crystal model of metal phase in the multiphase composite is obtained, and then a polycrystalline model with different grain sizes within a relatively narrow range is constructed by Vironoi method; compression simulation is carried out for polycrystalline models with different grain sizes in a relatively narrow range at a specific temperature, and stress-strain curves are drawn; yield strength and hardening coefficient can be read from the stress-strain curves; it is assumed that the hardening coefficient is independent of grain size; in order to reduce the error, hardening coefficient of polycrystalline models of all grain sizes can be expressed as an average of hardening coefficient of several different small-grain sizes; for the yield strength σ.sub.s of any model with different grain size, it can be obtained by fitting the yield strength data of polycrystalline model in a relatively narrow range read from the stress-strain curve by using the Hall-Petch relation;
σ.sub.s−=σ.sub.0+kd.sup.−1/2  (20) wherein, σ.sub.s represents the yield strength, σ.sub.0 represents lattice friction resistance when a single dislocation moves, k is a constant, and d represents desirable arbitrary grain size; in the fitting process, a calculation result of Peirls-Nabarro stress τ.sub.p is adopted, σ.sub.0≈τ.sub.p; a value of Peirls-Nabarro stress of the metal phase main slip system is used as the Peirls-Nabarro stress, and a calculation formula is as follows: $\begin{array}{cc}{\tau }_p=\frac{2G\left(T\right)}{1-v}\exp \left[-\frac{2\pi a}{\left(1-v\right)b}\right]& \left(21\right)\end{array}$ wherein, v represents Poisson's ratio, and takes the calculated result in step 1.1; a represents plane spacing of the sliding plane, and b represents atomic spacing in the sliding direction, which can be calculated by using the relaxed crystal structure in step 1.1 through the geometric relationship; G(T) is the shear modulus at a specific temperature, and calculated results in step 2.1 are taken; (4) building three-dimensional finite element model of real microstructure of multiphase composites: a dual beam focused ion beam system is used to obtain composite materials cuboids with micron side length; after image stack is obtained through focused ion beam experiment, an image processing is carried out, the image processing is as follows: first, FIB Stack Wizard module is applied to align and crop the image, as well as correct the shearing and shadow of the image, then median filter and edge-preserving smoothing Gaussian filter are used for denoising, finally, the image is binarized by threshold segmentation to distinguish the metal phase and ceramic phase; therefore, composite microstructure images are obtained which can be used to construct a geometric model in finite element model; then constructing the geometric model; (5) simulating the mechanical behavior of multiphase composites by using finite element method according to the Poisson's ratio obtained from step (1), Young's modulus and thermal expansion coefficient of metal and ceramic phase at particular temperature obtained in step (2), yield strength and hardening coefficient of metal phase with particular grain size at specific temperature obtained in step (3), a parameter model in finite element model of the ceramic and metal phase at different temperature is built; reading the geometric model in finite element model of composite material structure constructed in step (4); endowing the parameter model of ceramic phase and metal phase with different temperature to the geometric model of composite material microstructure for calculation.

2. The multi-scale simulation method for the mechanical behavior of multiphase composite material in claim 1, comprising: in step (4), the process of constructing the geometric model in finite element model of composite material microstructure is as follows: (a) extracting coordinates of all pixels of metal phase in a composite material microstructure image, which is the pixel coordinates of the metal phase; (b) in the finite element software, the geometric model in finite element model of the composite material microstructure image is established, which is a total geometric model in finite element model, it means both metal phase and ceramic phase are regarded as ceramic phase to establish the geometric model; obtaining coordinates of all elements in the total geometric model; (c) according to the coordinates of all elements in the total geometric model in (b) and the pixel coordinates of metal phase in step (a), an element that belongs to the metal phase in all elements of the total geometric model is determined, and material properties of such elements are modified to the material properties of metal phase to complete construction of the geometric model in finite element model of composite material microstructure.

3. The multi-scale simulation method for the mechanical behavior of multiphase composite material in claim 1, comprising: in step (5), defining boundary conditions and constraints; applying a load to a geometry model in finite element model of composite material microstructure and calculating complex stress; further quantitatively analyzing stress-strain relationship, stress distribution and its evolution, plastic deformation and other mechanical behaviors of the multiphase composite material under different temperature and complex stress conditions.

Description

BRIEF DESCRIPTION OF THE APPENDED DRAWINGS

(1) FIG. 1 is the flow chart of multi-scale simulation model and calculation method

(2) FIG. 2 is an example of a binarization microstructure image used to construct a geometric model in finite element model

(3) FIG. 3 shows the three-dimensional microstructure of a real material reproduced in the embodiment

(4) FIG. 4 shows the volume-energy relationship calculated by the embodiment

(5) FIG. 5 shows the comparison of stress-strain response curve of simulation results of the embodiment and experimental results

(6) FIG. 6 is the simulation result of stress distribution in the embodiment

PREFERRED EMBODIMENT

(7) The invention is further described in combination with embodiment, but the invention is not limited to the following embodiment. The following embodiments apply the multi-scale simulation calculation method of the present disclosure to study the stress-strain response and the distribution of stress and its evolution under the interaction between thermal stress and compression force for a composite material WC—Co.

Embodiment 1

(8) As shown in FIG. 2, WC—Co image stack is obtained by focused ion beam experiment, and the images are binarized to obtain microstructure images which can be used to construct the geometric model in finite element. Three-dimensional finite element model of WC—Co with real microstructure is obtained by using Avizo software for image processing, the size of WC—Co finite element model is 1.77 μm (axial)×2.35 μm (horizontal)×0.40 μm (thick), which is shown in FIG. 3.

(9) The Young's modulus, bulk modulus, shear modulus, Poisson's ratio and volume-energy relations of metal phase Co and ceramic phase WC in WC—Co composites at 0K are obtained by first-principles calculation. The volume-energy relationship is shown in FIG. 4.

(10) TABLE-US-00001 E (GPa) B (GPa) G (GPa) υ WC 679.6 378.0 283.1 0.200 Co 276.1 205.6 108.2 0.276

(11) By using the above Young's modulus, bulk modulus, shear modulus, Poisson's ratio and volume-energy relationship data at 0K which are calculated from the first principles, the Young's modulus after zero correction and thermal expansion coefficient of metal phase Co and ceramic phase WC at 300K, 700K and 1100K are obtained through thermodynamic calculation.

(12) TABLE-US-00002 thermal expansion coefficient Young's modulus E(T) (GPa) α(T) (10.sup.−5/K) WC E(300) 662 α(300) 0.58 E(700) 614 α(700) 0.74 E(1100) 560 α(1100) 0.82 Co E(300) 261 α(300) 0.75 E(700) 220 α(700) 0.99 E(1100) 166 α(1100) 1.24

(13) Three polycrystalline models of metal phase Co with grain sizes of 34 nm, 36 nm and 38 nm under uniaxial compression stress are simulated by means of molecular dynamics simulation under the conditions of 300K, 700K and 1100K to obtain the stress-strain curves. Values of yield strength and hardening coefficient are read from the stress-strain curves.

(14) TABLE-US-00003 Yield Grain size of polycrystalline Temperature strength Hardening models (nm) (K) (GPa) coefficient (GPa) 34 300 4.20 174.0 700 3.52 128.0 1100 2.93 126.4 36 300 3.53 162.0 700 3.43 134.8 1100 2.59 145.9 38 300 3.48 170.0 700 3.69 132.9 1100 2.69 145.3

(15) The Peirls-Nabarro stress of metal phase Coat 300K, 700K and 1100K are calculated.

(16) TABLE-US-00004 Temperature (K) Peirls-Nabarro stress (GPa) 300 0.237 700 0.200 1100 0.151

(17) Using the yield strength values read from stress-strain curves and the calculation results of Peirls-Nabarro stress, yield strength of metal phase Co with a grain size of 400 nm at 300K, 700K and 1100K is obtained by fitting the Hall-Petch relation. Assuming that the hardening coefficient of metal phase Co is independent of grain size, the hardening coefficient of metal phase Co at 400 nm is the average value of hardening coefficient calculated at 34 nm, 36 nm and 38 nm.

(18) TABLE-US-00005 Temperature (K) Yield strength (GPa) Hardening coefficient (GPa) 300 1287 168.7 700 1204 131.9 1100 927 139.2

(19) Poisson's ratio of metal phase Co and ceramic phase WC, and Young's modulus, thermal expansion coefficient, yield strength, hardening coefficient of Co at 300K, 700K, and 1100K, as well as Young's modulus and thermal expansion coefficient of WC at 300K, 700K, and 1100K are input into the parameter model in finite element model. The geometric model constructed by using real microstructure obtained by focusing ion beam experiment and image processing is read by finite element model. Using free boundary condition, the symmetrical boundary condition is set on three adjacent surfaces and free boundary condition is set on the other three opposite surfaces. The model is set from 1100K cooling down to 300K to simulate the cooling process after sintering. Stress-free state is set at 1100K, anelasto-plastic model with time variation is constructed. In the uniaxial compression simulation of sintered multiphase structure, the loading process is from 0 MPa to −3000 MPa, with an increase of 100 MPa for each load step, and the unloading process is from −3000 MPa to 0 MPa, with a decrease of 100 MPa for each load step.

(20) The internal statistical stress-strain response of Co phase in compression process is calculated and compared with the test results of neutron diffraction method in the literature, which is shown in FIG. 5. The simulation results of this model show a similar trend to the stress-strain curves of the experimental results, which verifies the rationality of this model. On this basis, the stress-strain law, stress evolution law and plastic deformation behavior in compression process can be quantitatively characterized and analyzed based on real microstructure. FIG. 6 shows the distribution of stress before compression, after loading and after unloading. Then the relationship between microstructure-plastic deformation behavior-macroscopic properties is studied to provide reliable guidance for design and control of microstructure of multiphase materials.