A.Yu.SMOLIN*，E.V.SHILKOS.V.ASTAFUROVI.S.KONOVALENKO，S.P.BUYAKOVAc，S.G.PSAKHIE，c，d

aInstitute of Strength Physics and Materials Science SB RAS，Tomsk 634021，RussiabTomsk State University，Tomsk 634050，RussiacInstitute of High Technology Physics，Tomsk Polytechnic University，Tomsk 634050，Russia

dSkolkovo Institute of Science and Technology，Skolkovo 143025，Russia

Received 7 June 2014；revised 14 August 2014；accepted 25 August 2014 Available online 4 December 2014

Modelingmechanical behaviors of compositesw ith various ratios of matrix-inclusion properties usingmovable cellular automatonmethod

A.Yu.SMOLINa，b，*，E.V.SHILKOa，b，S.V.ASTAFUROVa，b，I.S.KONOVALENKOa，S.P.BUYAKOVAa，b，c，S.G.PSAKHIEa，c，d

aInstitute of Strength Physics and Materials Science SB RAS，Tomsk 634021，Russia
bTomsk State University，Tomsk 634050，Russia
cInstitute of High Technology Physics，Tomsk Polytechnic University，Tomsk 634050，Russia

dSkolkovo Institute of Science and Technology，Skolkovo 143025，Russia

Received 7 June 2014；revised 14 August 2014；accepted 25 August 2014 Available online 4 December 2014

Tw o classes of com positematerials are considered：classicalmetal-ceram ic composites w ith reinforcing hard inclusions as well as hard ceramicsmatrix with softgel inclusions.Movable cellular automatonmethod is used formodeling themechanical behaviors of such different heterogeneousmaterials.Themethod is based on particle approach and may be considered as a kind of discrete elementmethod.Themain feature of the method is the use of many-body forces of inter-element interaction w ithin the formalism of simply deformable element approximation.Itwas shown that the strength of reinforcing particles and thew idth of particle-binder interphase boundaries had determining influence on the service characteristics ofmetal-ceram ic composite.In particular，the increasing of strength of carbide inclusionsmay lead to significant increase in the strength and ultimate strain of compositematerial.On the examp le of porous zirconia ceram ics itwas show n that the change in themechanical properties of pore surface leads to the corresponding change in effective elasticmodulus and strength lim it of the ceramic sample.The less is the pore size，themore is thiseffect.The increase in the elastic propertiesof pore surfaceof ceram icsmay reduce its fracture energy.

Composites；Metal ceramics；Zirconia ceramics；Gel；Modeling；Movable cellular automata；Many-body interaction

1.In troduction

Compositesarew idely used in industry and nature.Usually the composites consist ofmatrix，which constitutes themain part of the material，and a large number of small sized inclusions.Inclusions，as a rule，have higher physical and mechanical properties and are designed to improve the useful properties of composites.Typical examples are dispersionreinforced materials and，in particular，metal-ceram ic composites［1］.The availability of hard inclusions allow s producing the ultra-hard m aterials for cutting tools，and the antifriction and protecting materials for aircraft building［2］，etc.In case of porousmaterials，the void“inclusions”provide high heat-insulating propertiesand low specific density［3-5］.

Compositesconsisting of hardmatriceswith soft inclusions are not practically used in industry.As an exception，we can refer to metal-ceramics for self-lubricated bearings，which are produced by sintering of mixture of iron and carbon powders followed by oilsaturation.Muchmore suchmaterials exist in nature，first of all as fluid-saturated geologicalmedia and bone tissues of anim als and man.

It has to be noted that the properties of a com posite are often determined by the propertiesof the componentsand the fraction of inclusions.They also depend on the size and shape of inclusions as well as the properties of matrix-inclusion interface in a complicatedway［6-8］.That iswhy theproblem of compositemodeling is very difficult and complex.

An important trend inmechanicsof composite is to develop and apply the numericalmethods to study the dynamics and peculiarities ofmechanical behavior of heterogeneousmaterials and structures under complex loading conditions.The advantages of numerical m odeling in comparison w ith analytical solutions are associated w ith the ability to consider complex geometry of the system，to take into account heterogeneous structure of the material as well as to study in details the dynam ics of fracture and the related processes of redistribution and dissipation of elastic energy in the surroundings of the crack.At the same time，these advantages make strong demands for the formalism of applied numerical methods.

The most of research on computational mechanics are performed using finite elementmethod.Thismethod belongs to continuum concept in mechanics，i.e.，it assumes that the state controlling variables distribute in space continuously.A t the same time formodeling the severe distortion and failure of a material，converting finite elements into particles ismore effective and promising［9］.In Ref.［10］，so-called meshless（or mesh free）method were used for this purpose（namely smoothed particle hydrodynam ics and generalized particle algorithms）.Strictly speaking，meshlessmethods are not pure particle ones.They justuse thescattered nodesassociated w ith the centers of finite volumes to discretize the continuum mechanics equations in space，i.e.，they are really based on continuum approach.Concerning the problem of modeling heterogeneous media，it should be noted that at least 3-5 meshless particles are required to represent a matrix-inclusion interface.A true particle method can describe this interface using just one particle.The principal difference of particlemethods from computationalmethods in continuum mechanics is to replace the continuous representation of amaterial w ith an ensemble of interacting particles（atmicro-，meso-ormacroscopic scale）or pointmasses（at the atom ic scale w ithin the framework ofmolecular dynam ics or Monte-Carlo method）.This in turn determ ines the difference in governing equations aswell.In particular，the conventional partial differential equations of continuum m echanics are replaced by Newton-Euler equationsgoverning themotion of discrete ensemble.Constitutive laws for the considered material in tensor form，which conventionally describes the relationship between local stress and strain or their time derivatives，therew ith are replaced by the potentials/forces of inter-particle interaction.One of the most important consequences of these features of particle-based methods is an inherent capability of the discrete objects（particles）to change their surroundings（interacting neighbors）.This featuremakes“discrete”methods extremely attractive for directmodeling of complicated fracture-related processes including multip le fracture accom panied by formation and m ixing of large number of fragments［11-17］.

The key points determining the behavior of an ensemble of simple deformable discrete elementsare the structural form of the expressions for central and tangential interaction between elements and the relationship between these expressions and constitutive law of the modeled material.A conventional approach is based on the use of pair-w ise elastic interaction forces which can be treated as springs between elements. Corresponding value of the spring stiffness is derived on the assumption that strain energy stored in a unit cellof deform ed element ensem ble is equal to the associated strain energy of the equivalently deformed continuum［18，19］.An approximation of pair-w ise interaction has a number of important limitations，amongwhich are the follow ing：i）maximum value of Poisson's ratio of element ensemble depending on packing of elements；ii）packing-related artificial anisotropy of mechanical response of the ensemble；iii）fundamental problems in correct simulation of irreversible strain accumulation in ductilematerials.

