A simple Python code for computing effective properties of
2D and 3D representative volume element under periodic
boundary conditions
Fan Ye, Hu Wang*
State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,
Hunan University, Changsha 410082, P.R. China

Abstract Multiscale optimization is an attractive research field recently. For the most
of optimization tools, design parameters should be updated during a close loop.
Therefore, a simple Python code is programmed to obtain effective properties of
Representative Volume Element (RVE) under Periodic Boundary Conditions (PBCs).
It can compute the mechanical properties of a composite with a periodic structure, in
two or three dimensions. The computation method is based on the Asymptotic
Homogenization Theory (AHT). With simple modifications, the basic Python code
may be extended to the computation of the effective properties of more complex
microstructure. Moreover, the code provides a convenient platform upon the
optimization for the material and geometric composite design. The user may
experiment with various algorithms and tackle a wide range of problems. To verify the
effectiveness and reliability of the code, a three-dimensional case is employed to
illuminate the code. Finally numerical results obtained by the code agree well with the
available theoretical and experimental results.

Keywords

Representative Volume Element

Periodic Boundary Conditions

Asymptotic Homogenization Theory Multiscale

1. Introduction
Prediction of the mechanical properties of the composites has become an active
research area [1-3]. Except for experimental studies and macromechanical methods,
micromechanical methods are widely used to obtain overall properties of composites.
The methods usually contain the introduction of two scales, marcroscale and
microscale. Marcroscale is usually referred to a homogenized continuous medium,
*

Corresponding author, wanghu@hnu.edu.cn, Tel: +86-0731-88821417, Fax: +86-0731-88822051

and microscale is usually related to a statistically RVE [4]. The RVE must be selected
such that the microstructure should be composed as copies of RVEs, and without
overlapping boundaries or gaps between boundaries. The choice of RVE is usually not
unique, but it should be big enough to present the feature of composite and as small as
possible to reduce the computational cost. Additionally, the RVE should have the
same volume fraction as the composite [5]. Micromechanical methods provide overall
behaviors of the composites from known properties of their constituents (inclusion
and matrix) through an analysis of a RVE [6-7].
The above micromechanical models can be regarded as mechanical or
engineering models. AHT is proposed to predict overall behaviors of micromechanical
models [8-11]. Moreover, the AHT has explicitly used periodic boundary conditions
in modeling of linear and nonlinear composite materials [12]. On the other hand, the
theory for predicting overall behaviors of micromechanical models is laid on a
rigorous mechanical foundation by using strain energy equivalence principles in
conjunction with Finite Element (FE) analysis [13]. Explicit unified form of boundary
conditions for RVE is developed which satisfies the periodicity conditions, and is
suitable for any combination of multiaxial loads [14].
The AHT combined with the FE method is proven to be a powerful technique
[15-16], which can consider a more complex microstructure given by several
inclusions with different shapes, different orientations and different aspect ratios, and
even random distribution of inclusions. For this theory based on FEM, there is no
restriction on geometry and material, and the AHT based FE has been well
documented for the determination of the effective material properties of the
composites.
ABAQUS is a general purpose commercial software package widely used in the
analysis of the RVEs [17, 18]. Moreover, the Python code takes advantage of the
advanced FE analysis capacities of ABAQUS software and employs the ABAQUS
Scripting Interface (ASI) to communicate with ABAQUS. ASI is an extension of the
Python language. As a standard component of the ABAQUS software, ASI provides a
convenient interface to the models and results [19]. In this study, a Python code is
developed for computing effective properties of the RVE under PBCs based on the
secondary development of ABAQUS.

2. Periodic boundary conditions
Assumed that the composite is heterogeneous with a periodic microstructure, the
RVE is subjected to boundary conditions depending on macroscale. The displacement
field for the periodic structure can be expressed as [8]

ui ( x1 , x2 , x3 )   ik xk  ui* ( x1 , x2 , x3 )

(1)

where  ik is the macro (average) strain tensor of the periodic structure and the first
term on the right side represents a linear distributed displacement field. The second
term on the right side, ui* ( x1 , x2 , x3 ) , is a periodic function from one RVE to another.
It represents a modification to the linear displacement field due to the heterogeneous
structure of the composite.
Since the periodic array of RVEs represents a continuous physical body, two
continuities must be satisfied at the boundaries of neighboring RVEs. One is that
displacements must be continuous, and neighboring RVEs cannot be separated or
overlap each other at the boundaries after deformation. The second condition implies
that the traction distributions at opposite parallel boundaries of a RVE must be the
same. In this manner, the individual RVE can be assembled as a physically continuous
body.
Obviously, the assumption of the displacement field in the form of Eq. (2) meets
the first of above requirements. Unfortunately, it cannot be directly applied to the
boundaries since the second term on the right side of Eq. (1) is generally unknown.
For any RVE, its boundary surfaces must appear in parallel pairs, the displacements
on a pair of parallel opposite boundary surfaces can be written as

uij    ik xkj   ui

(2)

uij    ik xkj   ui

(3)

