This documentation is part of the IMTEK Mathematica Supplement (IMS)

Homogenization of a Composite Material with Spherical Inclusions

José Luis Gómez-Muñoz

Visit my home page


In this example we will calculate the effective thermal conductivity for a composite material. The inclusions are spherical and they form a periodic structure:


The Asymptotic Homogenization Method AHM (N. Bakhalov and G. Panasenko, Homogenisation Averaging Processes in Periodic Media, Kluwer Academic Publishers, 1989, pages 215 and 216) can be used to calculate the effective thermal conductivity, Overscript[k,^] of this composite material. It is given by:

Overscript[k,^]=∫_ (-0.5)^0.5∫_ (-0.5)^0.5∫_ (-0.5)^(0, 5)k(x,y,z)∂ M_1/∂ xdxdydz  

Here x, y and z expand the periodic cell, k(x,y,z) is the conductivity in each material,


k_i in the inclusion
k_m in the matrix

and, for isotropic materials, M_1 is a solution of the problem:

here [[u]] _CONTACT means the difference (contrast) between the values of u at the two sides of the contact between the matrix ant the inclusion. Therefore [[u]] _CONTACT=0 means that u is continuous from the fiber to the matrix (ideal contact).
FEM is going to be used in order to calculate the effective conductivity Overscript[k,^] for the case k_m=(
1 0 0
0 1 0
0 0 1
  and k_i=(
10 0 0
0 10 0
0 0 10
  and an inclusion concentration of 0.3

Outline of the calculation

First we will create a domain for a periodic celd of the composite, see the graph at the top of this document. This domain is an imsNexus made of some nodes connected by flat polygons (called domain segments). These nodes expand the boundaries and contact region of the periodic cell. The domain nexus is then used as an input for the mesher program Tetgen. The output of this program is read as another nexus, this time made of more nodes which are connected by tetrahedrons (called elements). These new nodes expand the three dimensional periodic cell. Then the FEM stiffness matrix is built from this last nexus and the boundary conditions. The corresponding linear system is solved. FEM is used again to calculate the flux, and the x-component of the flux is integrated in the periodic cell in order to obtain the effective conductivity.

Simulation data

We turn off the warning  about a new symbol that has a name that is similar to the name of an existing symbol, because several variables have similar names in this document.  


Off[General :: spell1] ;

Off[General :: spell] ;

First we define the physical and geometrical data we are going to work with:


inclusionConductivity = 10. ;

matrixConductivity = 1. ;

inclusionConcentration = 0.3 ;

Now we define the computational parameters. The value curveDensity is the number of points per unit of area in the frontier (contact region) between the inclusion and the matrix. If you do not want to see the intermidiate graphics of this calculation then set makeGraphics=False. The mesher program Tetgen will be used to generate the FEM mesh. The qSwitch especifies the maximum radius-edge ratio and the aSwitch the maximum volume of the tetrahedral elements generated by the program Tetgen.


curveDensity = 100 ;

qSwitch = "1.3" ;

aSwitch = "0.001" ;

makeGraphics = True ;

IMS Nexus for the contact region and the bounding box

The calculation status will appear at the lower left corner of this notebook


Needs["Imtek`ShowStatus`"] ;

imsShowStatus["IMS Nexus for the contact region and the bounding box"] ;

There are several different options to generate a sphere in Mathematica, for example ContourPlot3D[], ParametricPlot3D[] or even SurfaceOfRevolution[]. In this example, the Sphere[] coomand of the Graphics`Shapes package is used in order to generate the sphere (contact region).


Needs["Imtek`DomainElementLibrary`"] ;

Needs["Graphics`Shapes`"] ;

radius = 0.5 * (6 * inclusionConcentration)/π^(1/3) ;

sphereQuality = Ceiling[1. + 3.54491 * (-0.0795775 + curveDensity * radius^2)^(1/2)] ;

inclusionGraphics = RotateShape[Graphics3D[Sphere[radius, sphereQuality, sphereQuality]], 0, Pi/2, Pi/2] ;

If[makeGraphics, Show[inclusionGraphics]] ;


The bounding box is generated in a similar way. The bounding box and the spherical inclusion make the domainGraphics.


Needs["Imtek`DomainElementLibrary`"] ;

Needs["Graphics`Polyhedra`"] ;

boundingBoxGraphics = Polyhedron[Hexahedron, {0, 0, 0}, 1/2^(1/2)] ;

domainGraphics = Show[boundingBoxGraphics, inclusionGraphics, PlotRange→All, DisplayFunction→Identity] ;

If[makeGraphics, Show[domainGraphics, PlotRange→ {{-0.51, .51}, {-.51, .51}, {-.51, 0.01}}, Boxed→False, DisplayFunction→ $DisplayFunction]] ;


The command imsGraphics3DToNexus is used to transform the domainGraphics to an imsNexus made of imsNode and imsDomainSegment objects. In 3D, an imsDomainSegment is actually a flat polygon. The graph is stored in the unmarkedDomainNexusGraphics in order to use it later.   


unmarkedDomainNexus = imsGraphics3DToNexus[domainGraphics] ;


Inserting Region Markers in the Nodes and Segments


imsShowStatus["Inserting Region Markers in the Nodes and Segments"] ;

Here we define the constants and progams (newMarker[] and centerOfMassCoords[]) that will be used to insert region markers in the nodes and segments of the contact region and the bounding box.  


leftSide = 1 ;

rightSide = 2 ;

frontSide = 3 ;

backSide = 4 ;

bottomSide = 5 ;

topSide = 6 ;

contactRegion = 7 ;

Needs[ "Imtek`Polygon`" ] ;

Here we create the lists of marked nodes and segments. Remember that in 3D, an imsDomainSegment is actually a flat polygon connecting some nodes (incidents).    


Needs["Imtek`Graph`"] ;

unmarkedNodes = imsGetNodes[unmarkedDomainNexus] ;

markedNodes = Map[imsSetMarkers[#, newMarker[imsGetCoords[#]]] &, unmarkedNodes] ;

unmarkedSegments = imsGetElements[unmarkedDomainNexus] ;

markedSegments = Map[imsMakeDomainSegment[imsGetIds[#], imsGetIncidentsIds[#], newMarker[centerOfMassCoords[#, unmarkedDomainNexus]]] &, unmarkedSegments] ;

Creating the marked domain Nexus


imsShowStatus["Creating the marked domain Nexus"] ;

Here the nodes and segments, with the proper markers, are combined to make the domain nexus that will be used as an input for the mesher program


interiorNodes = Select[markedNodes, imsGetMarkers[#] == contactRegion&] ;

boundaryNodes = Select[markedNodes, imsGetMarkers[#] ≠  contactRegion&] ;

domainNexus = imsMakeNexus[boundaryNodes, interiorNodes, markedSegments] ;

Here we can see the domain nexus.




Material Markers for Tetgen


imsShowStatus["Material Markers for Tetgen"] ;

The materialDomainMarkers will let the mesher program "know" that inside the contact region we have one kind of material and outside we have another material. We can also have diferent maximum volume constrains for each material, but these will take effect only if the mesher program Tetgen is executed using its -a swith without a number after it.


matrixMarker = 123 ;

matrixMarkerPoint = {0.49, 0.49, 0.49} ;

matrixMaxElementVolume = 0.001 ;

inclusionMarker = 678 ;

inclusionMarkerPoint = {0.01, 0.01, 0.01} ;

inclusionMaxElementVolume = 0.001 ;

materialDomainMarkers = {{matrixMarkerPoint, matrixMarker, matrixMaxElementVolume}, {inclusionMarkerPoint, inclusionMarker, inclusionMaxElementVolume}} ;

Creating the input file for the mesher program Tetgen


imsShowStatus["Creating the input file for the mesher program Tetgen"] ;

The program Tetgen is going to be used to create the mesh. You can download Tetgen from

Here we define the filenames and the directory where the input and output will be stored. Notice that if the files already exist, then they are erased, so that old calculations do not interfere with this calculation.


SetDirectory[ $HomeDirectory ] ;

myFileName = "composite3D" ;

Map[If[MemberQ[FileNames[], #], DeleteFile[#]] &, myFileNameList] ;

Now we create the input file for Tetgen using the domainNexus and the materialDomainMarkers.


holes = {} ;

Needs["Imtek`Interfaces`TetgenInterface`"] ;

meshFile = imsToTetgenInputFile[ domainNexus , holes, materialDomainMarkers] ;

Export[  polyInputFileName, Chop[meshFile], "Table"] ;

Run the program Tetgen


imsShowStatus["Run the program Tetgen"] ;

Here we run the program Tetgen. It is assumed that the binary file tetgen is available. In a Unix system that means that it is located in your binary files directory (for example /usr/local/bin). On the other hand, in a Windows system it might be easier to have a copy of tetgen.exe in the current directory (the same folder where the input file was created). You can download Tetgen from . If everything goes right, then the result of the Run[] command must be zero. A number (aSwitch) is specified after the -a switch, therefore the maximum volume constrains that were specified in the materialDomainMarkers will not take effect. The aSwitch is the maximum Tetrahedron volume. The -A switch tells Tetgen to propagate the material markers to the corresponding tetrahedrons. The -p switch indicates that the input is a .poly file. The number (qSwitch) after the -q switch controls the quality of the tetrahedrons.


tetgenResult = Run[ "tetgen -pAa"<>aSwitch<>" -q"<>qSwitch<>" "<> polyInputFileName<>" > "<>logFileName ] ;

Print[Import[logFileName, "Text"]] ;


If[tetgenResult≠0, imsShowStatus["Tetgen Error"] ; Quit[]] ;

Now we verify that the output files (.1.ele, .1.node and .1.face) have been created by Tetgen




{composite3D.1.ele, composite3D.1.face, composite3D.1.node, composite3D.poly, composite3D.txt}

Reading Tetgen's output and inserting boundary conditions


imsShowStatus["Reading Tetgen's output and inserting boundary conditions"] ;

Tetgen has created more nodes inside the three-dimensional domain. Those nodes are not connected anymore by flat polygons (segments). Now they are connected by tetrahedral elements. This information is stored in the .node and the .ele file. Besides, the .face file has information about those nodes that belong to each face (either in the bounding box or the contact region). We will read all this information into a new nexus. First we read the information of the faces:


tetgenFaces = Drop[Cases[Import[faceOutputFileName, "Table"], {_Integer, ___}], 1] ;

leftSideNodesIds = Union[Flatten[Select[tetgenFaces, Last[#] == leftSide&][[All, {2, 3, 4}]]]] ;

rightSideNodesIds = Union[Flatten[Select[tetgenFaces, Last[#] == rightSide&][[All, {2, 3, 4}]]]] ;

frontSideNodesIds = Union[Flatten[Select[tetgenFaces, Last[#] == frontSide&][[All, {2, 3, 4}]]]] ;

backSideNodesIds = Union[Flatten[Select[tetgenFaces, Last[#] == backSide&][[All, {2, 3, 4}]]]] ;

bottomSideNodesIds = Union[Flatten[Select[tetgenFaces, Last[#] == bottomSide&][[All, {2, 3, 4}]]]] ;

topSideNodesIds = Union[Flatten[Select[tetgenFaces, Last[#] == topSide&][[All, {2, 3, 4}]]]] ;

contactRegionNodesIds = Union[Flatten[Select[tetgenFaces, Last[#] == contactRegion&][[All, {2, 3, 4}]]]] ;

The markerFunction will be used to assign the proper marker to each node of the new nexus, using the info that was obtained from the face file


The dataFunction and the valueFunction will be used in order to store the boundary conditions in the nodes of the new nexus.


defaultData = "defaultData" ;

defaultValue = 123456 ;

Next we read the output files into the nexus tetgenMesh. In the moment of reading, the marker function is used to mark those nodes that belong to a face, and the valueFunction and dataFunction store the boundary conditions in the corresponding nodes. Remember that in this new nexus the nodes are connected by tetrahedrons.


femMesh = imsReadTetgenOutput[ { nodeOutputFileName, eleOutputFileName}, markerFunction, valueFunction, dataFunction ] ;

Here you can see a wireframe representation of the tetrahedrons in the tetgenMesh nexus. They are separated based on the kind of material (the marker for an element shows the type of material, on the other hand, the marker for a node shows if the node belongs to a face):




Assembly of the global stiffness matrix


imsShowStatus[ "Assembly of the global stiffness matrix"] ;

Here we define the thermal conductivity. It is a function of the element marker, which identifies the type of material


Here we create the lists allEmptyMatrices, allElements and allElementNodes that will be used to build the stifness matrices.


elementDim = Length[ imsGetIncidentsIds[ imsGetElements[ femMesh, 1 ] ] ] ;

elementSMEmpty = Table[ 0., { elementDim }, { elementDim }] ;

elementRHSEmpty = Table[ 0., { elementDim }, { 1 }] ;

allElements = imsGetElements[femMesh] ;

allElementNodes = Map[imsGetNodes[femMesh, #] &, imsGetIncidentsIds[allElements]] ;

allRows = allCols = imsGetIds[allElementNodes] ;

Needs["Imtek`Assembler`"] ;

allEmptyMatrices = Map[{imsMakeElementMatrix[elementSMEmpty, #, #], imsMakeElementMatrix[elementRHSEmpty, #, {1}]} &, allRows] ;

Here we create all the element stiffness matrices


Needs["Imtek`FEMOperators`"] ;

Timing[allFilledMatrices = Thread[ imsNFEMDiffusion[allEmptyMatrices, allElements, allElementNodes, conductivity ]] ; ]


{4.947 Second, Null}

Finally the element stiffness matrices are assembled into the global stifness matrix


Timing[numberOfNodes = Length[imsGetNodes[femMesh]] ; stiffness = SparseArray[ {}, {numberOfNodes, numberOfNodes}, 0. ] ; imsAssemble[ Transpose[ allFilledMatrices ][[1]], stiffness ] ;]


{0.521 Second, Null}

Boundary conditions in the global stiffness matrix


imsShowStatus["Boundary conditions in the global stiffness matrix"] ;

Now we insert the boundary conditions in the stiffness matrix. The Neumann conditions are equal to zero, therefore we do not need to insert them. The continuity conditions are "ideal contact", therefore we do not need to insert them either. We only insert the Dirichlet conditions


Needs["Imtek`BoundaryConditions`"] ;

load = Table[0., {numberOfNodes}, {1}] ;

flattLoad = Flatten[load] ;

imsDirichlet[{stiffness, flattLoad}, pos, values, 1] ;

load = Partition[flattLoad, 1] ;

Solve the linear system


imsShowStatus["Solve the linear system"] ;

The variable solution holds the M_1values for every node:


Timing[solution = LinearSolve[ stiffness, load ] ; ]


{0.03 Second, Null}

Plot of the solution


imsShowStatus["Plot of the solution"] ;

First we load the required package


Needs["Imtek`UnstructuredPlot`"] ;

Here is a 3D plot of the solution:


<<Graphics`Legend` ;


Derivative of M_1


imsShowStatus["Derivative of M1"] ;

Now we need to calculate the partial derivative of M_1with respect to x. One way to do it is using shape functions for derivatives: M_1=Underscript[∑, j]N_jφ_j^*, this leads to the system ∫_ΩN_iN_jdΩ φ_j^*=∫_ΩN_iN_jdΩ M_1_j which we will solve with the imsNFEMConvection and imsNFEMTransientMatrix operators.


convectionX = Function[ { marker, x, y , z}, { { 1., 0. , 0.} } ] ;

Here we create all the element stiffness matrices



{8.773 Second, Null}

The element stiffness matrices are assembled into the global stifness matrix


Needs["Imtek`Assembler`"] ;


{1.251 Second, Null}

The global matrices massMatrix and stiffnessX build together a system of equations which is solved.


derivX = LinearSolve[massMatrix, Flatten[ stiffnessX . solution ] ] ;

Graph the derivative


imsShowStatus["Graph the derivative"] ;

We graph the derivative of M_1with respect to x_1




Integral of the x-component of the flux


imsShowStatus[Integral of the  x -component of the flux] ;

Finally we integrate the x-component of the flux in order to obtain the effective coefficient

Overscript[k,^]=∫_ (-0.5)^0.5∫_ (-0.5)^0.5∫_ (-0.5)^(0, 5)k(x,y,z)∂ M_1/∂ xdxdydz

First we load the required packages.


Needs["Imtek`ShapeFunctions`"] ;

Needs["Imtek`MeshElementLibrary`"] ;

(* Needs["Calculus`Integration`"] ; *)

We integrate inside each element and add. In order to integrate, we map to the unit element (also called master element) multiply by the conductivity and the Jacobian and integrate.



{6.589 Second, 2.07966}

Final answer


imsShowStatus[m = "e="<>ToString[Length[allElements]] <>" k="<>ToString[effectiveConductivity] ] ;

Our composite material is equivalent to a homogeneous material with conductivity:





This is an approximated answer, it can be improved if the mesh is refined by increasing the parameter curveDensity. Finally we restore the warning messages that were turned off at the begining of this notebook.


On[General :: spell1] ;

On[General :: spell] ;


[1] N. Bakhalov and G, Panasenko, Homogenisation: Averaging Processes in Periodic Media, Kluwer Academic Publishers, 1989
[2] Daryl L. Logan, A First Course in the Finite Element Method, Third Edition, Brooks/Cole Thomson Learning, 2002

Created by Mathematica  (April 28, 2006) Valid XHTML 1.1!