Our research shows thatmany of these problems can be solved by using many-particle interaction［20-22］.Note that the construction of such relations becom es possible due to using a hybrid com putational technique to com bine the mathematical formalisms of discrete elements and cellular automata［23］.Many-body formulation of inter-element forces is adopted from the Wiener-Rosenblueth model of cellular automaton interaction［23，24］.This hybrid technique is referred to as movable cellular automaton method（MCA）［21，22，25，27］.The proposed generalized expression for interelement forces is the ability to establish the relationship between vector parameters of the interaction and tensor parameters of the material constitutive law.It m akes possible to implementdifferentmodelsand criteria of elasticity，plasticity and fracturew ithin themathematical formalism of MCA.

This paper is devoted to the computational study of compositematerials consisting ofmatrix and equiaxial inclusions using MCA method.The main peculiarity of the considered composite is the ratio ofmechanical properties of thematrix and inclusions that varies in aw ide range.First，we study the classical dispersion-reinforcedmaterial based on the example of metal ceram ics.Then the porous ceram ics w ith soft gel inclusions are considered.In the last case，the inclusion rigidity ismuch less than them atrix one.Themain attention of the study is paid on the influence of properties and w idth of inter-phase boundaries on the effectivemechanical properties of the composites.

2.Method ofm ovable cellular automata

Within the frameof MCA，it isassumed thatanymaterial is composed by a certain amountof finitesizeelementary objects（automata）which interactamong each otherand can rotateand m ove from one location to another，thereby simulating a real deform ation process.The automaton motion is governed by the New ton-Euler equations

Note that the automata of the pairmay represent the parts of different bodies or of one consolidated body.Therefore its interaction is not always really contact one.That is why we put theword“contact”in quotationmarks.More of that，as it shown in Fig.1，the size of the automaton is characterized by one parameter di，but it does notmean that the shape of the automaton is spherical.Real shape of the automaton is determined by area of its“contacts”w ith neighbors.For example，ifwe use initial FCC packing，then the automata are shaped like a rhombic dodecahedron，but if we use cubic packing，then the automata are cube-shaped.

For locally isotropic media，the volum e-dependent component can be expressed in terms of the pressure Pjin thevolume of the neighboring automaton j as follows［21，22］.

where Sijis the area of interaction surface of automata i and j；and A is amaterial parameter.

Fig.1.Schematic diagram of determ ining of the spatial parameters of a pair of themovable cellular automata i and j.

The total force acting on automaton i can be represented as a sum of explicitly defined normal componentand tangential（shear）component

Using homogenization procedure for stress tensor in a particle described in Refs.［27，32］，the expression for the components of the average stress tensor in automaton i takes the form

whereαandβdenote the axes X，Y，Z of the laboratory coordinate system；Viis the current volume of automaton i；nij，αis theα-component of unit vector；and Fij，βisβ-component of the total force acting at the point of“contact”between automata i and j.

Tensor com ponents allow us to calculate the other tensor invariants in the automaton volume，in particular stress intensity

From Eqs.（1）-（3）it follow s that the specific form of the expressions foranddeterm ines the rheological behavior of amodelmedium.

For further convenience，the interaction parameters of movable cellular automata are considered in relative（specific）units.Thus，the central and tangential interactions of the automata i and j are characterized by the corresponding stressesηijand

Note that，in themostofpapersdevoted to thedescriptionand use of MCA，the equations and formulas of the method are w ritten for two-dimensional case.This study uses threedimensional version of the MCA method.That iswhy herein shear stressisavector in theplanewhich isnormal to n→ij.

To characterize the deformation of automaton i under the normal interaction of automata i and j，we can use the follow ing dimensionless parameter（normal strain）

In general case，each automaton of a pair represents differentmaterial，and the overlap of the pair is distributed between i th and j th automata

where symbolΔdenotes the incrementof a parameter per time stepΔt of numericalintegration of themotion equation（1）.The distribution ruleofstrain in thepair isintimatelyassociatedw ith the expression for computing the interaction forces of the automata.This expression for central interaction is sim ilar to Hooke's relations for diagonal stress tensor components

where K is the bulk modulus；G is the shear modulus of the material of i th automaton；and Piis the pressure of automaton i，which may be com puted using Eqs.（3）and（4）at previous time step or by predictor-corrector scheme.

To determ ine a parameter characterizing shear deformation in pair of automata i-j，we startw ith kinem atics formula for free motion of the pair as a rigid body

Besides such rotation of the pair as a whole（defined by the difference in translational velocities of the automata），each automaton rotates w ith its own rotational velocityω→i（Fig.2）.The difference betw een these rotational velocities produces a shear deform ation.Thus，the increm ent of shear deform ations of autom ata i and j per time stepΔt is defined by the relative tangential displacement at the contact pointdivided by the distance between the automata

Fig.2.Translational and rotational velocities of the automata.

The expression for tangential interaction of movable cellular automata is sim ilar to Hooke's relations for nondiagonal stress tensor components and is pure pairw ise.

The difference in automaton rotation leads also to the deformation of relative“bending”and“torsion”（only in 3D）of the pair（Fig.3）.It is obvious that the resistance to relative rotation in the pair cause the torque，which value is proportional to the difference between the automaton rotations

Eqs.（1）-（4），（6）-（9），（11）-（14）describe themechanical behavior of a linearly elastic body in the framework of MCA method.Note that Eqs.（8），（9），（12）and（13）are w ritten in increments，i.e.，in the hypoelastic form.In Ref.［33］it is shown that thismodel gives the same resultsas the numerical solving usual equation of continuum mechanics for isotropic linearly elastic medium by finite-difference method.That makes it possible to couple MCA method with the numerical methods of continuum mechanics.The results in Ref.［26］ show that the rotation allows themovable cellular automata correctly to describe the isotropic response ofmaterial.

We propose to use the plastic flow theory（namely von M ises model）for the description of elastoplastic behavior w ithin the MCA method.For this purposeWilkins'algorithm［34］is adapted to the MCA m ethod［27，40］.It is well known thatW ilkins'radial return algorithm consists in the solution of elastic problem in the increments and subsequent“drop”of components of deviator stress tensor Dαβ=σαβ-1/3σkkδαβto von M ises yield surface in the case that stress intensity exceeds it（Fig.4）

where M=σpl/σint，σintis stress intensity；andσplis current radius of von M ises yield circle.

This algorithm，as applied to the automaton i，can be w ritten in the following notation

By analogy w ith the elastic problem，the specific normal and tangential interaction forcesare corrected with the use of the current value of coefficient M by directly reforming Wilkins'algorithm for average stress（Eq.（15））

Fig.3.Deformation due to relative rotation between automata.

Fig.4.Schematic of functioning of W ilkins'radial return algorithm.

Thus，the rheological properties of thematerial of the automaton i are determined by defining a unified hardening curveis the average intensity of strain tensor that can be computed，similar to［27，32］）；this dependence is also named as automaton response function.