where indices on the left side identify the pair of two opposite parallel boundary
surfaces of a RVE.
It should be noted that ui* ( x1 , x2 , x3 ) is the same at two parallel boundaries
(periodicity), therefore, the difference between above two equations is

uij   uij    ik ( xkj   xkj  )   ik xkj .

(4)

Since xkj are constants for each pair of the parallel boundary surfaces, with
specified  ik , the right side becomes constants. Such equations can be easily applied

to FE analysis as nodal displacement constraint [14]. Moreover, if RVE is analyzed by
using a displacement-based FE method, the application of Eq. (4) can guarantee the
uniqueness of the solution and traction continuity conditions are automatically
satisfied instead of applying Eq. (1) directly as the boundary conditions.

3. Homogenization method
In the homogenization method, the RVE is modeled as a homogeneous
orthotropic medium with certain effective properties that describe the average material
properties of the composite [11]. To describe this macroscopically homogeneous
medium, macrostress and macrostrain are derived by averaging the stress and strain
tensor over the volume of the RVE:

 ij 
 ij 

1
VRVE
1
VRVE


V


V

( x, y, z )dV

(5)

( x, y, z )dV

(6)

ij

ij

The strain energy U * stored in the heterogeneous RVE of the volume VRVE is
U* 

1
 ij ij dV
2 VRVE

(7)

Instead the strain energy corresponding to the heterogeneous RVE, the
homogeneous RVE is used to determine the corresponding homogenized modulus. In
the homogeneous RVE, the total strain energy due to deformation is given by
1
U   ij ijVRVE
2

(8)

where U is homogeneous strain energy of RVE.  ij is average strains. VRVE
is volume of periodic domain.  ij is average stresses.
As shown in [13], Eq. (7) and Eq. (8) can be converted to a surface integral by
using Gauss theorem, thus:
U * U 

1
 ij (ui  ui )n j dS j
2 S j

(9)

where S j is the jth surface and n j is the unit outward normal. ui is the ith
displacment. ui is the ith average displacment .On the surface S j :

ui  ui
Thus:

(10)

U * U  0

(11)

The average stresses and strains quantities defined in Eq. (5) and Eq. (6) ensure
equivalence in strain energy between the equivalent homogeneous material and the
original heterogeneous material.
Alternatively, the average strains can be related to boundary displacements of the
RVE by using Gauss theorem. Average strains in Eq. (6) can be converted as:

 ij 

1
VRVE



S1

(ui n j  u j ni )dS1

(12)

where S1 denotes the outer boundary of the RVE. The relationship given by Eq. (12)
makes it possible to evaluate the volume averaged strains with using the boundary
displacements, and then avoiding the volume integration [13].
The average stresses  ij is obtained by the prescribed average strains  ij . The
tensors are the ratio of average stresses and strains and computed as follow:
Cij 

 ij
 ij

(13)

where Cij is stiffness modulus corresponding to apply deformation mode.
In the case of pure shear deformation, the tensorial shear strain  ij

 ij   ij   ji  2 ij

(14)

where  ij is the engineering shear strain. In this study, the tensorial shear strains are
suited to be the macro strain [20].

4. A three-dimensional case
To verify the homogenization method under PBCs, a 3-D RVE model is
considered. The model consists of a fiber reinforcement and matrix, with a volume
fraction of 47 %. The properties of the materials are given in Table 1.
Material
Fiber
Matrix

Table 1 Material properties of fiber and matrix
E (Gpa)
379.3
68.3

v
0.1
0.3

The unidirectional RVE is assumed to be orthotropic and linearly elastic. For 3D
RVE, from Eq. (13), the material constitutive relation of this RVE can be written as

   C  

(15)

where

 is the average strain matrix.  is the average stress matrix. C  is the

stiffness matrix as shown in Eq. (16) for the orthotropic and unidirectional RVE.
 11   C11 C12
  C
 22   12 C22
 33  C13 C23
 
0
 12   0
 13   0
0
  
0
 23   0

C13

0

0

C23

0

0

C33

0

0

0

C44

0

0

0

C55

0

0

0

0   11 
0   22 
0   33 
 
0   12 
0   13 
 
C66   23 

(16)

After obtaining  ij for given  ij from the computation of a RVE, Cij can be
obtained from Eq. (15). From Eq. (8) and Eq. (15), the relation between the
homogeneous strain energy and the average strain is
U

1
T
   VRVE
2

(17)

Thus
U
1
1
 C1112  C121 2  C131 3  C22 22  C23 2 3
VRVE 2
2

(18)

1
1
1
1
 C33 32  C44 42  C55 52  C66 62
2
2
2
2

Hence, the values of the stiffness matrix are computed from Eq. (18).

5. Python implement of the three-dimensional case
This section explains the basic form of the Python codes. The code is to be called
in ABAQUS/CAE through the menu command File-> Run Scipt. A flowchart of the
code is presented in Fig.1.

Part module

Geometry and local
coordinate system

Material and Section
module

Constituent material properties
and material Orientation

Assembly and Mesh
module

Reference points (PreStrain) , the size
of the mesh, and the type of the element

Step module

The field displayed, MPC(PBCs), and Load
case

