- Email: [email protected]
Software for estimation of second order effective material properties of porous samples with geometrical and physical nonlinearity accounted for
Advances in Engineering Software 86 (2015) 80–84
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
Software for estimation of second order effective material properties of porous samples with geometrical and physical nonlinearity accounted for A.V. Vershinin a, V.A. Levin a,⇑, K.M. Zingerman b, A.M. Sboychakov c, M.Ya. Yakovlev c a b c
Lomonosov Moscow State University, Moscow, Russian Federation National Research Nuclear University MEPhI, Moscow, Russian Federation Fidesys LLC, Moscow, Russian Federation
a r t i c l e
i n f o
Article history: Received 29 August 2014 Received in revised form 10 March 2015 Accepted 11 April 2015
Keywords: Effective properties Porous material Finite strain Physical nonlinearity Finite-element method Parallel computations
a b s t r a c t A method and an algorithm for numerical estimation of effective mechanical properties of porous materials are presented. The effective properties are sought in the form of the nonlinear relation between the second Piola–Kirchhoff stress tensor and the Green strain tensor for anisotropic materials with secondorder nonlinearities accounted for. The effective characteristics of test models are computed by means of a CAE Fidesys program module based on the proposed algorithm. The effective material properties as functions of porosity are examined. The finite element mesh that contained more than a million of elements was used while performing stress analysis of a specimen. To reduce computing time, assembly and solution of the global equation system was done in parallel using CUDA technology. The computations were carried out on NVIDIA Tesla C2050 graphics processors. Our results show that accounting for nonlinear effects is essential for correct estimation of effective properties of porous materials. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction The problem of estimating the effective properties of porous materials [1] is important for the analysis of specimen properties during geological exploration of oil fields. Such material properties are used for computational modeling of physical processes in the course of field development and operation. At present the existing mathematical models do not allow to calculate the effective properties of porous materials with the accuracy suitable for industrial use. The effective material properties are usually determined experimentally, by testing of multiple specimens. Nevertheless, an efficient development of an oil field is impossible without complicated numerical modeling that requires significant computer power. It is necessary to perform multiple computer solutions for the determination of effective averaged properties of porous samples. Computational procedure for estimation of averaged properties of porous materials is suitable for effective parallelization using large number of computational nodes, such as in cluster-type systems. The use of hybrid supercomputers (MPI, CUDA, Open MPbased) allows enormous increase of the speed of modeling
⇑ Corresponding author. E-mail address: [email protected] (V.A. Levin). http://dx.doi.org/10.1016/j.advengsoft.2015.04.007 0965-9978/Ó 2015 Elsevier Ltd. All rights reserved.
effective prototype properties thus decreasing time and cost of laboratory tests. 2. Problem statement Let us consider the representative volume, filled with porous material. We assume that an effective (averaged) material should satisfy to the following condition: average volume stresses in real and effective materials are equal under the same displacements of faces that bound the representative volume [2–4]. Here we describe a method for creating effective constitutive relationships for the porous material with the use of definition given above. We solve certain sequences of boundary value problems of nonlinear elasticity for the representative volume V 0 in its initial state (before deformation) [6–8]:
rr¼0
ð1Þ
with boundary conditions
ujC0 ¼ R ðF e IÞ
ð2Þ
It is shown [9] that the superposition of a rigid motion on the deformation of a porous body does not change the effective constitutive equations.
81
A.V. Vershinin et al. / Advances in Engineering Software 86 (2015) 80–84
The mechanical properties of a matrix material are described by Hooke’s law or by constitutive relations of Murnaghan [5]: 0
R ¼ kðE : IÞI þ 2GE þ 3C 3 ðE : IÞ2 I þ C 4 ðE2 : IÞI þ 2C 4 ðE : IÞE þ 3C 5 E2
3. Solution method
ð3Þ 0
Here r is the true total state stress tensor (Cauchy stress tensor); R is the second Piola–Kirchhoff stress tensor; E is the Green strain tensor; I is the second-order identity tensor; F is the corresponding dr ; R is the radius vector of a affinor (deformation gradient) [6], F ¼ dR particle in the initial state, r is the radius vector of a particle in the current state; r R ¼ u; u is the displacement vector from the initial state to the current state; k and G are the first-order elastic constants; C 3 ; C 4 , and C 5 are the second-order elastic constants. The index e corresponds to the effective material. The colon in equations is the sign of double tensor contraction. Each problem sequence corresponds to certain appearance of the Green strain tensor Ee of the effective material (and to certain deformation gradient F in boundary conditions). Also, different problems of one sequence differ by strain magnitude. We find stress tensor r by solving each problem of each sequence. Using r it is possible to calculate stress tensor re of effective material according to the following formula:
re ¼
1 V
Z
r dV ¼
V
1 V
Z
N rr dC
ð4Þ
C
where N is an external normal to boundary C and r is the radius vector. Last equality in (4) has been received with the help of the divergence theorem and the following observation from tensor analysis:
r ðrrÞ ¼ ðr rÞr þ rðrrÞ ¼ ðr rÞr þ r I ¼ r
ð5Þ
Knowing the deformation gradient (that was set by us) and true stress tensor, we can calculate Green strain tensor and second Piola–Kirchhoff stress of the effective material:
Ee ¼
1 e e ðF F IÞ 2
0
1
Re ¼ ðdet F e ÞðF e Þ
ð6Þ
re ðF e Þ
1
ð7Þ
We assume that the dependence of the Piola–Kirchhoff stress tensor components on a strain parameter q for each problem sequence (i.e. for each load type) can be described by a quadratic function: 0
Reij ¼ a0ij q þ a1ij q2
cases) in order to perform least squares method for determining coefficients in (8).
ð8Þ
Accordingly we will search for effective constitutive relations as 0
a nonlinear dependence of the Piola–Kirchhoff stress tensor Re on the Green strain tensor Ee :
In practice it is convenient to set the Green strain tensor and to employ it for computation of the deformation gradient using formula (6). Since the deformation gradient is a non-symmetric tensor of the second rank it cannot be uniquely determined from the symmetric Green strain tensor. That is why we set the deformation gradient as an upper triangular matrix. Then its six components are uniquely defined by six independent components of the Green strain tensor. After calculation of the deformation gradient, we apply boundary conditions (2) to the model, solve the boundary value problem of nonlinear elasticity and find stresses r. The following sequences of problems are solved: (1) (2) (3) (4) (5) (6) (7)
E11 ¼ q – uniaxial tension or compression in the X direction; E22 ¼ q – uniaxial tension or compression in the Y direction; E33 ¼ q – uniaxial tension or compression in the Z direction; E12 ¼ E21 ¼ q – shear in the XY plane; E13 ¼ E31 ¼ q – shear in the XZ plane; E23 ¼ E32 ¼ q – shear in the YZ plane; E11 ¼ E22 ¼ q – biaxial tension or compression in the XY plane;
(8) E011 ¼ E033 ¼ q – biaxial tension or compression in the XZ plane; (9) E011 ¼ E022 ¼ E033 ¼ q – uniform tension or compression; (10) E11 ¼ q; E12 ¼ E21 ¼ q – superposition of uniaxial tension or compression in the X direction and shear in the XY plane; (11) E22 ¼ q; E12 ¼ E21 ¼ q – superposition of uniaxial tension or compression in the Y direction and shear in the XY plane; (12) E33 ¼ q; E12 ¼ E21 ¼ q – superposition of uniaxial tension or compression in the Z direction and shear in the XY plane; (13) E11 ¼ q; E13 ¼ E31 ¼ q – superposition of uniaxial tension or compression in the X direction and shear in the XZ plane; (14) E22 ¼ q; E13 ¼ E31 ¼ q – superposition of uniaxial tension or compression in the Y direction and shear in the XZ plane; (15) E33 ¼ q; E13 ¼ E31 ¼ q – superposition of uniaxial tension or compression in the Z direction and shear in the XZ plane; (16) E11 ¼ q; E23 ¼ E32 ¼ q – superposition of uniaxial tension or compression in the X direction and shear in the YZ plane; (17) E22 ¼ q; E23 ¼ E32 ¼ q – superposition of uniaxial tension or compression in the Y direction and shear in the YZ plane; (18) E33 ¼ q; E23 ¼ E32 ¼ q – superposition of uniaxial tension or compression in the Z direction and shear in the YZ plane; (19) E12 ¼ E21 ¼ q; E13 ¼ E31 ¼ q – superposition of shears in the XY and XZ planes; (20) E12 ¼ E21 ¼ q; E23 ¼ E32 ¼ q – superposition of shears in the XY and YZ planes; (21) E13 ¼ E31 ¼ q; E23 ¼ E32 ¼ q – superposition of shears in the XZ and YZ planes;
0
Reij ¼ C 0ijkl Eekl þ C 1ijklmn Eekl Eemn
ð9Þ
The second term in the right-hand side is introduced in order to take into account nonlinear effects in the effective constitutive equations. The effective elastic moduli C 0ijkl and C 1ijklmn of the porous material depend on porosity. It is worth noting that finding the effective constitutive relations of the porous material in the nonlinear case above requires solving 21 sequences of cases (equal to the number of independent elastic modules in (9)). Each sequence must contain no less than three problems (i.e. for each load type – no less than three load
here q is a strain parameter. The effective stress tensor re is calculated with the help of aver0
aging (4). Using Eq. (7), the Piola–Kirchhoff stress tensor Re is determined. Then the dependence of the Piola–Kirchhoff stress tensor on the strain parameter q is generated for each problem sequence. These dependences are approximated by expressions (8). Coefficients a0ij and a1ij are found using the least squares method. The coefficients C 0ijkl from (9) are determined through the coefficients a0ij , and coefficients C 1ijklmn – through a1ij that are calculated
82
A.V. Vershinin et al. / Advances in Engineering Software 86 (2015) 80–84
for the corresponding problem sequence by solving a linear algebraic system composed of 21 equations. Taking into account the symmetry of the Green strain tensor, we assume that C 0ijkl ¼ C 0ijlk . Then it is sufficient to find the following coefficients C 0ijkl : C 0ij11
C 0ij22
C 0ij12
C 0ij23
C 0ij33
C 0ij13 These six coefficients are calculated from the results of first six problem sequences by solving a system of six linear algebraic equations. Taking into account the Green strain tensor symmetry and the appearance of effective constitutive relations we assume that C 1ijklmn ¼ C 1ijlkmn ; C 1ijklmn ¼ C 1ijklnm ; C 1ijklmn ¼ C 1ijmnkl . Then it is sufficient to find the following coefficients C 1ijklmn : C 1ij1111
C 1ij1212
C 1ij1313
C 1ij2222
C 1ij2323
C 1ij1112
C 1ij1213
C 1ij1322
C 1ij2223
C 1ij2333
C 1ij1113 C 1ij1122 C 1ij1123 C 1ij1133
C 1ij1222 C 1ij1223 C 1ij1233
C 1ij1323 C 1ij1333
C 1ij2233
C 1ij3333
These 21 coefficients are determined by solving a system of 21 linear algebraic equations. The additional averaging of effective modules along coordinate directions is performed. This averaging permits one to write the effective constitutive equations for porous materials in the form of Murnaghan model. 4. Numerical results The proposed algorithm was implemented as a program module (plug-in) in CAE Fidesys [10] with parallelization on the basis of CUDA technology [11] thus allowing the use of the computational power of graphics processing unit (GPU) as massively parallel equipment. A program package CAE Fidesys is dedicated to solving static and dynamic stress analysis problems for elastic and elastic–plastic solids under finite strains with the use of the finite element method (FEM), spectral element method (SEM) and the discontinuous Galerkin method (DG). The package has been developed in such a way that it can be executed on an ordinary desktop or notebook PC computer. If a PC computer has a powerful video card, computations can be made in parallel with CUDA technology giving considerable speedup (up to 30 times). In addition, the package is suitable for supercomputer usage with MPI-based parallelization. A uniqueness of CAE Fidesys package is related to a possibility of solving problems in which new boundary surfaces (cavities) are originated under loading or material properties are changed in some part of the solid that leads to the redistribution of finite strains in the computational domain [7,12]. The package allows direct modeling of the redistribution of finite strains. Unlike schemes with ‘‘element killing’’ or change of material properties inside elements (when boundary conditions are not taken into account precisely), this software uses an approach based on the exact expressions for the full analytical problem statement including a forced removal of the solid part. An approach based on the direct replacement of the system of nonlinear partial differential
equations by the system of nonlinear algebraic ones by means of FEM or SEM and their further solution is used in the package. CAE Fidesys includes some enhancements of computational schemes such as: 1. Enhancement of models, computational methods and algorithms, in particular: 1.1. Better precision for modeling of nonlinear effects arising under finite strains and their redistribution. 1.2. Possibility to take into account a forced change of the solid shape (structural element shape) under large deformations with the origination of new surfaces in a loaded solid. 1.3. Modeling a change of material properties (including phase transitions in solids) under large deformations including the analysis for nano-size templates. 1.4. Using modern variants of constitutive equations. 1.5. Solving different types of coupled problems including the case of finite strains. 1.6. Application of the newest computational schemes to solving problems of mathematical physics, for example, a fully explicit spectral element method and discontinuous Galerkin method. 2. Use of modern computational schemes, software and hardware: 2.1. Use of modern software libraries and tools that allow drastically simplify program development and speed-up computational processes. 2.2. Parallelization of computational processes with the use of shared memory systems, computer clusters and massively-parallel computers. 2.3. Application of newest technologies to the use of additional PC resources, for example, using CUDA technology for parallel computations of the GPU. 2.4. Cloud computing support by means of SAAS (software as a service) Fidesys-online [13] – a web-interface to CAE Fidesys for performing remote computations in a cloud. While performing numerical simulation in CAE Fidesys both assembly of the global finite element stiffness matrix [14,15] and equation solution [16] are performed in parallel. The assembly of the global stiffness matrix requires stiffness matrix calculation for all finite elements. In our program, each thread calculates its local element matrix, and saves it to the global stiffness matrix. The conjugate gradient method is selected for solving the global equation system. Its algorithm includes just computing sparse matrix–vector product and scalar products for vectors. These operations are relatively easy for parallelization. Profiling of a computational process demonstrated that in case of GPU the most time consuming part is sparse system solution (47% of a total time), while global stiffness matrix assembly takes about 23%. That is why it was decided to accelerate these parts by their parallelization on the GPU. The parallel version of the program based on CUDA technology demonstrated speedup of 6.8 compared to the sequential program when using a workstation with two NVIDIA Tesla C2050 graphics processors. Fig. 1 illustrates an example of a finite element mesh for the porous medium built in CAE Fidesys. Finite element mesh contains 2,323,969 tetrahedral elements and 401,946 nodes. As was mentioned in the description of the algorithm for estimation of effective properties, the computational model (the representative volume of porous material) is subjected to 21 load types. Fig. 2 shows analysis results for one of load types – displacement distribution in the computational model. For each load type we perform 3–4 solutions with different load levels. Dependences of the second Piola–Kirchhoff stress components on load magnitude are determined for each load type (Eq. (8)) using
83
A.V. Vershinin et al. / Advances in Engineering Software 86 (2015) 80–84
arrangements is performed with regard for their probabilities, which depend on porosity. Let P denote the probability of existence of a pore, whose center is at a certain node. It is assumed that P does not depend on the existence of pores whose centers are in other nodes. Taking into account these assumptions, one can obtain an expression for porosity ratio p:
p¼P
V pore h
3
¼P
4pq3 3h
3
ð10Þ
;
where q is the pore radius, and V pore ¼ ð4=3Þpq3 is the pore volume. From (10) one can express the probability P in terms of porosity: 3
P¼p
Fig. 1. Example of finite element mesh in CAE Fidesys.
3h : 4pq3
ð11Þ
Let a be an arrangement of pores, and there are k pores in this arrangement. With regard for the assumptions made, one can obtain the probability P a of this arrangement:
Pa ¼ Pk ð1 PÞMk : the least squares method. Then the coefficients of the effective constitutive relations (9) are computed. For the modeling of the porous medium, the ensemble averaging is used. The scheme of ensemble averaging is similar to the scheme that was proposed in [17] for 2D case. This scheme approximately simulates a statistically uniform and isotropic distribution of pores. We suppose that all the pores assume a spherical shape in the undeformed state and are uniform in size. Let the representative volume V 0 be a cube of side a, whose sides are parallel to the coordinate axes. The uniform mesh with nodes (xi ; yj ; zk Þ is constructed in V 0 :
xi ¼ ði þ 1=2Þh;
yj ¼ ðj þ 1=2Þh;
zk ¼ ðk þ 1=2Þh;
ð12Þ C 0ijkl ; C 1ijklmn
Let C ðaÞ be one of the moduli obtained for the arrangement a, and hC i be the correspondent modulus obtained by averaging by the arrangements. Then one can write
hC i ¼
X C ðaÞ Pa ;
ð13Þ
a
where summation is made over all arrangements. Using Eqs. (11)–(13), one can express the effective elastic moduli in terms of porosity. Note that the condition q < h=2 is necessary and sufficient to avoid intersections of boundaries of pores with one another and with the boundary of the representative volume. It follows from (10) that the limiting value of porosity for this 3
where i ¼ 0; . . . N 1; j ¼ 0; . . . N 1; k ¼ 0; . . . N 1; h ¼ a=N; N is a given integer. So, the mesh consists of M ¼ N 3 nodes. It is assumed that centers of pores are located exactly on the nodes of the mesh. There are 2M arrangements of pores satisfying this assumption. For each arrangement, the effective moduli C 0ijkl ; C 1ijklmn are computed using the scheme presented above. Then the averaging over all these
scheme is pmax ¼ 4pq3 =ð3h Þ. After the computation of moduli hC i0ijkl ; hC i1ijklmn , their averaging by all possible orientations of coordinate axes is performed [17]. This averaging permits one to represent the effective constitutive equations in a form (3) proposed by Murnaghan. The dependences of effective modules of the first and the second order on the porosity ratio are shown in Figs. 3 and 4. The computations are performed for N ¼ 2; q=h ¼ 0:42. The effective
Fig. 2. Displacement distribution in computational model of porous material sample.
84
A.V. Vershinin et al. / Advances in Engineering Software 86 (2015) 80–84
of a stiffness matrix and solution of a sparse equation system. The dependence of effective material characteristics on the porosity is examined. It was found that accounting for nonlinear effects is essential for correct determination of effective properties for porous materials undergoing large strains. Acknowledgments
Fig. 3. Dependence of the linear effective modules on the porosity ratio p.
The research for this article was financially supported by Russian Ministry of Education and Science (State contract 14.579.21.0007, project ID RFMEFI57914X0007). Investigations were performed in the FIDESYS company, a subsidiary of the Ministry of Education and Science of the Russian Federation. References
Fig. 4. Dependence of the second order effective modules on the porosity ratio p.
modules are related to the shear modulus of the solid material. The results show that the dependence of nonlinear modules on the porosity ratio (Fig. 4) is not always monotone. 5. Conclusion A computational procedure has been proposed that allows numerical estimation of effective mechanical properties for porous materials under finite strains. The effective properties are sought in the form of nonlinear dependence between the Piola–Kirchhoff stress tensor and the Green strain tensor for anisotropic materials. An algorithm for the estimation of effective properties for porous materials has been implemented in CAE Fidesys with the use of CUDA technology in order to perform parallelization of assembly
[1] Zingerman KM, Levin VA, Yakovlev MYa, Prokopenko AV, Terpyakov AA. Computation of effective elastic characteristics of porous and composite materials using the FIDESYS CAE-system. In: Proceedings of the 10th word congress on computational mechanics, Sao Paulo, Brazil; 2012. p. 290. [2] Levin VA, Zingermann KM. Effective constitutive equations for porous elastic materials at finite strains and superimposed finite strains. Trans ASME: J Appl Mech 2003;70:982–5. [3] Levin VA, Zingerman KM. Interaction and microfracturing pattern for successive origination (introduction) of pores in elastic bodies: finite deformation. Trans ASME: J Appl Mech 1998;65:431–5. [4] Tsukrov I, Kachanov M. Stress concentrations and microfracturing patterns in a brittle-elastic solid with interacting pores of diverse shapes. Int J Solids Struct 1997;34:2887–904. [5] Murnaghan FD. Finite deformation of an elastic solid. New York: Willey; 1951. [6] Lurie AI. Nonlinear theory of elasticity. North-Holland: Amsterdam; 1990. [7] Levin VA. Theory of repeated superposition of large deformations. Elastic and viscoelastic bodies. Int J Solids Struct 1998;35:2585–600. [8] Eshelby JD. The continual theory of dislocations. IL: Moscow; 1963. [9] Levin VA, Zingerman KM. On the construction of effective constitutive relations for porous elastic materials subjected to finite deformations including the case of their superposition. Dokl Phys 2002;47(2):136–40. [10] http://www.cae-fidesys.com. [11] Kirk DB, Hwu WW. Programming massively parallel processors: a hands-on approach. Burligton: Morgan Kaufmann; 2010. [12] Levin VA, Vershinin AV. Non-stationary plane problem of the successive origination of stress concentrators in a loaded body. Finite deformations and their superposition. Commun Numer Methods Eng 2008;24:2229–39. [13] http://online.cae-fidesys.com. [14] Zienkiewicz IC, Taylor RL. The finite element method. The basis, vol. 1. London: McGraw-Hill; 1971. [15] Zienkiewicz OC, Taylor RL. The finite element method. Solid mechanics, vol. 2. London: McGraw-Hill; 1971. [16] Barrett R et al. Templates for the solution of linear systems: building blocks for iterative methods. Philadelphia: SIAM; 1994. [17] Levin VA, Lokhin VV, Zingerman KM. Effective elastic properties of porous materials with randomly dispersed pores: finite deformation. Trans ASME: J Appl Mech 2000;67:667–70.
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
Software for estimation of second order effective material properties of porous samples with geometrical and physical nonlinearity accounted for A.V. Vershinin a, V.A. Levin a,⇑, K.M. Zingerman b, A.M. Sboychakov c, M.Ya. Yakovlev c a b c
Lomonosov Moscow State University, Moscow, Russian Federation National Research Nuclear University MEPhI, Moscow, Russian Federation Fidesys LLC, Moscow, Russian Federation
a r t i c l e
i n f o
Article history: Received 29 August 2014 Received in revised form 10 March 2015 Accepted 11 April 2015
Keywords: Effective properties Porous material Finite strain Physical nonlinearity Finite-element method Parallel computations
a b s t r a c t A method and an algorithm for numerical estimation of effective mechanical properties of porous materials are presented. The effective properties are sought in the form of the nonlinear relation between the second Piola–Kirchhoff stress tensor and the Green strain tensor for anisotropic materials with secondorder nonlinearities accounted for. The effective characteristics of test models are computed by means of a CAE Fidesys program module based on the proposed algorithm. The effective material properties as functions of porosity are examined. The finite element mesh that contained more than a million of elements was used while performing stress analysis of a specimen. To reduce computing time, assembly and solution of the global equation system was done in parallel using CUDA technology. The computations were carried out on NVIDIA Tesla C2050 graphics processors. Our results show that accounting for nonlinear effects is essential for correct estimation of effective properties of porous materials. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction The problem of estimating the effective properties of porous materials [1] is important for the analysis of specimen properties during geological exploration of oil fields. Such material properties are used for computational modeling of physical processes in the course of field development and operation. At present the existing mathematical models do not allow to calculate the effective properties of porous materials with the accuracy suitable for industrial use. The effective material properties are usually determined experimentally, by testing of multiple specimens. Nevertheless, an efficient development of an oil field is impossible without complicated numerical modeling that requires significant computer power. It is necessary to perform multiple computer solutions for the determination of effective averaged properties of porous samples. Computational procedure for estimation of averaged properties of porous materials is suitable for effective parallelization using large number of computational nodes, such as in cluster-type systems. The use of hybrid supercomputers (MPI, CUDA, Open MPbased) allows enormous increase of the speed of modeling
⇑ Corresponding author. E-mail address: [email protected] (V.A. Levin). http://dx.doi.org/10.1016/j.advengsoft.2015.04.007 0965-9978/Ó 2015 Elsevier Ltd. All rights reserved.
effective prototype properties thus decreasing time and cost of laboratory tests. 2. Problem statement Let us consider the representative volume, filled with porous material. We assume that an effective (averaged) material should satisfy to the following condition: average volume stresses in real and effective materials are equal under the same displacements of faces that bound the representative volume [2–4]. Here we describe a method for creating effective constitutive relationships for the porous material with the use of definition given above. We solve certain sequences of boundary value problems of nonlinear elasticity for the representative volume V 0 in its initial state (before deformation) [6–8]:
rr¼0
ð1Þ
with boundary conditions
ujC0 ¼ R ðF e IÞ
ð2Þ
It is shown [9] that the superposition of a rigid motion on the deformation of a porous body does not change the effective constitutive equations.
81
A.V. Vershinin et al. / Advances in Engineering Software 86 (2015) 80–84
The mechanical properties of a matrix material are described by Hooke’s law or by constitutive relations of Murnaghan [5]: 0
R ¼ kðE : IÞI þ 2GE þ 3C 3 ðE : IÞ2 I þ C 4 ðE2 : IÞI þ 2C 4 ðE : IÞE þ 3C 5 E2
3. Solution method
ð3Þ 0
Here r is the true total state stress tensor (Cauchy stress tensor); R is the second Piola–Kirchhoff stress tensor; E is the Green strain tensor; I is the second-order identity tensor; F is the corresponding dr ; R is the radius vector of a affinor (deformation gradient) [6], F ¼ dR particle in the initial state, r is the radius vector of a particle in the current state; r R ¼ u; u is the displacement vector from the initial state to the current state; k and G are the first-order elastic constants; C 3 ; C 4 , and C 5 are the second-order elastic constants. The index e corresponds to the effective material. The colon in equations is the sign of double tensor contraction. Each problem sequence corresponds to certain appearance of the Green strain tensor Ee of the effective material (and to certain deformation gradient F in boundary conditions). Also, different problems of one sequence differ by strain magnitude. We find stress tensor r by solving each problem of each sequence. Using r it is possible to calculate stress tensor re of effective material according to the following formula:
re ¼
1 V
Z
r dV ¼
V
1 V
Z
N rr dC
ð4Þ
C
where N is an external normal to boundary C and r is the radius vector. Last equality in (4) has been received with the help of the divergence theorem and the following observation from tensor analysis:
r ðrrÞ ¼ ðr rÞr þ rðrrÞ ¼ ðr rÞr þ r I ¼ r
ð5Þ
Knowing the deformation gradient (that was set by us) and true stress tensor, we can calculate Green strain tensor and second Piola–Kirchhoff stress of the effective material:
Ee ¼
1 e e ðF F IÞ 2
0
1
Re ¼ ðdet F e ÞðF e Þ
ð6Þ
re ðF e Þ
1
ð7Þ
We assume that the dependence of the Piola–Kirchhoff stress tensor components on a strain parameter q for each problem sequence (i.e. for each load type) can be described by a quadratic function: 0
Reij ¼ a0ij q þ a1ij q2
cases) in order to perform least squares method for determining coefficients in (8).
ð8Þ
Accordingly we will search for effective constitutive relations as 0
a nonlinear dependence of the Piola–Kirchhoff stress tensor Re on the Green strain tensor Ee :
In practice it is convenient to set the Green strain tensor and to employ it for computation of the deformation gradient using formula (6). Since the deformation gradient is a non-symmetric tensor of the second rank it cannot be uniquely determined from the symmetric Green strain tensor. That is why we set the deformation gradient as an upper triangular matrix. Then its six components are uniquely defined by six independent components of the Green strain tensor. After calculation of the deformation gradient, we apply boundary conditions (2) to the model, solve the boundary value problem of nonlinear elasticity and find stresses r. The following sequences of problems are solved: (1) (2) (3) (4) (5) (6) (7)
E11 ¼ q – uniaxial tension or compression in the X direction; E22 ¼ q – uniaxial tension or compression in the Y direction; E33 ¼ q – uniaxial tension or compression in the Z direction; E12 ¼ E21 ¼ q – shear in the XY plane; E13 ¼ E31 ¼ q – shear in the XZ plane; E23 ¼ E32 ¼ q – shear in the YZ plane; E11 ¼ E22 ¼ q – biaxial tension or compression in the XY plane;
(8) E011 ¼ E033 ¼ q – biaxial tension or compression in the XZ plane; (9) E011 ¼ E022 ¼ E033 ¼ q – uniform tension or compression; (10) E11 ¼ q; E12 ¼ E21 ¼ q – superposition of uniaxial tension or compression in the X direction and shear in the XY plane; (11) E22 ¼ q; E12 ¼ E21 ¼ q – superposition of uniaxial tension or compression in the Y direction and shear in the XY plane; (12) E33 ¼ q; E12 ¼ E21 ¼ q – superposition of uniaxial tension or compression in the Z direction and shear in the XY plane; (13) E11 ¼ q; E13 ¼ E31 ¼ q – superposition of uniaxial tension or compression in the X direction and shear in the XZ plane; (14) E22 ¼ q; E13 ¼ E31 ¼ q – superposition of uniaxial tension or compression in the Y direction and shear in the XZ plane; (15) E33 ¼ q; E13 ¼ E31 ¼ q – superposition of uniaxial tension or compression in the Z direction and shear in the XZ plane; (16) E11 ¼ q; E23 ¼ E32 ¼ q – superposition of uniaxial tension or compression in the X direction and shear in the YZ plane; (17) E22 ¼ q; E23 ¼ E32 ¼ q – superposition of uniaxial tension or compression in the Y direction and shear in the YZ plane; (18) E33 ¼ q; E23 ¼ E32 ¼ q – superposition of uniaxial tension or compression in the Z direction and shear in the YZ plane; (19) E12 ¼ E21 ¼ q; E13 ¼ E31 ¼ q – superposition of shears in the XY and XZ planes; (20) E12 ¼ E21 ¼ q; E23 ¼ E32 ¼ q – superposition of shears in the XY and YZ planes; (21) E13 ¼ E31 ¼ q; E23 ¼ E32 ¼ q – superposition of shears in the XZ and YZ planes;
0
Reij ¼ C 0ijkl Eekl þ C 1ijklmn Eekl Eemn
ð9Þ
The second term in the right-hand side is introduced in order to take into account nonlinear effects in the effective constitutive equations. The effective elastic moduli C 0ijkl and C 1ijklmn of the porous material depend on porosity. It is worth noting that finding the effective constitutive relations of the porous material in the nonlinear case above requires solving 21 sequences of cases (equal to the number of independent elastic modules in (9)). Each sequence must contain no less than three problems (i.e. for each load type – no less than three load
here q is a strain parameter. The effective stress tensor re is calculated with the help of aver0
aging (4). Using Eq. (7), the Piola–Kirchhoff stress tensor Re is determined. Then the dependence of the Piola–Kirchhoff stress tensor on the strain parameter q is generated for each problem sequence. These dependences are approximated by expressions (8). Coefficients a0ij and a1ij are found using the least squares method. The coefficients C 0ijkl from (9) are determined through the coefficients a0ij , and coefficients C 1ijklmn – through a1ij that are calculated
82
A.V. Vershinin et al. / Advances in Engineering Software 86 (2015) 80–84
for the corresponding problem sequence by solving a linear algebraic system composed of 21 equations. Taking into account the symmetry of the Green strain tensor, we assume that C 0ijkl ¼ C 0ijlk . Then it is sufficient to find the following coefficients C 0ijkl : C 0ij11
C 0ij22
C 0ij12
C 0ij23
C 0ij33
C 0ij13 These six coefficients are calculated from the results of first six problem sequences by solving a system of six linear algebraic equations. Taking into account the Green strain tensor symmetry and the appearance of effective constitutive relations we assume that C 1ijklmn ¼ C 1ijlkmn ; C 1ijklmn ¼ C 1ijklnm ; C 1ijklmn ¼ C 1ijmnkl . Then it is sufficient to find the following coefficients C 1ijklmn : C 1ij1111
C 1ij1212
C 1ij1313
C 1ij2222
C 1ij2323
C 1ij1112
C 1ij1213
C 1ij1322
C 1ij2223
C 1ij2333
C 1ij1113 C 1ij1122 C 1ij1123 C 1ij1133
C 1ij1222 C 1ij1223 C 1ij1233
C 1ij1323 C 1ij1333
C 1ij2233
C 1ij3333
These 21 coefficients are determined by solving a system of 21 linear algebraic equations. The additional averaging of effective modules along coordinate directions is performed. This averaging permits one to write the effective constitutive equations for porous materials in the form of Murnaghan model. 4. Numerical results The proposed algorithm was implemented as a program module (plug-in) in CAE Fidesys [10] with parallelization on the basis of CUDA technology [11] thus allowing the use of the computational power of graphics processing unit (GPU) as massively parallel equipment. A program package CAE Fidesys is dedicated to solving static and dynamic stress analysis problems for elastic and elastic–plastic solids under finite strains with the use of the finite element method (FEM), spectral element method (SEM) and the discontinuous Galerkin method (DG). The package has been developed in such a way that it can be executed on an ordinary desktop or notebook PC computer. If a PC computer has a powerful video card, computations can be made in parallel with CUDA technology giving considerable speedup (up to 30 times). In addition, the package is suitable for supercomputer usage with MPI-based parallelization. A uniqueness of CAE Fidesys package is related to a possibility of solving problems in which new boundary surfaces (cavities) are originated under loading or material properties are changed in some part of the solid that leads to the redistribution of finite strains in the computational domain [7,12]. The package allows direct modeling of the redistribution of finite strains. Unlike schemes with ‘‘element killing’’ or change of material properties inside elements (when boundary conditions are not taken into account precisely), this software uses an approach based on the exact expressions for the full analytical problem statement including a forced removal of the solid part. An approach based on the direct replacement of the system of nonlinear partial differential
equations by the system of nonlinear algebraic ones by means of FEM or SEM and their further solution is used in the package. CAE Fidesys includes some enhancements of computational schemes such as: 1. Enhancement of models, computational methods and algorithms, in particular: 1.1. Better precision for modeling of nonlinear effects arising under finite strains and their redistribution. 1.2. Possibility to take into account a forced change of the solid shape (structural element shape) under large deformations with the origination of new surfaces in a loaded solid. 1.3. Modeling a change of material properties (including phase transitions in solids) under large deformations including the analysis for nano-size templates. 1.4. Using modern variants of constitutive equations. 1.5. Solving different types of coupled problems including the case of finite strains. 1.6. Application of the newest computational schemes to solving problems of mathematical physics, for example, a fully explicit spectral element method and discontinuous Galerkin method. 2. Use of modern computational schemes, software and hardware: 2.1. Use of modern software libraries and tools that allow drastically simplify program development and speed-up computational processes. 2.2. Parallelization of computational processes with the use of shared memory systems, computer clusters and massively-parallel computers. 2.3. Application of newest technologies to the use of additional PC resources, for example, using CUDA technology for parallel computations of the GPU. 2.4. Cloud computing support by means of SAAS (software as a service) Fidesys-online [13] – a web-interface to CAE Fidesys for performing remote computations in a cloud. While performing numerical simulation in CAE Fidesys both assembly of the global finite element stiffness matrix [14,15] and equation solution [16] are performed in parallel. The assembly of the global stiffness matrix requires stiffness matrix calculation for all finite elements. In our program, each thread calculates its local element matrix, and saves it to the global stiffness matrix. The conjugate gradient method is selected for solving the global equation system. Its algorithm includes just computing sparse matrix–vector product and scalar products for vectors. These operations are relatively easy for parallelization. Profiling of a computational process demonstrated that in case of GPU the most time consuming part is sparse system solution (47% of a total time), while global stiffness matrix assembly takes about 23%. That is why it was decided to accelerate these parts by their parallelization on the GPU. The parallel version of the program based on CUDA technology demonstrated speedup of 6.8 compared to the sequential program when using a workstation with two NVIDIA Tesla C2050 graphics processors. Fig. 1 illustrates an example of a finite element mesh for the porous medium built in CAE Fidesys. Finite element mesh contains 2,323,969 tetrahedral elements and 401,946 nodes. As was mentioned in the description of the algorithm for estimation of effective properties, the computational model (the representative volume of porous material) is subjected to 21 load types. Fig. 2 shows analysis results for one of load types – displacement distribution in the computational model. For each load type we perform 3–4 solutions with different load levels. Dependences of the second Piola–Kirchhoff stress components on load magnitude are determined for each load type (Eq. (8)) using
83
A.V. Vershinin et al. / Advances in Engineering Software 86 (2015) 80–84
arrangements is performed with regard for their probabilities, which depend on porosity. Let P denote the probability of existence of a pore, whose center is at a certain node. It is assumed that P does not depend on the existence of pores whose centers are in other nodes. Taking into account these assumptions, one can obtain an expression for porosity ratio p:
p¼P
V pore h
3
¼P
4pq3 3h
3
ð10Þ
;
where q is the pore radius, and V pore ¼ ð4=3Þpq3 is the pore volume. From (10) one can express the probability P in terms of porosity: 3
P¼p
Fig. 1. Example of finite element mesh in CAE Fidesys.
3h : 4pq3
ð11Þ
Let a be an arrangement of pores, and there are k pores in this arrangement. With regard for the assumptions made, one can obtain the probability P a of this arrangement:
Pa ¼ Pk ð1 PÞMk : the least squares method. Then the coefficients of the effective constitutive relations (9) are computed. For the modeling of the porous medium, the ensemble averaging is used. The scheme of ensemble averaging is similar to the scheme that was proposed in [17] for 2D case. This scheme approximately simulates a statistically uniform and isotropic distribution of pores. We suppose that all the pores assume a spherical shape in the undeformed state and are uniform in size. Let the representative volume V 0 be a cube of side a, whose sides are parallel to the coordinate axes. The uniform mesh with nodes (xi ; yj ; zk Þ is constructed in V 0 :
xi ¼ ði þ 1=2Þh;
yj ¼ ðj þ 1=2Þh;
zk ¼ ðk þ 1=2Þh;
ð12Þ C 0ijkl ; C 1ijklmn
Let C ðaÞ be one of the moduli obtained for the arrangement a, and hC i be the correspondent modulus obtained by averaging by the arrangements. Then one can write
hC i ¼
X C ðaÞ Pa ;
ð13Þ
a
where summation is made over all arrangements. Using Eqs. (11)–(13), one can express the effective elastic moduli in terms of porosity. Note that the condition q < h=2 is necessary and sufficient to avoid intersections of boundaries of pores with one another and with the boundary of the representative volume. It follows from (10) that the limiting value of porosity for this 3
where i ¼ 0; . . . N 1; j ¼ 0; . . . N 1; k ¼ 0; . . . N 1; h ¼ a=N; N is a given integer. So, the mesh consists of M ¼ N 3 nodes. It is assumed that centers of pores are located exactly on the nodes of the mesh. There are 2M arrangements of pores satisfying this assumption. For each arrangement, the effective moduli C 0ijkl ; C 1ijklmn are computed using the scheme presented above. Then the averaging over all these
scheme is pmax ¼ 4pq3 =ð3h Þ. After the computation of moduli hC i0ijkl ; hC i1ijklmn , their averaging by all possible orientations of coordinate axes is performed [17]. This averaging permits one to represent the effective constitutive equations in a form (3) proposed by Murnaghan. The dependences of effective modules of the first and the second order on the porosity ratio are shown in Figs. 3 and 4. The computations are performed for N ¼ 2; q=h ¼ 0:42. The effective
Fig. 2. Displacement distribution in computational model of porous material sample.
84
A.V. Vershinin et al. / Advances in Engineering Software 86 (2015) 80–84
of a stiffness matrix and solution of a sparse equation system. The dependence of effective material characteristics on the porosity is examined. It was found that accounting for nonlinear effects is essential for correct determination of effective properties for porous materials undergoing large strains. Acknowledgments
Fig. 3. Dependence of the linear effective modules on the porosity ratio p.
The research for this article was financially supported by Russian Ministry of Education and Science (State contract 14.579.21.0007, project ID RFMEFI57914X0007). Investigations were performed in the FIDESYS company, a subsidiary of the Ministry of Education and Science of the Russian Federation. References
Fig. 4. Dependence of the second order effective modules on the porosity ratio p.
modules are related to the shear modulus of the solid material. The results show that the dependence of nonlinear modules on the porosity ratio (Fig. 4) is not always monotone. 5. Conclusion A computational procedure has been proposed that allows numerical estimation of effective mechanical properties for porous materials under finite strains. The effective properties are sought in the form of nonlinear dependence between the Piola–Kirchhoff stress tensor and the Green strain tensor for anisotropic materials. An algorithm for the estimation of effective properties for porous materials has been implemented in CAE Fidesys with the use of CUDA technology in order to perform parallelization of assembly
[1] Zingerman KM, Levin VA, Yakovlev MYa, Prokopenko AV, Terpyakov AA. Computation of effective elastic characteristics of porous and composite materials using the FIDESYS CAE-system. In: Proceedings of the 10th word congress on computational mechanics, Sao Paulo, Brazil; 2012. p. 290. [2] Levin VA, Zingermann KM. Effective constitutive equations for porous elastic materials at finite strains and superimposed finite strains. Trans ASME: J Appl Mech 2003;70:982–5. [3] Levin VA, Zingerman KM. Interaction and microfracturing pattern for successive origination (introduction) of pores in elastic bodies: finite deformation. Trans ASME: J Appl Mech 1998;65:431–5. [4] Tsukrov I, Kachanov M. Stress concentrations and microfracturing patterns in a brittle-elastic solid with interacting pores of diverse shapes. Int J Solids Struct 1997;34:2887–904. [5] Murnaghan FD. Finite deformation of an elastic solid. New York: Willey; 1951. [6] Lurie AI. Nonlinear theory of elasticity. North-Holland: Amsterdam; 1990. [7] Levin VA. Theory of repeated superposition of large deformations. Elastic and viscoelastic bodies. Int J Solids Struct 1998;35:2585–600. [8] Eshelby JD. The continual theory of dislocations. IL: Moscow; 1963. [9] Levin VA, Zingerman KM. On the construction of effective constitutive relations for porous elastic materials subjected to finite deformations including the case of their superposition. Dokl Phys 2002;47(2):136–40. [10] http://www.cae-fidesys.com. [11] Kirk DB, Hwu WW. Programming massively parallel processors: a hands-on approach. Burligton: Morgan Kaufmann; 2010. [12] Levin VA, Vershinin AV. Non-stationary plane problem of the successive origination of stress concentrators in a loaded body. Finite deformations and their superposition. Commun Numer Methods Eng 2008;24:2229–39. [13] http://online.cae-fidesys.com. [14] Zienkiewicz IC, Taylor RL. The finite element method. The basis, vol. 1. London: McGraw-Hill; 1971. [15] Zienkiewicz OC, Taylor RL. The finite element method. Solid mechanics, vol. 2. London: McGraw-Hill; 1971. [16] Barrett R et al. Templates for the solution of linear systems: building blocks for iterative methods. Philadelphia: SIAM; 1994. [17] Levin VA, Lokhin VV, Zingerman KM. Effective elastic properties of porous materials with randomly dispersed pores: finite deformation. Trans ASME: J Appl Mech 2000;67:667–70.