The equations of motion，i.e.，Eq.（1），for the system of movable cellular automata are numerically integrated with the use of velocity Verlet algorithm modified by introducing a predictor for estim ation ofat the current time step［35］.

A pairof elementsmightbe considered asa virtualbistable automaton having two stable states（bound and unbound），which perm its simulation of fracture and coupling of fragments（or crack healing）by MCA.These capabilitiesare taken into accountbymeansof corresponding change of the state of the pair of automata（Fig.5（a））.A fracture criterion depends on the physical mechanisms of material deformation.An important advantage of the formalism described above is that itmakes possible direct application of conventional fracture criteria（Huber-M ises-Hencky，Drucker-Prager，Mohr-Coulom b，Podgorski，etc.），which are w ritten in tensor form.

To use a conventional fracture criterion for breaking interautomaton bond，the local stress tensor at the area of interaction（“contact”area）of the considered pair i-j（hereinafter

In the local coordinate system X＇Y＇of the pair（Fig.5（b）），the componentsandfor the pair i-j are numerically equal to the specific forces of central（σij）and tangential（τij）interaction of the autom ata（these forces are applied to the“contact”area Sij）.The other componentsandof local stress tensor in the local coordinate system Xare defined on the basis of linear interpolation of the corresponding valuesfor the automata i and j to the area of interaction

whereσcis the corresponding threshold value for the pair（strength of cohesion/adhesion）；a=σc/σtis the ratio of compressive strength（σc）to tensile strength（σt）of the pair bond；andandare the corresponding invariantsof the stress tensor

When the explicit scheme of integration ofmotion equations isused，thevalueof time stepΔt is lim ited from above by a quantity associated w ith the time of sound propagation through the bulk of elem ent.Normally the time step isequal to or less than a quarter of this lim it.In such a situation the conventionalmodel of bond breaking during one time stepΔt is an idealized condition because it virtually suggests that the spatial separation of atom ic layers occurs uniform ly over the whole surface of interaction of elements.Herein the follow ing approach is suggested for more accurate description of the crack grow th dynam ics.It is assumed that the breaking of a bond（linked→unlinked transition of the pair state）is a process distributed in time and space.To express it numerically，we propose to use the dimensionless factorwhich has themeaning of the portion of linked part of the contact area Sij.In this case，the linked part of the contactarea in the pair i-j isexpressed as Thus， the dynamics of bond breaking is determined by the dependence（t），wheredecreases from the initial value 1（totally linked pair）to the final value 0（totally unlinked pair）. Depending on the automaton size and features of the internal structure of the material the stable or unstable crack grow th models can be applied to describe the pair bond breaking.

Fig.5.Schematic diagram of the interacting pair of automata i and j.

In the first model（stable crack grow th），the process of fracture develops in accordance with the predefined dependence

In thesecondmodel（unstablecrackgrow th），itisassumed that，if the fracture criterion threshold isexceeded，then acrack begins to grow spontaneously according to thesomehow specified law

where t is time.In the simplest case，this dependence can be considered asa constant，whichmeans that the crack advances through thepair interaction area at the constantvelocity Vcrack. The value of Vcrackis a predefined（model）parameter，which reflects the rheology of the interface material between the interacting elements.In particular，for brittlematerials，Vcrackcould be close to Raleigh wave speed，while for ductile materials，its value should obviously be significantly smaller.

The distinctive features of interaction of“unlinked”（i.e.，contacting）automata i and j，among others，are（i）lack of resistance to tension（a pair is considered as interacting one if σij≤0 only）；and（ii）limited magnitude of the force of tangential interaction.Maximum value of the tangential force（τij）allowed between“unlinked”automata is determined by the model of friction applied（for example，Amonton's law，m odel of Dieterich and so on）.

3.Study of par ticle-reinforced m etal-ceram ic com posites

Metal-ceramic composites（MCC）are the advanced representativesof dispersion-reinforcedmaterials，which have the enhanced values of mechanical and service characteristics，such as strength，stiffness-to-weight ratio，crack grow th resistance，wear resistance，fracture energy，ratio of thermal conductivity to thermalexpansion coefficient，thermalstability and so on.Thismakes them very attractive forw idespread use in various industries，such asmaterials for extreme operating conditions［2，36］.A t the present time，the working parts of m achines and mechanisms operating under the conditions of shock loading，abrasion，high temperature and corrosive environment are made mostly of metal-ceramic composites on the basisof very hard and refractory compounds（carbides，nitrides，carbonitrides）w ith metallic binder（nickel and iron alloys）［1，2，37］.These materials are fabricated from the powder m ixtures of the compounds by powder metallurgy methods［2，38］.Themechanical and physicalpropertiesof the sinteredmetal-ceramic compositearedeterm ined，in addition to phase composition，by a number of structural factors（volume fraction，dispersion，geometry and faultiness of reinforcing particles，structural-phase state ofmetallic binder，etc.［2，39，40］）.

Numerical simulation is of great importance for studying the mechanical properties of composite materials and their dependence on geometric and mechanical characteristics of inclusion-matrix interfaces.In this paper，MCA method is applied to study the influence of reinforced ceram ic inclusions and geometric（w idth）properties of inclusion-binder interphase boundaries on the peculiarities ofmechanical response of metal-ceramic composite under dynamic loading.TiC-particle-reinforced Ni-Cr matrix composite（50 vol.%of TiC inclusionsw ith average size of3μm）is chosen asamodel system classified as metal-ceram ic composite material in which the inclusions aremuch stiffer than the binder.

Note that，for correct simulation ofmechanical response of such complex heterogeneousmaterials，it is necessary to take into accountmain features of their internal structure.For this purpose，a two-dimensional structural and rheologicalmodel of MCC in the framework of MCA method was proposed in the paper.In thismodel，each of the composite constituents is simulated by an ensemble ofmovable cellular automataw ith appropriate rheological parameters（thus the cellular automaton simulatesa domain/fragmentof an inclusion，or a binder，or an interface zone）.As an example of MCA-based structural model，F(xiàn)ig.6（a）show s the structure of an idealized m etal--ceram ic composite w ith TiC particles having an equiaxed shape and the average size DTiCof 3μm.The components of themetal-ceramic compositeshave very different rheological characteristics（in case of considered composite these are elastic-brittle high-strength TiC inclusions and elastic-plastic NiCr binder）.For correctmodeling of deformation and fracture of such complex systems，themathematical formalism of constructing many-particle interaction forces for cellular automata（discrete elements）w ith different rheological characteristicswas developed［27，33，41］.