Job module

Abaqus Solver

Post-processing
module

Effective Properties

Fig. 1 The work flowchart of the Python code for the computation of RVE

5.1. Header
As a common programming practice, a header exists at the beginning of the
script to import a number of modules and functions. For the four modules, abaqus,
abaqusConstants, caeModules, and odbAccess, they are standard components of the
ABAQUS software to employ the ASI. The numpy and math modules are used to
import arithmetic functions.
from abaqus import *
from abaqusConstants import *
from caeModules import *
from odbAccess import *
from numpy import *
import math

5.2. Geometric parameters of the RVE (Create part)
The shape of the RVE is decided by the microstructure, which can be an accurate
reflection of reality or an idealisation. While it remains at the user’s discretion which
idealisation is adopted, considerations should be made so as to have the idealised
microstructure as representative of reality as possible. In this study, in the case of
uniderectionally fiber-reinforced composites with fibers distributed at the center of the
RVE, the square packing is selected to be idealistic representive element of the
composites.
This function runs the ABAQUS to creat the shape of the RVE. The geometric
parameters of the RVE are defined into Wx  Wy  Wz  W  d  1mm as shown in
Fig.2. Moreover, the fiber volume fraction of the model is designed to be 47%.
##Function of creating RVE part
def CreatePart(RVE_para,RVEModel ):
#Initial parameter
# Create RVE sketchs
RVEModel.ConstrainedSketch(name='', sheetSize='')
#Define two point to creat a rectangle
RVESketch.rectangle(point1='', point2='')
#Define the property of the part
RVEPart= RVEModel.Part(name='',dimensionality='',type='')
#Assign sketch to the part
RVEPart.BaseSolidExtrude(RVESketch, depth=)
#Create fiber sketchs
RVEt = RVEPart.MakeSketchTransform(sketchPlane='', sketchUpEdge='')
RVESketch1= RVEModel.ConstrainedSketch(name='', sheetSize=,transform=)
#Define two point to creat a Circle
RVESketch1.CircleByCenterPerimeter(center=, point1=)
#Assign fiber sketch to the part
RVEPart.PartitionFaceBySketch(sketchUpEdge=,faces=, sketch=)
#Partition the fiber
RVEPart.PartitionCellByExtrudeEdge(line=,cells=,edges=,sense=)
return RVEPart

5.3. Definition of constituent materials (Create materials and section)
This function allows for the definition of each phase of the composite,
inclusion/reinforcement and matrix. In this study, two material types are available for
the composite. The matrix and inclusion are each defined in a local coordinate system
to account for the material orientations as appropriate. The material properties to be
entered are in terms of engineering constants, Young’s modulus and Poisson’ratio for

elastic behaviours. The material properties of the RVE are given in Table 1. Moreover,
the principal direction of the RVE is z-direction.
##Function of creating materials for RVE
def CreateMaterial(RVE_para,RVEModel ):
#Define the name of two materials
RVEModel.Material(name='')
RVEModel.Material(name='')
#Define the property of materials
RVEModel.materials[''].Elastic(table=)
RVEModel.materials[''].Elastic(table=)
##Function of creating sections for RVE
def CreateSection(RVE_para,RVEModel ,RVEPart):
#Initial parameter
#Define the name of the Section
RVEModel.HomogeneousSolidSection(name='',material='',thickness=)
RVEModel.HomogeneousSolidSection(name='',material='',thickness=)
#find the face
#Assign the property of the section1
RVEPart.SectionAssignment(region=, sectionName='')
RVEPart.MaterialOrientation(region=, orientationType=, axis=,additionalRotationType=,
localCsys=,fieldName='', stackDirection=)
#find the other face
#Assign the property of the section2
RVEPart.SectionAssignment(region=, sectionName='')
RVEPart.MaterialOrientation(region=, orientationType=, axis=,additionalRotationType=,
localCsys=,fieldName='', stackDirection=)

5.4. Key degrees of freedom and meshing (Create assembly)
Macroscopix/average strains ( x0 ,  y0 ,  z0 ,  xy0 ,  xz0 ,  yz0 ) appearing in boundary
conditions are physical entities which can be assigned by independent node numbers
and treated as ordinary nodes or Degrees of Freedoms (DOFs). The nodal
‘displacements’ (dimensionless) at these special DOFs give the corresponding
macroscopic strains directly, and eliminate the need to obtain them by averaging
strains from all elements. The macroscopic/average strains ( x0 ,  y0 ,  z0 ,  xy0 ,  xz0 ,  yz0 )
can be assigned to the ‘Boundary constraints’ which are related to the boundary
displacement of the RVE.
In the terminology of ABAQUS, these key DOFs have been introduced as
‘reference points’. The key DOFs representing the macroscopic strains appear in the
PBCs that they become physical entities as part of the RVE.
The FE analysis requires reasonable meshes for numerical convergence

