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 the mechanical behavior of a multiphase composite material, comprising the following technical schemes and 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 0K 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 in the multiphase composite material under the condition of 0K; 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 a 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).Math.\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 the elastic constants of single crystal materials 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 u represents Poisson's ratio; Poisson's ratio is approximately unchanged with temperature; the bulks modulus, the shear modulus and the Young's modulus of the body will change with the change of temperature, its' calculation method is described in step (2); 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, which is compliance coefficient S of the material; the compliance matrix S and the 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 0K; cell volumes of relaxed metal and ceramic phases in step 1.1 are respectively shrunk or expanded in the 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 the quasi-harmonic Debye model, to obtain the 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 total energy of the cell calculated 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.Math.T\left[\frac{9}{8}.Math.\frac{\theta }{T}+3.Math.\ln \left(1-e^{-\theta .Math.\text{/}.Math.T}\right)-D\left(\theta .Math.\text{/}.Math.T\right)\right]& \left(11\right)\end{array}$ wherein, D(θ/T) represents the Debye function, $D\left(\theta .Math.\text{/}.Math.T\right)=\frac{3}{{\left(\theta .Math.\text{/}.Math.T\right)}^3}.Math.{\int }_0^{\theta .Math.\text{/}.Math.T}.Math.\frac{x^3}{e^x-1}.Math.\mathrm{dx}.Math..Math.\left(x=\theta .Math.\text{/}.Math.T\right),$ n represents the number of atoms in each cell, k.sub.B represents Boltzmann constant, and the Debye temperature can be expressed as: $\begin{array}{cc}\theta ={\frac{\hslash }{k_B}\left[6.Math.{\pi }^2.Math.V^{1.Math.\text{/}.Math.2}.Math.n\right]}^{1.Math.\text{/}.Math.3}.Math.f\left(v\right).Math.\sqrt{\frac{B_s}{M}}& \left(12\right)\end{array}$ wherein, M represents the cell mass, υ represents the Poisson's ratio, ℏ represents the reduced Planck constant, and f(υ) represents the Poisson's ratio function $f\left(v\right)={\left\{{3\left[2.Math.{\left(\frac{2}{3}.Math.\frac{1+v}{1-2.Math.v}\right)}^{3.Math.\text{/}.Math.2}+{\left(\frac{1}{3}.Math.\frac{1+v}{1-v}\right)}^{3.Math.\text{/}.Math.2}\right]}^{-1}\right\}}^{1.Math.\text{/}.Math.3};$ B.sub.s represents the adiabatic elastic modulus, which can be obtained by using the 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.Math.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), the thermodynamic state of system can be obtained, that is, the universal equation of system state can be obtained, and the relationship between p, T and V of the system 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 0K calculated by formula (16) and the Young's modulus E at 0K 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.Math.\left(1+v\right)};$ 2.2 calculating the thermal expansion coefficients of metal and ceramic phases at different temperatures: on the other hand, heat capacity at constant pressure Cv and Gruneisen parameters γ can be described as: $\begin{array}{cc}{C}_V=3.Math.{\mathrm{nk}}_B\left[4.Math.D\left(\theta .Math.\text{/}.Math.T\right)-\frac{3.Math.\theta .Math.\text{/}.Math.T}{e^{\theta .Math.\text{/}.Math.T}-1}\right]& \left(17\right)\\ \gamma =-\frac{d.Math..Math.\ln .Math..Math.\theta }{d.Math..Math.\ln .Math..Math.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 .Math..Math.C_V}{B_T.Math.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, the single crystal model of metal phase in the multiphase composite is obtained by LAMMPS or other software, and then the polycrystalline model with different grain sizes (between 30 and 40 nm) within a relatively narrow range is constructed by Vironoi method in AtomEye or other software; 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 the grain size; in order to reduce the error, hardening coefficient of polycrystalline models of all grain sizes can be expressed as the average of hardening coefficient of several different small-grain sizes; for the yield strength c 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 the lattice friction resistance when a single dislocation moves, k is a constant, and d represents the desirable arbitrary grain size; in the fitting process, the calculation result of Peirls-Nabarro stress τ.sub.p is adopted, σ.sub.0≈τ.sub.p; the value of Peirls-Nabarro stress of the metal phase main slip system is used as Peirls-Nabarro stress, and the calculation formula is as follows: $\begin{array}{cc}{\tau }_p=\frac{2.Math.G\left(T\right)}{1-v}.Math.\exp \left[-\frac{2.Math.\pi .Math..Math.a}{\left(1-v\right).Math.b}\right]& \left(21\right)\end{array}$ wherein, υ 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 the calculated results in step 2.1 are taken; (4) building three-dimensional finite element model of real microstructure of multiphase composites: dual beam focused ion beam system is used to obtain composite materials cuboids with micron side length, for example, a cube with a side length of about 15 μm is scanned at an interval of 20 nm; after image stack is obtained through the focused ion beam experiment, the image processing is carried out by Avizo or other software, the process is as follows: first, the FIB Stack Wizard module in Avizo 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 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), 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. A multi-scale simulation method for the mechanical behavior of multiphase composite material in claim 1, comprising: In step (4), the process of constructing geometric model in finite element model of composite material microstructure is as follows: (a) extracting the coordinates of all pixels of metal phase in the composite material microstructure image by using Matlab or other software, which is the pixel coordinates of the metal phase; (b) in the finite element software, geometric model in finite element model of composite material microstructure image is established, which is the total geometric model in finite element model, it means both metal phase and ceramic phase are regarded as ceramic phase to establish geometric model; obtaining the coordinates of all elements in total geometric model; (c) according to the coordinates of all elements in total geometric model in (b) and the pixel coordinates of metal phase in step (a), the element that belongs to the metal phase in all elements of total geometric model is determined, and the 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. A 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 load to geometry model in finite element model of composite material microstructure and calculating the complex stress; further quantitatively analyzing stress-strain relationship, stress distribution and its evolution, plastic deformation and other mechanical behaviors of the composite under different temperature and complex stress conditions.

Description

BRIEF DESCRIPTION OF THE APPENDED DRAWINGS

[0052] FIG. 1 is the flow chart of multi-scale simulation model and calculation method

[0053] FIG. 2 is an example of a binarization microstructure image used to construct a geometric model in finite element model

[0054] FIG. 3 shows the three-dimensional microstructure of a real material reproduced in the embodiment

[0055] FIG. 4 shows the volume-energy relationship calculated by the embodiment

[0056] FIG. 5 shows the comparison of stress-strain response curve of simulation results of the embodiment and experimental results

[0057] FIG. 6 is the simulation result of stress distribution in the embodiment

PREFERRED EMBODIMENT

[0058] 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

[0059] 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.

[0060] 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.

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

[0061] 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.

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

[0062] 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.

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

[0063] The Peirls-Nabarro stress of metal phase Coat 300K, 700K and 1100K are calculated.

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

[0064] 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.

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

[0065] 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.

[0066] 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.