In the follow ing calculations，the plane stress approximation is used.The elastic constants and diagram of uniaxial loading of thematerial are used as the input parametersof the model to characterize the inter-element interaction（the mechanical response function of movable cellular automaton）. The response function of automaton modeling NiCr is considered as a stress-strain diagram w ith linear hardening（Curve 1 in Fig.6（b））.Thisdiagram isan approximation of the experimentaldiagram of uniaxial compression ofmacroscopic samp les of the alloy.The mechanical properties of the automata that simulate the high-strength brittle inclusions meet the real properties of TiC particles in the ceram ic phase（polycrystalline particles，Curve 2 in Fig.6（b））.Twoparameter Drucker-Prager fracture criterion in Eq.（17）is used to describe the fractures ofmetallic binder and carbide particles in the developed model［27，33，41］.Thus，in this work，thestrength characteristicsof the compositematerialare determined by the tensile（σt）and compressive（σc）strength of its components.

The developedmodelwas used for investigating the influence of the m ain features of internal structure on the basic characteristics ofmechanical behavior of the m etal-ceram ic composite.For this purpose，a three-pointbending test［42，43］of themodel samplesofMCCwas simulated.The structure ofthe model set-up is shown in Fig.7.Dynamic loading of 20μm cylindricalmandrel atconstant velocity Vload=0.5m/s wasapplied to the sample of 24×130μm in dimensions.The value of themaximal resistance of the sample to bending and the corresponding displacementof themandrel，aswell as the dynam ics and pattern of the sample fracturewere analyzed.

Fig.6.Themodelofmetal-ceram ic composite and the response functions ofmovable cellular automata of its components.M ain components of the structure are plastic nickel-chrom ium binder（1）and high-strength brittle titanium carbide inclusions（2）.

Fig.7.Scheme of the simulated three-point bending test of themodel composite sample.

It is well known that the strength of reinforcing particles largely determ ines the integral characteristics of the composite.Cracking of particles，containing significant（large）internal defects and damages，leads to the nucleation of internal cracks at lower applied loads and furthermore accelerates the crack grow th.As a consequence，such integrated mechanical characteristics of the composite as strength and deformation ability could substantially reduce.Therefore，the influence of the strength of the TiC particles，which characterizes their initial defect structure or density of damages，on the integral mechanical response of thematerialwasanalyzed in the paper. Itshould be noted that the compressive and tensile strengthsof brittle m aterials can differ by an order of magnitude（σc＞＞σt）.The physical reason of this difference is that the damages，which initially exist and emerge in the process of loading，function in different ways under the action of compressive and tensile stresses.In conditionsof compression（including uniaxial），the surfaces of a damage could occlude. This leads to their relative resistance to compression and shear（because of adhesion and dry friction of the damage surfaces）. A spatial diversity of the damage surfaces occurs under the conditions of tension.Since the edges of damages are the stress concentrators，the absence of plastic deformation provides the conditions for the fast grow ing of damages and cracks.Thus，the key strength characteristic，which determines the integral strength of TiC particles in the composite，is its tensile strength（σt）.In the paper，the value ofσtvaried from 700 MPa that corresponds to the tensile strength of the metallic binder NiCr to 2000MPa.The choice of this interval isdetermined by theknown factthat the tensile strength of real samples w ith typical dimensions of 1÷10 m icrons are generally 15÷50 tim es less than the theoretical strength of material（44÷88 GPa for TiC in the case of absence of large damages）.The value ofσcfor TiC particles in all calculations was equal to 10，000MPa.

Note that the so-called model of“narrow”particle-binder interfaces［7］is used in the calculation presented below.In this model，it is assumed that the width of the interphase boundary ismuch smaller than the size of movable cellular automaton.Moreover，the assumption of ideal adhesion of the composite components（TiC inclusions and NiCr binder）was used.This assumes that the strength characteristics of the interface are equal to the corresponding values for the least durable component of the composite material（in this case Ni-Cr alloy）.Use of this approach allows to study the influence of the strength characteristicsof carbide inclusionson the mechanical response of the composite in the raw.

Fig.8 shows the dynamic loading of the simulated composite samplesw ith different valuesof tensile strength of TiC inclusions.The diagram is plotted in terms of normalized stressσ/versus bending angle.Here the stressσis defined as resistance force of composite applied on mandrel divided by the area of the upper surface of the mandrel，=700MPa is the tensile strength of NiCr binder.It could be seen from Fig.8 that the strength and deformation characteristicsof the composite（Curves3 and 4）are decreased at tensile strength of TiC inclusions）less than 1100 MPa. In particular，the decrease into the value of tensile strength ofmetallic binder leads to a twofold decrease in ultimate strain.

As could be seen from Fig.8，the essential change in the mechanical parameters of the composite occurs in a fairly narrow strength interval of TiC inclusions and has a threshold character.Analysis of the fracture dynam ics of the model samples shows that this is due to the active involvement ofcracking mechanism of low strength＜1100 MPa）TiC particles in the processof composite fracture.Fig.9 shows the m ain stages of fracture of the composites w ith different strength valuesof TiC inclusions.Itcould be seen from Fig.9 that，in the composites w ith high-strength inclusions，the cracks are nucleated and propagated in thematrix by passing the reinforcing particles along the interphase boundaries=2000MPa，top row in Fig.9）.For the compositeswith particleshaving a tensile strength approaching to the threshold value，the cracks are generated and propagated in the binder，but they“cut”them=1000MPa，middle row in Fig.9）when they reach to TiC particles.Further decrease in the particle strength leads to a change in fracture mechanism of the com posite.As could be seen from Fig.9（bottom row），fracture beginsat the bottom partof the sample（in the area of maximum tensile stress）by means of nucleation of cracks in ceramic particles.Then，w ith the increase in theapplied strain，these cracks consequently connect w ith the cracks passing through thebinder into a singlemain crack.Thischange in the fracture mechanism is accompanied by the significant decrease in integral deformation characteristics of the composite.Thus，the strength of the reinforcing inclusions is one of themost important factors determ ining a number of service characteristicsofmetal-ceram ic composites，such as strength，critical deformation，fracture toughness，and others.

Fig.8.The loading diagrams（normalized stressσ/versusbending angle）for themodel samples of metal-ceramic composite w ith different strengthof TiC inclusions under three-pointbending test.

It should be noted that one of the key elements of the internal structure of the metal-ceramic composites is the interface between the particles of refractory compounds and the metallic binder.The change in the technological peculiaritiesofmetal-ceramic composite fabrication（in particular，applying additional heat treatmentof the composite）can vary the geometry（w idth）of interphase boundaries and，consequently，theirmechanical properties［7］.Therefore，the influence of this factor on the mechanical characteristics of MCC was investigated in the paper.The investigation was carried out using amesoscopicmodel of“w ide”interphase boundary（transition zone）.In thism odel，it isassumed that thew idth of the interface is comparable to or greater than the size of the movable cellular automaton［7］.Here the particle-binder interface is regarded as an area of variable composition of chemical elements（Ti，Ni，Cr，C）and modeled by several layers of cellular automata.In this area，the volume fractions of TiC and NiCr vary w ith a distance from the particle surface to thebulk of binderaccording to agiven law.This leads to an appropriate change in the physical andmechanical（including response function and strength parametersσcandσt）characteristics of these“transition”autom ata.This model can be efficiently used in numerical simulation of composites in the case of thew idth of the transition zone higher than the size of the cellular automaton.