considerations. For RVEs, extra restrictions must be imposed on meshes. To impose
periodic boundary conditions on a RVE, the paired faces corresponding to periodic
boundary conditions must be tessellated identically. This can usually be achieved by
copying one tessellated surface to another before generation of mesh inside the RVE
between paired faces. The mesh of the RVE may be appropriately formulated based
on rational considerations of boundary symmetries. In this study, the infusion is at the
center of the RVE, and it can be easy to keep the mesh of the faces periodicity through
the Python script without user involvement.

Fig.2 FE mesh of RVE model for the unidirectional composite

##Function of creating assembly for RVE
def CreateAssembly (RVE_para,RVEModel ,RVEPart):
#Initial parameter
RVEAssembly =RVEModel.rootAssembly
#Create a independent mesh Instance
RVEAssembly.Instance(name='', part=, dependent=)
# Creat nine reference points, and set a set
RVEAssembly.ReferencePoint(point=)
RVEAssembly.Set(referencePoints=, name=)
#Define two types of the mesh
elemType1 = mesh.ElemType(elemCode=,elemLibrary=)
elemType2 = mesh.ElemType(elemCode=,elemLibrary=)
#Assign two types of the mesh to the fiber and matrix respectively
RVEAssembly.setElementType(regions=, elemTypes=)
#creat mesh seed
RVEAssembly.seedEdgeBySize(edges=, size=)
#mesh the partInstances
RVEAssembly.generateMesh(regions=)
return RVEAssembly

5.5. The displayed field (Create step)
This function runs the ABAQUS to creat the analysis step of the RVE. In this
study, the outdata fields are strain (E), displacement (U), and stress (S).
##Create RVE step
def CreateStep (RVE_para,RVEModel ):
#Define a static step
RVEModel.StaticStep(name='', previous='')
#Define the outdata
RVEModel.fieldOutputRequests[''].setValuesInStep(stepName='',variables=)

5.6. Imposition of periodic boundary conditions (Create MPCs)
z
E

F

Wz
H

G
A

B

Wx
D

x

C

Wy

y

Fig.3 Prescribed PBC of 3-D case
There are three types of node sets: faces, edges and vertices. To avoid redundant
constraints on edges and vertices, the nodes on edges and vertices are excluded from
faces, and nodes on vertices are also excluded from edges when defining node sets in
order to impose PBCs. A vertex node set contains only one node, while each node set
for edges and faces could contain multiple nodes. The equations of all faces, edges,
and vertices are as follows:
The faces:
uEFGH  u ABCD  0
uHGCD  uEFBA  Wx  x0 uBCGF  u ADHE  0


0
vHGCD  vEFBA  0
vBCGF  v ADHE  Wy  y vEFGH  vABCD  0
w


0
 HGCD  wEFBA  0  wBCGF  wADHE  0  wEFGH  wABCD  Wz  z

(19)

The edges:
0
0
0
0 u
uCG  u AE  Wx  x0  Wy  xy
uHG  u AB  Wx x0  Wz  xz
FG  u AD  Wy  xy  Wz  xz




0
0 
0
0 
0
0
vCG  v AE  Wx yx  Wy  y vHG  v AB  Wx  yx  Wz  yz vFG  v AD  Wy  y  Wz  yz (20)

0
0 
0
0
0
0
 wCG  wAE  Wx  zx  Wy  zy  wHG  wAB  Wx  zx  Wz  z wFG  wAD  Wy  zy  Wz  z
0
0
0
0 u
uHD  uFB  Wx x0  Wy  xy
uDC  uEF  Wx  x0  Wz  xz
BC  u EH  Wy  xy  Wz  xz




0
0 
0
0 
0
0
vHD  vFB  Wx yx  Wy  y vDC  vEF  Wx  yx  Wz  yz vBC  vEH  Wy  y  Wz  yz (21)



0
0
0
0
0
0
 wHD  wFB  Wx zx  Wy  zy  wDC  wEF  Wx zx  Wz  z  wBC  wEH  Wy  zy  Wz  z

The vertices:
0
0 
0
0
uG  u A  Wx  x0  Wy  xy
 Wz  xz
uF  uD  Wx x0  Wy  xy
 Wz  xz



0
0
0 
0
0
0
vG  v A  Wx yx  Wy  y  Wz  xz vF  vD  Wx  yx  Wy  y  Wz  xz

0
0
0
0
0
0
 wG  wA  Wx  zx  Wy  zy  Wz  z  wF  wD  Wx  zx  Wy  zy  Wz  z

(22)

0
0 
0
0
uH  uB  Wx x0  Wy  xy
 Wz  xz
uC  uE  Wx  x0  Wy  xy
 Wz  xz



0
0
0 
0
0
0
vH  vB  Wx yx  Wy  y  Wz  xz vC  vE  Wx  yx  Wy  y  Wz  xz

0
0
0
0
0
0
 wH  wB  Wx zx  Wy  zy  Wz  z  wC  wE  Wx  zx  Wy  zy  Wz  z

(23)

In the terminology of ABAQUS, these equations have been introduced as
‘EQUATION’. The following is the steps in applying PBC on a model [20]:
(1) Group the node sets on faces, edges, and vertices of the RVE. The
*EQUATION keyword allows only those nodes on opposite faces which have
matching coordinates. Hence, regular or identical mesh on each the opposing faces are
required. In case of irregular mesh, the outermost layer of the RVE geometry need be
remeshed such that this condition is satisfied.
(2) The prescribed loading (displacement) is applied on a reference point situated

