Finite element modeling of effective properties of nanoporous thermoelastic composites with surface effects

Finite element modeling of effective properties of nanoporous thermoelastic composites with surface effects VII International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2017 M. Papadrakakis, E. O˜nate and B. Schrefler (Eds) FINITE ELEMENT MODELING OF EFFECTIVE PROPERTIES OF NANOPOROUS THERMOELASTIC COMPOSITES WITH SURFACE EFFECTS ANDREY V. NASEDKIN, ANNA A. NASEDKINA AND ALEXANDR S. KORNIEVSKY Institute of Mathematics, Mechanics & Computer Science Southern Federal University Miltchakova str. 8a, 344090 Rostov on Don, Russia e-mail: [email protected] Key words: Porous Thermoelastic Composite, Nanosized Pore, Stochastic Porosity, Surface Stress, Highly Thermal Conducting Interface, Effective Moduli, Finite Element Method Abstract. This investigation concerns to the determination of the material properties of nanoscale thermoelastic composites of an arbitrary anisotropy class with stochastically distributed porosity. In order to take into account nanoscale level at the borders between material and pores, the GurtinMurdoch model of surface stresses and the highly conducting model are used. Finite element package ANSYS was used to simulate representative volume and to calculate the effective material properties. This approach is based on the theory of effective moduli of composite mechanics, modeling of representative volumes and the finite element method. Here, the contact boundaries between material and pores were covered by the surface membrane elastic and thermal shell elements in order to take the surface effects into account. 1 INTRODUCTION As it is well known from experiments, a scale effect can be observed for nanoscale bodies, which results in the change of effective stiffness and other material moduli compared to the corresponding macroscale bodies. Among various approaches that explain this phenomenon, the models of theory of elasticity with surface stresses are widely used now. The idea of surface stresses in solids has been formulated long time ago [27]. However, significant development of this idea was done later in [10, 13, 26]. As it was shown further, the theory of surface stresses can be considered as a particular case of the models with imperfect interface boundaries. At present the theory of surface stresses, commonly referred to as the model of Gurtin– Murdoch, has been become widely used for describing scale effects at nanolevel, which can 1 1140 Andrey V. Nasedkin, Anna A. Nasedkina and Alexandr S. Kornievsky be seen, for example, from overviews given in [3, 28]. In a range of papers this theory was applied for modelling of thermoelastic nanoscale composites. For example, in [1, 2, 16, 17] thermomechanical properties of composites with spherical nanoinclusions (nanopores) and fiber nanocomposites were studied in the frames of the theory of thermal stresses with surface effects. The methodology of finite element approximations for thermoelastic materials with surface effects was demonstrated in [11]. Models of lowly and highly thermal conducting interfaces [18, 19] are well known for modeling of effective thermal conductivity of composites with imperfect interface boundaries. The model of high conductivity with continuous thermal field when passing through the phase interface is similar to the Gurtin-Murdoch model for elastic fields. The model of lowly conducting interface, which includes Kapitsa contact thermal resistance, allows discontinuous temperature field. Generalizations of these models for a more general case of thermoelastic interface boundaries were presented [14], and a related review was given in [12]. The problems on the determination of effective thermal conductivity moduli for composite materials with imperfect boundaries, including microand nanoscale, were studied in [6, 15, 18, 19, 20, 30] and others. This paper considers anisotropic thermoelastic materials with randomly located nanopores. In order to take into account nanoscale level at the borders between material and pores, the Gurtin-Murdoch model of surface stresses and the highly conducting model are used. The paper is organized as follows. Section 2 presents the mathematical statement of a homogenization problem for two-phase composites with special conditions for stresses and heat flux discontinuities at the phase interfaces. Both composite phases are assumed to be anisotropic thermoelastic materials. The boundary value problem statements, their weak formulations and the resulting formulas for determination of the full set of effective constants for a two-phase composite with arbitrary types of phase anisotropy and surface properties are also described. We note that homogenization procedures for porous composites with surface stresses and heat fluxes can be regarded as special cases of the corresponding procedures for two-phase composites with imperfect interface boundaries under negligibly small stiffnesses and thermal stresses for nanoinclusions. The finite element approximations of the considered homogenization problems are given in Section 3. We note that homogenization problems for the composites under investigation can be solved with the help of known finite element software, using shell finite elements with membrane stresses options and plate thermal elements in order to take into account interphase surface stresses and heat fluxes. Following [25] Section 4 describes an implementation of the proposed approaches in the finite element software ANSYS. We suggest an algorithm for automatic determination of interphase boundaries and location of shell and plate elements on them, which will work for various sizes of representative volumes built in forms of cubic lattice of hexahedral thermoelastic and thermal finite elements. As an example, in [25] we consider the models of porous material of cubic crystal system for various values of surface moduli, porosity and number of pores. We note the influence of the magnitude of the area of interphase 2 1141 Andrey V. Nasedkin, Anna A. Nasedkina and Alexandr S. Kornievsky boundaries on the values of the effective moduli for porous material with nanoscale structure. 2 EFFECTIVE MODULI METHOD FOR HOMOGENIZATION OF THERMOELASTIC MIXED TWO-PHASE NANOCOMPOSITES Let Ω be a representative volume of thermoelastic two-phase composite body with nanodimensional inclusions; Ω = Ω(1) ∪Ω(2);Ω (1) is the volume occupied by the main materials of the first phase (matrix); Ω(2) is the set of the volumes occupied by the materials of the second phase (inclusions); Γ = ∂Ω is the external boundary of the volume Ω; Γsis the set of frontier surfaces of materials with different phases (Γs=∂Ω(1) ∩∂Ω(2)); νiare the components of the external unit normal vector νto the boundary, outward with respect to the region Ω(1) occupied by the material of the matrix; x={x1,x 2,x 3}is the vector of the spacial coordinates. We assume that the volumes Ω(1) and Ω(2) are filled with different anisotropic thermoelastic materials. Then in the framework of linear static theory of thermoelasticity we have the following system of differential equations σij,j =0,σ ij =cijklεkl −βij θ, εij =(ui,j +uj,i)/2,(1) qj,j =0,q i=−kijθ,j ,(2) where σij are the components of the stress tensor σ;εij are the components of the strain tensor ε;uiare the components of the displacement vector u;θis the temperature increment from natural state, cijkl are the components of the forth rank tensor of elastic stiffness moduli; βij are the thermal stress coefficients; qiare the components of the heat flux vector q;kij are the components of the tensor kof thermal conductivities; cijkl =c(m) ijkl, βij =β(m) ij ,σij =σ(m) ij ,x∈Ω(m), etc. In accordance with Gurtin–Murdoch model for surface stresses we will assume that on nanosized interphase boundaries Γsthe following equation is satisfied νi[σij]=∂s iσs ij,x∈Γs,(3) where [σij]=σ(1) ij −σ(2) ij ;∂s i=∂i−νi(νl∂l) are the components of the surface gradient operator; σs ij are the components of the surface stress tensor σs. We adopt that the surface stresses σs ij are related to the surface strains εs ij and the temperature θby the formulas σs ij =cs ijklεs kl −βs ij θ, ε s kl =(∂s kumAml +Akm∂s lum)/2,A ml =δml −nmnl,(4) where cs ijkl are the components of the forth rank tensor of elastic surface stiffness moduli; βs ij are the surface thermal stress coefficients; δml is the Kronecker delta. Similarly, for interphase boundaries Γs, we accept the equation of highly thermal conducting boundaries ni[qi]=∂s iqs i,q s i=−ks ij∂s jθ, x∈Γs,(5) 3 1142 Andrey V. Nasedkin, Anna A. Nasedkina and Alexandr S. Kornievsky where ks ij are the surface thermal conductivities. Setting the appropriate boundary conditions at Γ = ∂Ω, we can find the solutions of the problems (1)–(5) for heterogeneous medium in the representative volume Ω. Then the comparison of the solution characteristics averaged over Ω (such as stresses, heat flux, etc.) with analogous values for homogeneous comparison medium will permit to determine the effective moduli for the composite material. We note that for anisotropic media in order to determine the full set of the effective moduli it is necessary to solve several problems of the considered types for different boundary conditions. Here the main difficulties consist in the choice of the representative volume and boundary problems for the heterogeneous medium and the comparison medium, as well as the technologies for solving the problems for heterogeneous media. According to the previously developed methods of modeling the thermoelastic composite materials of ordinary sizes [23, 24], we consider analogous approaches for the problems of thermoelasticity with surface effects [25]. For thermoelastic homogeneous comparison medium we adopt that the same equations (1)–(5) are satisfied with constant moduli ceff ijkl,βeff ij ,keff ij , which are to be determined. Note that thermal problem (2), (5) is independent, and so the moduli ceff ijkl,βeff ij and the moduli keff ij can be found from separate problems. For determination of the moduli ceff ijkl,βeff ij let us assume that at the boundary Γ the following boundary conditions take place ul=xkε0kl,θ=θ0,x∈Γ,(6) where ε0kl =ε0lk,θ0are some values that do not depend on x. Then ul=xkε0kl,εkl =ε0kl, θ=θ0,σij =σ0ij =ceff ijklε0kl −βeff ij θ0will give the solution for the problem (1)–(6) in the volume Ω for the homogeneous comparison medium. Note, that for θ=θ0= const the equations (2), (5) are satisfied identically, because qi=q0i= 0, qs i= 0, and this pure thermal problem is not actually used for solution of mechanical problem with thermal stresses (1), (3), (4), (6). Let us solve now problem (1)–(6) for heterogeneous medium (or (1), (3), (4), (6) because for θ=θ0the equations (2), (5) are satisfied indentically) and assume that for this medium and for the comparison medium the averaged stresses are equal ⟨σij⟩=⟨σ0ij⟩, where hereinafter the angle brackets ⟨(...)⟩denote the averaged by the volume Ω and by the surfaces Γsvalues ⟨(...)⟩=1 |Ω|(∫Ω (...)dΩ+∫Γs (...)sdΓ).(7) Therefore we obtain that for the effective moduli of the composite the equation σij = σ0ij =ceff ijklε0kl −βeff ij θ0=⟨σij⟩is satisfied, where ε0kl and θ0are the given values from the boundary conditions (6). Hence, even in the assumption of the anisotropy of the general form for the comparison medium, all the stiffness moduli ceff ijkl and thermal stress coefficients βeff ij can be computed. Indeed, setting in (6) ε0kl =ε0(δkmδln +δlmδkn)/2, 4 1143 Andrey V. Nasedkin, Anna A. Nasedkina and Alexandr S. Kornievsky ε0= const, θ0= 0, where m,nare some fixed indexes, we get the computation formulas for the elastic moduli ceff ijmn:ceff ijmn =⟨σij⟩/ε0. If in (6) we set ε0kl = 0, θ0= 0, than from the boundary problem (1), (3), (4), (6) the thermal stress effective moduli can be obtained: βeff ij =−⟨σij⟩/θ0. As mentioned above, in order to determine the effective coefficients of the tensor kij it is sufficient to consider thermal conductivity equation (2) and interface relation (5). For the formulation of the corresponding boundary-value problem we adopt the boundary conditions in the following form θ=xjG0j,x∈Γ,(8) where G0jare the components of some constant vector that does not depend on x. It is obvious that θ=xjG0j,Gj=∂jθ,Gj=G0j,qi=q0i=−keff ij G0jwill give the solution of the problem (2), (5), (8) in the volume Ω for the homogeneous comparison medium. Having solved the problem (2), (5), (8) for heterogeneous medium, we can set that for this medium and for the comparison medium the averaged heat fluxes are equal ⟨qi⟩=⟨q0i⟩. As a result we get the equation for the effective moduli of the composite keff ij G0j=−⟨qi⟩, where G0jare the components of the vector known from the boundary conditions (8). Then for the comparison medium with anisotropy of general form it is not difficult to obtain computation formulas for thermal conductivity moduli keff ij . Indeed, setting in (8) G0j=G0δjl,G0= const, where l=1,2,3 is some fixed index, we get computation formulas for the moduli keff il :keff il =−⟨qi⟩/G0. The approaches described above are associated with the averaging of the moduli cijkl, βij,kij. Note that the boundary value problems (1), (3), (4), (6) and (2), (5), (8) differ from the usual problems of linear thermoelasticity by the presence of the interface boundary conditions (3)–(5) which are typical for the Gurtin–Murdoch model of surface stresses and the model of highly thermal conducting boundaries for nanosized bodies. For the numerical solution of the problems (1), (3), (4), (6) and (2), (5), (8) we derive their weak or generalized statements. Previously we introduce the space of the functions θand the vector functions u, defined on Ω. On the set of vector functions u∈C1which satisfy the first homogeneous boundary condition (6), i.e. ul= 0 on Γ, we introduce the scalar product (v,u)H1 u=∫Ω εij(v)εij(u)dΩ+∫Γs εs ij(v)εs ij(u)dΓ. The closure of this set of vector functions uin the norm generated by the indicated scalar product will be denoted by H1 u. For functions θ∈C1which satisfy second homogeneous boundary condition (6) or (8), i.e. θ= 0 on Γ, we introduce the scalar product (η,θ)H1 θ=∫Ω ∂iη∂ iθdΩ+∫Γs ∂s iη∂s iθdΓ. 5 1144 Andrey V. Nasedkin, Anna A. Nasedkina and Alexandr S. Kornievsky The closure of this set of functions φin the norm generated by the indicated scalar product will be denoted by H1 θ. In order to formulate the generalized or weak solution we scalar multiply the first equations (1) by arbitrary components viof the vector function v∈H1 u, and we multiply the first equation (2) by some function η∈H1 θ. By summing over iand by integrating the obtained equations on Ω, and by using the standard technique of the integration by parts with Eqs. (1)–(5), we obtain the following integral relations c(v,u)−β(v,θ)=0,(9) k(η,θ)=0,(10) where c(v,u)=cΩ(v,u)+cΓs(v,u),β(v,θ)=βΩ(v,θ)+βΓs(v,θ),(11) cΩ(v,u)=∫Ω cijklεij (v)εkl(u)dΩ,c Γs(v,u)=∫Γs cs ijklεs ij(v)εs kl(u)dΓ,(12) βΩ(v,u)=∫Ω εij(v)βij θdΩ,β Γs(v,u)=∫Γs εs ij(v)βs ij θdΓ,(13) k(η,θ)=kΩ(η,θ)+kΓs(η,θ),(14) kΩ(η,θ)=∫Ω kij∂iη∂ jθdΩ,k Γs(η,θ)=∫Γs ks ij∂s iη∂s jθdΓ.(15) Further, we present the solution {u,θ}of the problem (1), (3), (4), (6) in the form uj=udl +ubl,θ=θ0,(16) where udl satisfies homogeneous boundary mechanical conditions and ad hoc fitted functions ubl satisfy the inhomogeneous boundary conditions on Γ, i.e. udl =0,u bl =xkε0kl,x∈Γ,(17) and therefore, ud∈H1 u. By using (16) we can rewrite Eq. (9) in the form c(v,ud)=Lu(v),L u(v)=β(v,θ 0)−c(v,ub).(18) Now we can define the generalized or weak solution of the problem with thermal stresses (1), (3), (4), (6) using introduced functional space. Namely, the functions u,θin the form (16), (17) are the weak solution of the problem (1), (3), (4), (6), if Eq. (18) with (11)–(13) is satisfied for ∀v∈H1 u. Analogously, we will find the solution θof the purely thermal problem (2), (5), (8) in the form θ=θd+θb,(19) 6 1145 Andrey V. Nasedkin, Anna A. Nasedkina and Alexandr S. Kornievsky where θdsatisfies homogeneous boundary thermal conditions and θbis known function, satisfying the inhomogeneous boundary conditions on Γ, i.e. θd∈H1 θ θd=0,θ b=xjG0j,x∈Γ.(20) By using (19) we can also rewrite (10) in the form k(η,θd)=Lθ(η),L θ(η)=−k(η,θb).(21) Then we can introduce the generalized or weak solution of the thermal problem (2), (5), (8) as the function θin the form (19), (20), for which Eq. (21) with (14), (15) is satisfied for ∀η∈H1 θ. So far we have been discussing the two-phase composites. However, we can note that the presented models also describe homogenization procedures for porous composites with surface effects, if we put the stiffness and thermal stresses moduli negligible, and set the thermal conductivities equal to the coefficient of thermal conductivity of air. 3 FINITE ELEMENT SOLUTION For solving problems (18) and (21) for thermoelastic body with surface effects in weak forms we will use classical finite element approximation techniques. Let Ωhbe a region of the corresponding finite element mesh composed of volume elements, Ωh≈Ω, Ωh= Ω(1) h∪Ω(2) h,Ω (j) h≈Ω(j),Ω h=∪kΩek, where Ωek is a separate volume finite element with number k. Let also Γs hbe a finite element mesh of surface elements conformable with the volume mesh Ωh,Γ s h=∂Ω(1) h∩∂Ω(2) h,Γ s h≈Γs,Γ s h=∪mΓs em,Γ s em is a separate surface finite element with number m, and the elements Γs em are the faces of the suitable volume elements Ωek located on the interface boundaries. We will use the classic Lagrangian or serendipity volume finite elements with nodal degrees of freedom of displacements and temperature. Note that due to the structure of surface mechanical and thermal fields (4), (5), for the elements Γs em we can use standard shell or plate elements with elastic membrane stresses options, i. e. only with nodal degrees of freedom of displacements, and the standard thermal shell elements. For these elements we can take a fictitious unit thickness so that the surface moduli from (4), (5) can be determined by the product of specially defined volume moduli and shell thickness. On these finite element meshes we will find the approximation to the weak solutions {uh≈u,θh≈θ}for static thermoelasic problem in the form uh(x)=N∗ u(x)·U,θ h(x)=N∗ θ(x)·Θ,(22) where N∗ uis the matrix of the shape functions for displacements, N∗ θis the row vector of the shape functions for temperature, U,Θare the global vectors of nodal displacements and temperature, respectively. Here, the surface shape functions are the reduction on the boundaries Γs hof the volume shape functions. 7 1146 Andrey V. Nasedkin, Anna A. Nasedkina and Alexandr S. Kornievsky According to conventional finite element technique, we approximate the continuous weak formulations of the thermoelasticity problems by the corresponding problems in finite-dimensional spaces. Substituting (22) and similar representations for projection functions into integral relations (9), (10) for Ωh, we obtain the following finite element system Kuu ·U−Kuθ ·Θ=0,(23) Kθθ ·Θ=0,(24) where Kuu =KuuΩ+KuuΓ,Kuθ =Kuθ Ω+Kuθ Γ,Kθθ =Kθθ Ω+Kθθ Γ,(25) KuuΩ=∫Ωh B∗ u·c·BudΩ,KuuΓ=∫Γs h Bs∗ u·cs·Bs udΓ,(26) Kuθ Ω=∫Ωh B∗ u·βN∗ θdΩ,Kuθ Γ=∫Γs h Bs∗ u·βsN∗ θdΓ,(27) Kθθ Ω=∫Ωh B∗ θ·k·BθdΩ,Kθθ Γ=∫Γs h Bs∗ θ·ks·Bs θdΓ,, (28) B(s) u=L(s)∗(∇)·A·N∗ u,B(s) θ=∇(s)N∗ θ,(29) L(s)∗(∇)=   ∂(s) 1000∂(s) 3∂(s) 2 0∂(s) 20∂(s) 30∂(s) 1 00∂(s) 3∂(s) 2∂(s) 10   ,∇(s)=     ∂(s) 1 ∂(s) 2 ∂(s) 3      .(30) In (22)-(30), we use the following vector-matrix notation: cis the 6×6 matrix of elastic moduli; cαζ =cijkl,α, ζ =1, ..., 6, i,j,k,l =1,2,3 with correspondence law α↔(ij), β↔(kl), 1 ↔(11), 2 ↔(22), 3 ↔(33), 4 ↔(23) ∼(32), 5 ↔(13) ∼(31), 6 ↔(12) ∼ (21); β={β11,β22,β33,β23,β13,β12};(...)∗is the transpose operation; and (...)·(...) is the scalar product operation. Also, we can represent the finite element solutions in the another variants considering the main boundary conditions: uh=udh +ubh,θh=θdh +θbh,udh ≈ud,udh =N∗ ud ·Ud, ubh ≈ub,ubh =N∗ ub ·Ub,θdh ≈θd,θdh =N∗ θd ·Θd,θbh ≈θb,θbh =N∗ θb ·Θb,Nu= {Nud,Nub},Nθ={Nθd,Nθb}, Kuu =[Kdd uu Kdb uu Kbd uu Kbb uu ],Kuθ =[Kdd uθ Kdb uθ Kbd uθ Kbb uθ ],Kθθ =[Kdd θθ Kdb θθ Kbd θθ Kbb θθ ], U={Ud Ub},Θ={Θd Θb}, where Ub,Θbare the vectors, known from the main boundary conditions. So, after using these expressions we can transform Eq. (23) for the problem with thermal stresses in the form 8 1147 Andrey V. Nasedkin, Anna A. Nasedkina and Alexandr S. Kornievsky Kdd uu ·Ud=Fd u,Fd u=Kdd uθ ·Θd+Kdb uθ ·Θb−Kdb uu ·Ub,(31) where the vector Θdis also known, and its component values are equal to θ0from (6). Analogously, for the thermal problem (21) we can rewrite Eq. (24) as the system relative to unknown vector Θd: Kdd θθ ·Θd=Fd θ,Fd θ=−Kdb θθ ·Θb.(32) Thus, the homogenizing problems for thermoelastic composite with surface stresses and with highly conducting porous boundaries can be solved by finite element approaches. The resulting finite element systems (31), (32) differ from similar systems for the bodies of usual sizes by the matrices KuuΓ,KuθΓ,KθθΓin (26)–(28). These matrices arise due to the surface mechanical and thermal effects. 4 DISCUSSION AND CONCLUSION For automated coating of internal boundaries of pores in the cubic representative volume the following algorithm was used [25]. At the beginning, as a result of the formation of the porous structure, the finite element mesh from octanodal cubic elements was created, some of which had the material properties of thermoelastic material, and the other part of the elements had the material properties of the pores (with negligible elastic stiffness moduli). Further, only the finite elements with thermoelastic material properties were selected. The resulting elements on the outer boundaries were covered by four nodal target contact elements. Then, the contact elements, which were located on the external surfaces of the full representative volume, were removed, and the remaining contact elements were replaced by the four nodal membrane elastic elements. As a result, all the facets of the contact of thermoelastic structural elements with pores were coated by membrane finite elements. The next step consisted in solving the static problems for obtained representative volume with the main boundary conditions which were conventional for effective moduli method. Further, in the ANSYS postprocessor the averaged stresses were calculated, both on the volume finite elements and on the surface finite elements. Finally, the effective moduli of porous composite with surface effects were calculated from the corresponding formulas of the effective moduli method by using the estimated average characteristics. In the results of computational experiments, the following features were observed [25]. If we compare two similar bodies, one of which has usual dimensions and the other is a nanoscale body, then for the nanosized body due to the surface stresses the effective stiffness will be greater than for the body with usual sizes. Furthermore, for the porous body of the usual size the effective elastic stiffness decreases with increasing porosity. Meanwhile, the effective stiffness of nanocomposite porous body with the same porosity may either decrease or increase depending on the values of surface moduli, dimensions and 9 1148