The existing experimental data on the structure of particlebinder interfaces show that these areas contain a significant amount of secondary dispersed particles of titanium carbide（w ith the sizeof 50÷100 nm）［1，7，38］.Thevolume contentof this particles decreases w ith the increase in distance from particle surface to the bulk of binder.In thepaper，this factor is considered using so-called multiscale approach.In the framework of this approach，the elements of the internal structure of the composite at the mesoscopic scale（binder，primary particles TiC and w ide interphase boundaries）are taken into account explicitly.Taking into account structural elements of lower scales（in particular，the presence of secondary titanium carbide nanoparticles in the interphase boundaries）is carried outby considering thestructuralmodels of lower scales（in this case the submicron scale）w ith explicit prestoring concentrations，sizes and spatial distribution of nanoscale elements.On the base of analysis of simulation results ofmechanical tests of submicron-sized samples their integral m echanical（including rheological）properties are determined.In the subsequent，they are used as input data for cellular automatamodeling composite response at themesoscopic scale.Thus，using of themultiscale approach allows to obtain the dependence ofmechanical properties of the interphase boundaries on local concentration of the nanoparticles of the carbide phase.

The influence of thew idth of interphase boundaries on the integral mechanical characteristics of metal-ceram ic composite is investigated hereafter.Fig.10 show s the loading diagrams of the composite samp les w ith different w idth of particle-binder interfaces.It could be seen from Figs.10 and 11 that the increase in the w idth of the interphase boundaries leads to the increase in the strength of the composite，as well as the increase in the value of the ultimate strain（critical value of bending angle in the considered test）.As follows from the analysis of dependencies presented in Fig.11，the main effect of increasing thewidth of the interphase boundaries up to 1.6 microns is manifested in the increase in the critical strain of thematerial（up to 2 times）.

Fig.9.Dynam ics of fracture in metal-ceramic composite samples characterized by different strength values of TiC inclusions.Hereαis the bending angle.

Analysis of the results of computer simulation show that the increase in strength and value of ultim ate strain of the compositematerialw ith the increase in w idth of the interphase boundaries is due to the significant expansion of the region of the stress reducing from the high value in the reinforcing particles（which are the stress concentrators in the composite）to the significantly lower value in theplastic binder.VonM ises stress distributions for the sampleswhich are characterized by differentw idthsof interphase boundaries are shown in Fig.12. It can be seen from Fig.12 that the formation of w ide transition zones around the reinforcing particles，which are characterized by a smooth change of themechanical characteristics w ith the distance from the surface of the ceram ic inclusion to the bulk of the binder，leads to“sm earing”the stress field and，consequently，to reducing thestressgradientatthe interface.Thismeans that the formation ofw ide interphase boundaries among reinforcing particles and metallic binder interfaces in metal-ceramic composites provides a relatively low level of stress in the transition zones.This leads to increasing thedeformation ability and strength of themodified surface layers.

It is necessary to note that the results obtained herein by simulation are in good agreement w ith the existing experimental data.As pointed out in Ref.［40］，the constraint im posed on matrix p lastic deform ation by the ceram ic reinforcements induces large tensile hydrostatic stresses in the matrix.This enhances the load carried by the reinforcements and hence the composite flow stress，and also triggers the early development of internal damage in the form of particle fracture，interface decohesion，and/or matrix void grow th.New experimental techniques，such as automated serial section and X-ray computed microtomography（XCT），provide detailed information on the relationship between the damage nucleation or grow th and specific features of the three-dimensional m icrostructure［40，42］.The experimentsw ithm odelmaterials，presented in Ref.［40］，show that the damage of the model compositematerialm ade of a softmatrix ism ainly attributed to decohesion along theparticle/matrix interface，while for the materialw ith the same structure butw ith a hardermatrix，the damagemechanism changes to particle fracture.In Ref.［42］，for anotched glass fiber/epoxy cross-ply laminatesubjected to three-pointbending，the onsetand evolution of the damage in three dimensionswere studied by XCT.Itwas found that the damage began by formation of intra-ply cracks in the 90°plies followed by intra-ply cracking in the 0°plies.

Our simulation results also exp lain that such a perspective way can improve the m echanical characteristics of m etal--ceram ic composites if the high-strengthm ultiscale therm ally stable structures are formed in their surface layers under the impact of concentrated streams of charged particles.For example，the electron-beam irradiation on thesurface layersof metal-ceramic alloy containing 50 vol.%of TiC and 50 vol.% of NiCr leads to a considerable increase in themechanical and operating characteristics of thematerial［44，45］.The experimental results［44，45］showed that the increase in the operating characteristicsofmodified composite isdue to a number of factors.Among them are the significant increase in the w idth of the interphase boundary between the matrix and reinforcing particles and the formation of a significantamount of secondary nanosized particlesof TiC among these transition zones.Furthermore，the superfast heating of surface layer leads to the fragmentation and partial dissolution of most damaged carbide inclusions，and to the healing of small defectswhich are contained in the reinforcing TiC particles.The direct consequence of this is a significant increase in the strength（including tensile strength）of ceram ic inclusions and the integral operational properties of thematerial.

Fig.13 shows the change of temporal durability at cutting and the m icro-hardness of MCC TiC-NiCr after pulsed electron beam irradiation in nitrogen-containing gas discharge plasma.It can be seen from Fig.13 that the operating characteristicsof the composite can be considerably improved by the electron-beam treatment of the surface layers of the material.

4.Study of porous ceram icsw ith gel

The next part of the paper is devoted to modeling the porous ceram ics filled w ith gel，which matrix stiffness and strength are much higher than those of the inclusions. ZrO2（Y2O3）based ceram icswas used as amatrix［3，4］.There are three maxima in pore size distribution bar chart for this material：1μm；2μm（comparable w ith grain size）；6μm（cells）.Gel-forming composition GALKA，developed for enhancing oil recovery［46］，was used as the filler of the composite.In Ref.［47］，gel filling was performed through spontaneous soaking of ceramics samples by initial gelforming liquid composition during 24 h.Soaking occurs due to capillary effect.Then，thesoaked samplewasbaked to form the gel at 80°C for 1 h；the samp le cooled in a kiln.M echanical testing of the soaked ceram ics showed that their mechanical properties depends on the structure of the initial ceramics powder.Corresponding experimental data are shownin Fig.14.Thus，the effective elasticmodulus and strength of the ceramics obtained from nanocrystalline powders significantly decreased after gel filling（Fig.14（a）），while the properties of coarse-crystalline ceramics increased（Fig.14（b））［47］.Such strange effectmay be explained by changing the interface between thematrix and gel inclusions.