outside the RVE domain.
##Create RVE boundary conditons
def CreateBoudary (RVE_para,RVEModel ,RVEPart,RVEAssembly):
#Initial parameter
#find nodes, define MPC
RVEAssembly.Set(nodes=,name='')
RVEAssembly.Set(nodes=,name='')
RVEModel.Equation(name='',terms=)

5.7. Load case generation (Create load)
Making use of key DOFs, six macroscopic strains ( x0 ,  y0 ,  z0 ,  xy0 ,  xz0 ,  yz0 ) can
be prescribed in any way, individually or in any combination.
In the code, six macroscopic stresses are prescribed individually as six load cases
within a single step of analysis, so that uniaxial strain states are obtained, as required
in Eq. (23). The user has the option of generating an additional combined load case by
specifying the macroscopic strain components involved. Thus, three stiffness tensors

(C12 , C13 , C23 ) are obtained by the combined load case.
##Create the load
def CreateLoad (RVE_para,RVEModel ,RVEAssembly):
region1 = RVEAssembly.sets['']
RVEModel.DisplacementBC(name='', createStepName='',region=,u1=)

5.8. Generating a analysis work (Create job)
This function runs the ABAQUS to creat the analytical task of the RVE. The
display of the post-processing is set up, and the precision of the outdata is also set up
in this function. To save computation time, the number of the CPU may be added to
accelerate the computation.
##Create the job
def CreateJob (RVE_para):
# Inital parameter
RVEJob = mdb.Job(name=, model=)
# Submint job
RVEJob.submit()
#Wait the job to complete
RVEJob.waitForCompletion()

5.9. The computation of the stiffness tensor (Create post-processing)
This function outputs the analysis results of the RVE. With the macroscopic
strains being expressed in terms of ‘displacements’ applied to RVEs, all effective
properties of the RVEs are obtained in terms of the key DOFs. As the post-process of
averaging stresses and strains are not usually available directly from commercial FE
codes, the ABAQUS offers extremely useful means of simplifying the post-processing
of FE analysis for RVEs.
def CreateResult (RVE_para,stif):
# Inital parameter
odb = openOdb(path=)
history = odb.steps[''].historyRegions[''].historyOutputs['']
subHistory = history.data
# Averaged stress
#Total energy
RVEenergy = subHistory[len(subHistory)-1][1]
RVEstress=2*RVEenergy/(RVEvolume*RVEstrain)

5.10. Main Program (Create the initial parameter)
The main program defines the global execution flow. The main program locates
at the end of the source code since the Python interpreter needs the definition of all
other functions before they are called in the main program [19].
In the current main program, this program is followed by the function body that
can be partitioned into three blocks: input initialization, iteration loop, and compliance
matrix.
(1) Input initialization block: the initial design parametric model is introduced in
the following. Since the information in the design part, such as the strain information,
is frequently used in the program. Thus, all parameters are predefined referring to the
objects: the mesh seed, the length of the RVE, the fiber volume fraction, the type of
the element, the material properties of the fiber and matrix, and the prescribed strain.
A model database object is created referring to the input design model and assigned to
a structure variable RVE_para.
(2) Iteration loop: for the stiffness tensor, there are nine unknown parameters in
the matrix for the unidirectional composite. Hence, nine jobs are created to obtain the
unknown stiffness tensor parameters as shown in Figure 4. Except for MPCs and
Loads, other modules such as the geometry, material properties, and mesh are the

same to be created in this study. The MPCs and Loads modules are created based on

( x0  1;  y0 ,  z0 ,  xy0 ,  xz0 ,  yz0  0) . The load case

different load cases such as

( x0  1;  y0 ,  z0 ,  xy0 ,  xz0 ,  yz0  0) is to compute the stiffness tensor C11 .
(3) Compliance matrix: the compliance tensor is obtained by the inverse of the
stiffness tensor as shown in Eq. (24)
 S11
S
 12
S
 S    013

0

 0

S12

S13

0

0

S22

S23

0

0

S23

S33

0

0

0

0

S44

0

0

0

0

S55

0

0

0

0

0  C11 C12
0  C12 C22
0  C13 C23

0   0
0
0   0
0
 
S66   0
0

C13

0

0

C23

0

0

C33

0

0

0

C44

0

0

0

C55

0

0

0



0 
0 
0 
 (24)
0 
0 

C66 

Thus, the effective properties can be obtained as follows from the analysis of the
RVEs:

E1 

1
S11

v12  

S12
S11

G12 

1
S44

E2 

1
S22

v13  

S13
S11

G13 

1
S55

E3 

1
S33

v23  

S23
S 22

G23 

1
S66

(25)

Due to the unidirectional boron/aluminum RVE, predicted elastic properties are

E1  E2 , G13  G23 and v13  v23 . Moreover, it is noted that nine indendent material
constants must be determined for the RVE. Therefore, in this study, in total there are
nine ABAQUS analytical tasks to obtain the nine indendent material constants. Three
Poisson’ratios v12 , v13 ,and v23 are computed after the computation of E1 , E2 ,and E3 as
shown in Figure 4. To improve the computation time of obtaining each of the above
properties, there is a parallel method to simultaneously calculate the multiple material
constants. On the other hand, calculating them individually from the RVE analysis
and finding out the relationships among material constants can serve as a valuable
check on the RVE. In particular, the correct application of all the boundary conditions
for the RVE is important and a more sophisticated study.

Geometry, material
properties, and mesh

MPC,
Load

MPC,
Load

 C11 C12
C
C22
 12
C13 C23
 0
0

0
 0

0
 0

MPC,
Load

E1 

1
S11

v12  

S12
S11

G12 

1
2 S 44

E2 

1
S 22

v13  

S13
S11

G13 

1
2 S55

E3 

1
S33

v12  

S 23
S 22

G23 

1
2 S66

 S11
S
 12
 S13
0

0

0

S12

S13

0

0

S 22

S 23

0

0

S 23

S33

0

0

0

0

S 44

0

0

0

0

S55

0

0

0

0

Fig. 4 The computation flow of the code

0 
0 

0 
0 

0 

S66 

C13

0

0

C23

0

0

C33

0

0

0

C44

0

0

0

C55

0

0

0

0 
0 

0 
0 

0 

C66 

Obtained by
the inverse of
the stiffness
tensor

##========Main Program=======
##Set parameters and input
RVE_para={}
stif=zeros([6,6])
RVE_para["mesh"]=0.05
RVE_para["length"]=1.0000
RVE_para["volumeFraction"]=0.47
RVE_para["elementCode"]=C3D20
RVE_para["matrixProperties"]=(68.3,0.3)
RVE_para["fiberProperties"]=(379.3,0.1)
RVE_para["strain"]=0.0001
RVE_para["prestrain"]=[0.0000,0.0000,0.0000,0.0000,0.0000,0.0000]
for n in range(len(RVE_para["prestrain"])):
RVE_para["prestrain"]=[0.0000,0.0000,0.0000,0.0000,0.0000,0.0000]
RVE_para["number"]=n
RVE_para["name"]='RVE'+str(n+1)
RVEModel = mdb.Model(RVE_para["name"])
if (n<6):
RVE_para["prestrain"][n]=RVE_para["strain"]
else:
RVE_para["prestrain"]=[RVE_para["strain"],RVE_para["strain"],RVE_para["strain"],
0.0000,0.0000,0.0000]
k=8-n
RVE_para["prestrain"][k]=0.0000
RVEPart=CreatePart(RVE_para,RVEModel)
CreateMaterial(RVE_para,RVEModel)
CreateSection(RVE_para,RVEModel,RVEPart)
RVEAssembly=CreateAssembly (RVE_para,RVEModel,RVEPart)
CreateStep (RVE_para,RVEModel )
CreateBoudary (RVE_para,RVEModel,RVEPart,RVEAssembly)
CreateLoad (RVE_para,RVEModel,RVEAssembly)
CreateJob (RVE_para)
if (n<6):
stif[n][n]=CreateResult(RVE_para,stif)
if (5<n<8):
stif[0][n-5]=CreateResult(RVE_para,stif)
stif[n-5][0]=CreateResult(RVE_para,stif)
if (n==8):
stif[n-7][n-6]=CreateResult(RVE_para,stif)
stif[n-6][n-7]=CreateResult(RVE_para,stif)

print stif, print "change"
stif1=mat(stif)
print stif1
compliance=stif1**(-1)
print "compliance matrix",print compliance

6. Results
The visualization module of ABAQUS can be used to output von Mises’ stress as

shown in Figure 5. It is observed that the deformations are generally similar to the
results in the literature [18].

(a) ( x  1;  y ,  z ,  xy ,  xz ,  yz  0)

(b) ( y  1;  x ,  z ,  xy ,  xz ,  yz  0)

(c) ( z  1;  x ,  y ,  xy ,  xz ,  yz  0)

(d) ( xy  1;  x ,  y ,  z ,  xz ,  yz  0)

(e) ( xz  1;  x ,  y ,  z ,  xy ,  yz  0)

(f) ( yz  1;  x ,  y ,  z ,  xy ,  xz  0)

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

Figure 5 The von Mises’ Stress field corresponding to different loads.
The orthotropic elastic matrix is obtained by applying PBC to the RVE. The

results of the stiffness tensor are given as follows by using ABAQUS.
 11   0.006962
  -0.001779
 22  
 33  -0.000906
 
0
12  
 13  
0
  
0
 23  

  11 
  
-0.000906
0
0
0
  22 
  33 
0.004653
0
0
0
 
0
0.021864
0
0
  12 
  13 
0
0
0.018437
0
 
0
0
0
0.018434   23 

-0.001779 -0.000906
0.006962
-0.000906
0
0
0

0

0

0