To ascertain the nature of changing mechanical properties of ceramics after filling its poresw ith gel，herewemodel the uniaxial compression of plane（p lain strain）ceram ics samp les w ith em pty pores as well as w ith pores filled w ith gel.It is assumed that fills only the 2μm and 6μm sized pores are filled w ith gel，which haveequiaxialshape.Thedimensionsof sample are 300×300μm.Response function of the automata modeling ceramics［3，4］corresponds to elastic-brittlematerial（like curve 2 in Fig.6（b）），response function of the automata modeling gel corresponds to elastic-plastic material w ith bilinear hardening（like curve 1 in Fig.6（b））.Large pores of the ceram icsaremodeled explicitly in theMCA model，while the 1μm and lesssized poresare taken into account implicitly（via the parameters of response function）.The fracture criterion of ceram ics automata is based on the threshold value of von M ises stress，but for pairs“gel-ceramics”a twoparametric criterion of Drucker-Prager is used.

Total porosity of themodeled samples is 0.30（this corresponds to 3D porosity about40%）.Eachmaxima of pore size distribution have fraction of 50%，25%and 25%of the total porosity.

Fig.12.Von Mises stress distributions in the central partof themodelmetal-ceramic composite samplesw ith“narrow”（0.1microns）and“w ide”（1micron）interphase boundaries（the value of bending angle is 3.5°）.

Fig.13.Increase in strength of MMC after electron-beam irradiation in nitrogen-containing gas-discharge p lasma.

The size of cellular automata corresponds to the averaged grain size，i.e.is equal to 2μm.Pore structure is built using two procedures.First procedure assumes only cells of 6μm sized pores；thesecond one generates2μm sized poresaswell as cells of pores.The pores are generated by removing the random ly selected single automata（2μm sized pores）and their six closed neighbors（cells of pores）from initial closed packing of automata.Thus，in the first procedure，the exp licit porosity is of 7.5%and the imp licit one is of 22.5%；in the second procedure，the corresponding valuesare 15%and 15%（Fig.15）.

Analysisof the calculation data shows that the filling pores w ith gel results in a slight increase in their elastic and strength properties.Young'smodulus，Poisson's ratio and strength for the samples in Fig.15（a）increases by 0.8%，2.4%and 3.7%，respectively；and for the samples in Fig.15（b），they increase by 1.5%，4.3%and 1.3%，respectively.Thus，these changesare sm aller for strength and larger for elastic properties.

The fracture patterns of the sim ulated samples，asnetworks of inter-automata bonds，are shown in Fig.16.In all the cases，the samples failsw ith the system ofmacrocracks propagating along the lines inclined to the direction of loading.For the samples generated using the first procedure（Fig.16（a）and（b）），the cracksare linearand wellexpressedwhen thevalueof explicit porosity is minimal（7.5%）.This is because the structure isclosed to homogeneousbody.For the sampleswith exp licit porosity equaling to 15%（Fig.16（c）and（d）），the cracks form a“blurred”picture and their path is curved and complicated which can be explained by large numberof stress concentrators in the sample.

Fig.14.Experimental data on the impact of gel filling on mechanical properties of ceram ics.

Fig.15.Initial structure of samples obtained w ith different procedures of pore generation.

In all the considered cases，the filling of porous samples w ith gel leads to the change of distribution of stress concentrators in the samples，in particular makes this distribution more homogenous.Thus the generation and propagation of cracks in thesamplesw ith filled poresdiffer significantly from those of the same sam ple w ith empty pores.A t that，the fracture pattern of the samp lew ith filled pores is closer to that of hom ogeneousm aterial［48，49］.

In themodeling mentioned above，itwas assumed that the gelwas fully filled in a pore and had perfectly contacted with ceramics surface.Itwas not possible to describe the experimental evidence in the framework of thatmodel［47］.Therefore another model，based on the assumption that the gel changes the local elastic and strength properties of the interface“gel-ceramics”，was developed at the next step.

The total porosity of the samples in thismodelwas varied from 0.15 to 0.27.Herewith the porosity C1corresponding to 1μm and less sized small poreswas equal to 0.1 for all the samples and taken into account implicitly.The size of large poreswas equal to 2μm or 6μm；the poreswere taken into account explicitly，the corresponding porosity C2was varied from 0.05 to 0.17.Sample dimension was300×300μm，and automaton sizewas2μm.Threshold value of von M ises stress was used as a fracture criterion.

Three groups of samp les w ith total porosity of 0.15，0.21 and 0.27 are considered.The porosities C2corresponding to large poresare equal to 0.05；0.11；and 0.17，thereby formingthree subgroups.Each subgroup contains six samples w ith individual spatial distribution of pores.

Fig.16.Networksof inter-automata bonds in the samplesgenerated using different procedures show ing fracture patterns in thematerial.

Fig.17.Initial structures of themodeled sampleswith different pore size D and propertiesof the“pore-matrix”interface for explicit porosity C2=0.11.

Fig.18.Normalized effective elastic modulus of ceramics samples w ith modified elastic property Elocalof“pore-matrix”interface for different pore size D versus fraction of pores w ith modified properties.

To model the different effect of a gel on the properties of the interface“gel-ceram ics”，wevary the elasticity（Elocal）and strength（σlocal）parameters of the response function of the corresponding automata.For pure ceram ics，these parameters are E0=98 GPa andσ0=730 MPa.Due to gel effect，the elastic properties was assum ed to be improved tw ice：Elocal=0.5E0=49 GPa or Elocal=2E0=196 GPa；the strength propertyσlocal=0.6σ0=438 MPa or σlocal=1.4σ0=1022 MPa（i.e.to be improved by±40%）. Typical sample structures for thismodel are shown in Fig.17.

Analysisof the calculated results shows that the change of local elastic properties for the“pore-matrix”interface has the most impact on the effective elastic modulus.Fig.18 shows the effectiveelasticmodulus Eeffof the ceram icssamplesw ith modified elastic property of“pore-matrix”interface，which is normalized to the m odulus Eeff_0of the corresponding pure ceramics samples versus fraction of pores w ith modified properties.One can see that changing Elocaltw ice leads to change Eeffby+10%/-7%for C2=0.05，and by+35%/ -25%for C2=0.17 for the samplesw ith poresize D=6μm. For the samplesw ith smaller pore size（D=2μm），thiseffect is stronger（due to larger value of specific surface of pores）：+25%/-23%for C2=0.05 and+85%/-48%for C2=0.17.

The local strength properties of the interface automata are changed to cause a variation of the effective strength of the composite and its fracture pattern；it does not influence on the effective elastic modulus.Fig.19 show s the effective strength limitσeffof the ceramics samples w ith modified strength propertyσlocalof“pore-matrix”interface，which isnormalized to the strength limitσeff_0of the corresponding pure ceram ics samples versus fraction of pores w ith modified properties. Thus，changingσlocalby±40%leads to changeσeffby+14%/ -22%for C2=0.05，and by+19%/-37%for C2=0.17 for the samplesw ith pore size D=6μm（Fig.19）.For the samples w ith smaller pore size（D=2μm），this effect is stronger：+27%/-32%for C2=0.05 and+40%/-41%for C2=0.17.