It is seen from examining Table 3 that the predicted properties are generally in
good agreement with the results in the literature, and the experimental values.
.Table 3 Results and comparison for unidirectional boron/aluminum RVE (Vf =0.47)*
material constants
E3(Gpa) E2(Gpa) G23(Gpa) G12(Gpa)
V23
Present
213
143
53.8
45.4
0.194
Ref.13
214
143
54.2
45.7
0.195
Ref.14
215
144
57.2
45.9
0.19
Ref.21
214
135
51.1
0.19
Ref.22
214
156
62.6
43.6
0.20
Ref.23
215
123
53.9
0.19
Ref.24
215
135.2
53.9
52.3
0.195
Ref.25
216
140
52
0.29

V12
0.256
0.253
0.29
0.31
0.295
-

7. Conclusions
This study presents a Python code for computing effective properties of the RVE
under PBCs. With simple modifications, the basic Python code may be extended to
the computation of the effective properties of the more complex microstructure such
as several inclusions with the different shapes, different orientations and different
aspect ratios, and even the random distribution of inclusions. More importantly, the
code provides a convenient platform upon which further extensions such as the closed
loop optimization for the material and geometric parameter design of the composite
microstructure, and different optimizers could be also easily built. In doing so, the
user may experiment with various algorithms and tackle a wide range of problems.
The Python code is also presented for educational and engineering practice purposes.
Moreover, with extensions such as those presented in this study, and with the
ABAQUS FEA/modelling power, the Python code is capable of solving further
engineering design problems by using efficient optimization methods.
*

The principal direction of the RVE is z-direction in this study, and the principal direction of the RVE is
x-direction in the references.

Acknowledgements
This study was funded by Project of the Program of National Natural Science
Foundation of China (NSFC) grant number 11572120 and 61232014; Program for
New Century Excellent Talents in University under the grant number NCET-11-0131.

Appendix
1. A two-dimensional case
A 2-D RVE model is created. The model consists of a fiber reinforcement and
matrix, with a volume fraction of 50 %. The properties of the materials are given in
Table 4.
Material
Fiber
Matrix

Table 4 Material properties of fiber and matrix
E (Gpa)
72.5
2.6

v
0.22
0.4

The unidirectional RVE is assumed to be orthotropic and linearly elastic. For 2D
RVE, from Eq. (13), the constitutive relation of this effective material can be written
as

   C  

(26)

where C  is the stiffness matrix.

C11 C12
C   C12 C22
 0
0

0 
0 
C66 

(27)

After obtaining  ij for given  ij from the computation of a RVE, Cij can be
obtained from Eq. (26). The relation between the homogeneous strain energy and the
average strain is
U
1
1
1
 C1112  C121 2  C22 22  C66 62
VRVE 2
2
2

(28)

2. Imposition of periodic boundary conditions
C

D
Matrix
Y

Fiber

Wy

X

Wx

A

B

Fig. 6 Prescribed PBC of 2-D case
Periodic boundary conditions:
The edges:
0
uCD  u AB  Wy xy0

uBC  u AD  Wx x 

0 
0
v

v

W


BC
AD
x
yx


vCD  vAB  Wy y

(29)

0
0
0
0

uC  u A  Wx x  Wy xy 
uB  uD  Wx x  Wy xy

0
0 
0
0
v

v

W


W



C
A
y
y
x
yx

vB  vD  Wy y  Wx yx

(30)

The vertices:

3. Python implement of the two-dimensional case
##========Main Program=======
##Set parameters and input
RVE_para={}
stif=zeros([3,3])
RVE_para["mesh"]=0.05
RVE_para["length"]=1.0000
RVE_para["volumeFraction"]=0.5
RVE_para["elementCode"]='CPS3'
RVE_para["matrixProperties"]=(2.6,0.4)
RVE_para["fiberProperties"]=(72.5,0.22)
RVE_para["strain"]=0.0001
RVE_para["prestrain"]=[0.0000,0.0000,0.0000]
for n in range(len(RVE_para["prestrain"])+1):
RVE_para["prestrain"]=[0.0000,0.0000,0.0000]
RVE_para["number"]=n
RVE_para["name"]='RVE'+str(n+1)
RVEModel = mdb.Model(RVE_para["name"])
if (n<3):
RVE_para["prestrain"][n]=RVE_para["strain"]
else:
RVE_para["prestrain"]=[RVE_para["strain"],RVE_para["strain"],0.0000]
RVEPart=CreatePart(RVE_para,RVEModel)
CreateMaterial(RVE_para,RVEModel)
CreateSection(RVE_para,RVEModel,RVEPart)
RVEAssembly=CreateAssembly (RVE_para,RVEModel,RVEPart)
CreateStep (RVE_para,RVEModel )
CreateBoudary (RVE_para,RVEModel,RVEPart,RVEAssembly)
CreateLoad (RVE_para,RVEModel,RVEAssembly)
CreateJob (RVE_para)
if (n<3):
stif[n][n]=CreateResult(RVE_para,stif)
if (n==3):
stif[n-3][n-2]=CreateResult(RVE_para,stif)
stif[n-2][n-3]=CreateResult(RVE_para,stif)
print stif
print "change"
stif1=mat(stif)
print stif1
compliance=stif1.I
print "compliance matrix"
print compliance