Analysis of specific fracture energy of themodel samples（Efr）shows that the increase in local strength propertyσlocalleads to the corresponding increase in Efr，whereas the increasing/decreasing of local elastic property Elocalleads to decreasing/increasing of Efr（Figs.20 and 21）.As described in the foregoing，this change is stronger for smaller pore size. Thus，the possible influence of gel soaking of ceramicson its dissipative properties is ambiguous，because the fracture energy is defined by both elastic and strength properties of the m aterial.In particular，if the soaking causes the increase in local elastic properties of pore surface，it may result in decreasing the fracture energy，but the increase in the strength properties always results in increasing the fracture energy.

Fig.19.Normalized effective strength lim itof ceram icssamp lesw ithmodified strength propertyσlocalof“pore-matrix”interface for different pore size D versus fraction of poresw ithmodified properties.

Fig.20.Normalized specific fracture energy of ceramics samples w ith modified elastic property Elocalof“pore-matrix”interface for different pore size D versus fraction of pores w ith modified properties.

Fig.21.Normalized specific fracture energy of ceramics samples w ith modified strength propertyσlocalof“pore-matrix”interface for different pore size D versus fraction of poresw ithmodified properties.

5.Conclusions

Based on movable cellular automaton method，themodels for computer-aided studying mechanical behavior of composites w ith special properties of inclusions and inclusionmatrix interphasewere developed.

The influences of strength of reinforcing ceram ic particles as well as the mechanical and geometric properties of interphase boundaries on the peculiarities ofmechanical response of TiC-reinforced Ni-Cr metal-ceram ic composite under dynam ic loading were studied using the proposed model.It was shown that the strength of reinforcing particles and the w idth of particle-binder interphase boundaries had a determining influence on the service characteristics of themetal--ceramic composite.In particular，the increase in strength of carbide inclusions may lead to a significant increase in strength and valueof ultimate strain of the compositematerial. And use of the technologies form ing thick interphase boundaries among the carbide inclusions and the matrix during production and processing of the composites allow to enhance their service characteristics.

Based on the example of porous zirconia ceramics，changing the elastic properties of the pore surface mainly leads to the corresponding change in effective elasticmodulus of the ceram ic samples.But changing the strength properties of the pore surface leads to the change in effective strength limit of the samples only.The less the pore size is，themore this effect is.It is interesting that the fracture energy of the ceram ics w ith increased elastic properties of the pore surface may reduce.This effectmay explain the peculiarities ofm echanical behavior of gel-filled zirconia ceramics.

Acknow ledgm en ts

The investigation has been carried out at partial financial supportof the Projects Nos.III.23.2.3（I.S.Konovalenko，S.P. Buyakova）and III.23.2.4（S.G.Psakhie）of the Basic Scientific Research Program of State Academ ies of Sciences for 2013-2020，the RFBR Project No.12-01-00805-a（A.Yu. Smolin，E.V.Shilko），and the grant No.14-19-00718 of the Russian Science Foundation（A.Yu.Smolin，E.V.Shilko，S.V. Astafurov）.

［1］Yu B，Ovcharenko VE，Psakhie SG，Lapshin OV.Electron-beam treatment of tungsten-free TiC/NiCr cermet II：structural transformations in the subsurface layer.JM ater Sci Technol 2006；22（4）：511-3.

［2］Chaw la N，Chaw la KK.Metalmatrix composites.New York：Springer；2006.

［3］Kulkov SN.Structure，phase composition and mechanical properties of ZrO2-based nanosystems.PhysM esomech 2008；11（1-2）：29-41.

［4］Buyakova SP，WeiH，Dunmy L，Haiyun C，Sablina TYu，Mel'nikov AG，Kul'kov SN.Mechanical behavior of porous zirconium dioxide under active deformation by compression.Tech Phys Lett 1999；25（9）：695-6.

［5］Shen J，Ren X.Experimental investigation on transm ission of stress waves in sandw ich samples made of foam concrete.Def Technol 2013；9（2）：110-4.

［6］Konovalenko IgS，Smolin AYu，Korostelev SYu，Psakhe SG.Dependence of themacroscopic elastic properties of porousmedia on the parameters of a stochastic spatial pore distribution.Tech Phys 2009；54（5）：758-61.

［7］Psakhie S，Ovcharenko V，Yu B，Shilko E，Astafurov S，Ivanov Yu，Byeli A，M okhovikov A.Influence of features of interphase boundaries on mechanical properties and fracture pattern in metal-ceram ic composites.JMater Sci Technol 2013；29（11）：1025-34.

［8］Dimaki AV，Shil'ko EV，Psakh'e SG.Simulating the propagation of exothermic reactions in heterogeneous media.Combust Explos Shock Waves 2005；41（2）：151-7.

［9］O~nate E，Rojek J.Combination of discrete element and finite element methods for dynamic analysis of geomechanics problems.Comput Method Appl Mech Eng 2004；193：3087-128.

［10］Johnson GR，Stryk RA.Conversion of 3D distorted elements into meshless particles during dynamic deformation.Int J Impact Eng 2003；28：947-66.

［11］Munjiza AA，Knight EE，Rougier E.Computationalmechanics of discontinua.Chichester：W iley；2012.

［12］Cundall PA.A discontinuous future for numerical modeling in geomechanics.Geotechn Eng 2001；149（1）：41-7.

［13］Psakhie SG，Smolin AYu，Tatarintsev EM.Discrete approach to study fracture energy absorption under dynam ic loading.Comput M ater Sci 2000；19：179-82.

［14］¨Osterle W，Dmitriev AI.Functionality of conventional brake friction materials-perceptions from findingsobserved atdifferent length scales. Wear 2011；271（9-10）：2198-207.

［15］Morris JP，Rubin MB，Blair SC，Glenn SA，Heuze FE.Simulations of underground structures subjected to dynamic loading using the distinct elementmethod.Eng Comput2004；21：384-408.

［16］Dimaki AV，Popov VL.Themethod of reduction of dimensionality and its application to simulation of elastomer friction under complex dynamic loads.PhysMesomech 2012；15（5-6）：319-23.

［17］Kuznetsov VP，Smolin IYu，Dmitriev AI，Konovalov DA，Makarov AV，Kiryakov AE，Yurovskikh AS.Finite element simulation of nanostructuring burnishing.Phys M esomech 2013；16（1）：62-72.

［18］Griffiths DV，Mustoe GGW.Modelling of elastic continua using a grillage of structural elements based on discrete element concepts.Int J Numer M ethods Eng 2001；50：1759-75.

［19］Wang G，Al-Ostaz A，Cheng AH-D，Mantena PR.Hybrid lattice particle modeling：theoretical considerations for a 2D elastic spring network for dynamic fracture simulations.Comput Mater Sci2009；44：112634.