References
1.

Ogierman, W., & Kokot, G. (2015). Modeling of constitutive behavior of anisotropic
composite material using multi-scale approach. Mechanics, 21(2), 118-122.

2.

Sen, O., Davis, S., Jacobs, G., & Udaykumar, H. S. (2015). Evaluation of convergence
behavior of metamodeling techniques for bridging scales in multi-scale multimaterial
simulation. Journal of Computational Physics, 294(C), 585-604.

3.

Zhou, X. Y., Gosling, P. D., Pearce, C. J., Ullah, Z., & Kaczmarczyk, L. (2015).
Perturbation-based stochastic multi-scale computational homogenization method for woven
textile composites. International Journal of Solids & Structures, 80, 368-380.

4.

Kouznetsova, V., Brekelmans, W. A. M., & Baaijens, F. P. T. (2001). An approach to
micro-macro modeling of heterogeneous materials. Computational Mechanics, 27(1), 37-48.

5.

Michel, J. C., Moulinec, H., & Suquet, P. (1999). Effective properties of composite materials
with periodic microstructure: a computational approach. Computer methods in applied
mechanics and engineering, 172(1), 109-143.

6.

Würkner, M., Berger, H., & Gabbert, U. (2014). Numerical investigations of effective
properties of fiber reinforced composites with parallelogram arrangements and imperfect
interface. Composite Structures, 116, 388-394.

7.

Xia, Z., Zhou, C., Yong, Q., & Wang, X. (2006). On selection of repeated unit cell model and
application of unified periodic boundary conditions in micro-mechanical analysis of
composites. International Journal of Solids and Structures, 43(2), 266-278.

8.

P. Suquet. (1987). Elements of homogenization theory for inelastic solid mechanics.
Sanchez-Palencia E, Zaoui A, Homogenization techniques for composite media, Berlin:
Springer-Verlag.

9.

Bakhvalov, N. S., & Panasenko, G. P. (1984). Homogenization in periodic media,
mathematical problems of the mechanics of composite materials.

10. Sánchez-Palencia, E. (1980). Non-homogeneous media and vibration theory.
11. Bensoussan, A., Lions, J. L., & Papanicolaou, G. (2011). Asymptotic analysis for periodic
structures (Vol. 374). American Mathematical Soc.
12. Tchalla, A., Belouettar, S., Makradi, A., & Zahrouni, H. (2013). An ABAQUS toolbox for
multiscale finite element computation. Composites Part B: Engineering, 52, 323-333.
13. Sun, C. T., & Vaidya, R. S. (1996). Prediction of composite properties from a representative
volume element. Composites Science and Technology, 56(2), 171-179.
14. Xia, Z., Zhang, Y., & Ellyin, F. (2003). A unified periodical boundary conditions for
representative volume elements of composites and applications. International Journal of
Solids and Structures, 40(8), 1907-1921.
15. Xia, Z., Ju, F., & Sasaki, K. (2007). A general finite element analysis method for balloon
expandable stents based on repeated unit cell (RUC) model. Finite Elements in Analysis and
Design, 43(8), 649-658.
16. Xu, K., & Xu, X. W. (2008). Finite element analysis of mechanical properties of 3D
five-directional braided composites. Materials Science and Engineering: A, 487(1), 499-509.
17. Lubineau, G., & Ladeveze, P. (2008). Construction of a micromechanics-based intralaminar

mesomodel, and illustrations in ABAQUS/Standard. Computational Materials Science, 43(1),
137-145.
18. Yuan, Z., & Fish, J. (2008). Toward realization of computational homogenization in practice.
International Journal for Numerical Methods in Engineering, 73(3), 361-380.
19. Zuo, Z. H., & Xie, Y. M. (2015). A simple and compact Python code for complex 3D
topology optimization. Advances in Engineering Software, 85, 1-11.
20. Shahzamanian, M. M., Tadepalli, T., Rajendran, A. M., Hodo, W. D., Mohan, R., Valisetty, R.,
& Ramsey, J. J. (2014). Representative volume element based modeling of cementitious
materials. Journal of Engineering Materials and Technology, 136(1), 011007.
21. Chamis, C. C. (1983). Simplified composite micromechanics equations for hygral, thermal
and mechanical properties.
22. Sun, C. T., & Chen, J. L. (1991). A micromechanical model for plastic behavior of fibrous
composites. Composites Science and Technology, 40(2), 115-129.
23. Riley, M. B., & Whitney, J. M. (1966). Elastic properties of fiber reinforced composite
materials. Aiaa Journal, 4(9), 1537-1542.
24. Hashin, Z., & Rosen, B. W. (1964). The elastic moduli of fiber-reinforced materials. Journal
of applied mechanics, 31(2), 223-232.
25. Kenaga, D., Doyle, J. F., & Sun, C. T. (1987). Characterization of boron/aluminum composite
in the nonlinear range as an orthotropic elastic-plastic material. J. Compos. Mater.;(United
States), 21.