［20］Psakhie SG，Ostermeyer GP，Dmitriev AI，Shilko EV，Smolin AYu，Korostelev SYu.Method ofmovable cellularautomata asa new trend of discrete computational mechanics.I.Theoretical description.Phys Mesomech 2000；3（2）：5-12.

［21］Psakhie SG，Horie Y，Ostermeyer G-P，Korostelev SYu，Smolin AYu，Shilko EV，et al.Movable cellular automatamethod for simulatingmaterials w ith mesostructure.Theor App l Fract Mech 2001；37（1-3）：311-34.

［22］Popov VL，Psakhie SG.Theoreticalprinciples ofmodeling elastoplastic media by movable cellular automatamethod.I.Homogeneousmedia. Phys M esomech 2001；4（1）：15-25.

［23］M ikhailov AS.Foundations of synergetics I.Distributed active systems. Berlin：Springer；1994.

［24］Wiener N，Rosenblueth A.Themathematical formulation of the problem of conduction of impulses in a network of connected excitable elements，specially in cardiac muscle.A rch Inst Cardiol M ex 1946；16（3-4）：205-65.

［25］Psakhie SG，Horie Y，Korostelev SYu，Smolin AYu，Dmitriev AI，Shilko EV，et al.Method of movable cellular automata as a tool for simulation w ithin the framework of mesomechanics.Russ Phys J 1995；38（11）：1157-68.

［26］Smolin AYu，Roman NV，Dobrynin SA，Psakhie SG.On rotation in the movable cellular automaton method.Phys Mesomech 2009；12（3-4）：124-9.

［27］Psakhie S，Shilko E，Smolin A，Astafurov S，Ovcharenko V.Developmentof a formalism ofmovable cellular automatonmethod for numerical modeling of fracture of heterogeneous elastic-plastic materials. Fratt Int Strutt 2013；24：26-59.

［28］Cundall PA，Strack ODL.A discrete numericalmodel for granular assemblies.Geotechnique 1979；29（1）：47-65.

［29］Jing L，Stephansson O.Fundamentals of discrete element method for rock engineering：theory and applications.Ox ford：Elsevier；2007.

［30］Sibille L，Nicot F，Donze FV，Darve F.Material instability in granular assemblies from fundamentally different models.Int J Numer Anal Methods Geomech 2007；31（3）：457-81.

［31］M artin CL，Bouvard D.Study of the cold compaction of composite powders by the discrete element method.Acta Mater 2003；51（2）：373-86.

［32］Potyondy DO，CundallPA.A bonded-particlemodel for rock.Int JRock M ech M in Sci 2004；41（8）：1329-64.

［33］Psakhie SG，Shilko EV，Smolin AYu，Dimaki AV，Dmitriev AI，Konovalenko IgS，Astafurov SV，Zavshek S.Approach to simulation of deformation and fracture of hierarchically organized heterogeneous media，including contrast media.Phys Mesomech 2011；14（5-6）：224-48.

［34］Wilkins ML.Computer simulation of dynamic phenomena.Berlin：Springer-Verlag；1999.

［35］Smolin AYu，Roman NV，Konovalenko IgS，Erem ina GM，Buyakova SP，Psakhie SG.3D simulation of dependence ofmechanical properties of porous ceramics on porosity.Eng FractMech 2014；130：96-115.

［36］Kainer KU.In：Kainer KU，editor.Metalmatrix composites：custommade materials for automotive and aerospace engineering.Weinheim：Wiley-VCH Verlag；2006.

［37］Pramanik A，Zhang LC，Arsecularatne JA.Prediction of cutting forcesin machining of metal matrix composites.Int J Mach Tool M anuf；46（14）：1795-1803.

［38］Vityaz PA，Grechikhin LI.Nanotechnology for producing titanium-based cermets.Phys Mesomech 2004；7（5-6）：51-6.

［39］Ayyar A，Chaw la N.M icrostructure-based modeling of crack grow th in particle reinforced composites.Compos Sci Technol 2006；66：1980-94.

［40］M ortensen A，Llorca J.M etalmatrix composites.Annu Rev M ater Res 2010；40：243-70.

［41］Psakhie SG，Horie Y，Shilko EV，Smolin AYu，Dmitriev AI，Astafurov SV.Discrete element approach to modeling heterogeneouselastic-plastic materials and media.Int J Terrasp Sci Eng 2011；1：93-125.

［42］Xing MZ，Wang YG，Jiang ZX.Dynamic fracture behaviors of selected aluminum alloys under three-point bending.Def Technol 2013；9（4）：193-200.

［43］Sket F，Seltzer R，Molina-A ldareguía JM，Gonz＇alez C，Llorca J.Determination of damagemicromechanisms and fracture resistance of glass fiber/epoxy cross-ply laminate by means of X-ray computed microtomography.Compos Sci Technol 2012；72：350-9.

［44］Ovcharenko VE，Yu B，Psakhie SG.Electron-beam treatment of tungsten-free TiC/NiCr cermet.I：influence of subsurface layermicrostructure on resistance to wear during cutting of metals.JMater Sci Technol 2005；21（3）：427-9.

［45］Ovcharenko VE，M okhovikov AA，Ignat'ev AS.Influence of surface nanostructure on the life of cermet in metal cutting.Steel Transl 2013；43（6）：348-50.

［46］Altunina LK，Kuvshinov VA，Stasieva LA.Thermoreversible polymer gels for enhanced oil recovery.Chem Sustainable Dev 2011；19：121-30.

［47］Developmentof scientific principles for the construction of a new class of contrast materials w ith nanocomposite gels as a functional filler. Tomsk：Institute of Petroleum Chemistry SB RAS；2013.Report on interdisciplinary integration projectNo 66 of the Siberian Branch of the Russian Academy of Sciences，［in Russian］.

［48］Skripnyak VA，Skripnyak EG，Kozulin AA，Skripnyak VV，Korobenkov MV.Computer simulation of the relation between mechanical behavior and structural evolution of oxide ceram ics under dynamic loading.Russ Phys J 2009；52（12）：1300-8.

［49］Kostandov YuA，Makarov PV，Eremin MO，Smolin IYu，Shipovskii IE. Fracture of compressed brittle bodies w ith a crack.Int JAppl Mech 2013；49（1）：95-101.

.Institute of Strength Physics and M aterials Science SB RAS，Tomsk 634021，Russia.

E-mail address：asmolin@ispms.tsc.ru（A.Yu.SMOLIN）.

Peer review under responsibility of China Ordnance Society.

http：//dx.doi.org/10.1016/j.dt.2014.08.005

2214-9147/Copyright?2014，China Ordnance Society.Production and hosting by Elsevier B.V.All rights reserved.

Copyright?2014，China Ordnance Society.Production and hosting by Elsevier B.V.A ll rights reserved.