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

Homogenization of a Composite Material with Spherical Inclusions

By

José Luis Gómez-Muñoz

Visit my home page

Introduction

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, of this composite material. It is given by:

=k(x,y,z)dxdydz

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

k(x,y,z)=

in the inclusion | |

in the matrix |

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

here means the difference (contrast) between the values of u at the two sides of the contact between the matrix ant the inclusion. Therefore =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 for the case =(

1 | 0 | 0 |

0 | 1 | 0 |

0 | 0 | 1 |

10 | 0 | 0 |

0 | 10 | 0 |

0 | 0 | 10 |

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.

In[1]:=

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

In[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 http://tetgen.berlios.de/ 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.

In[6]:=

IMS Nexus for the contact region and the bounding box

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

In[10]:=

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).

In[12]:=

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

In[18]:=

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.

In[23]:=

Inserting Region Markers in the Nodes and Segments

In[25]:=

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.

In[26]:=

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).

In[36]:=

Creating the marked domain Nexus

In[41]:=

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

In[42]:=

Here we can see the domain nexus.

In[45]:=

Material Markers for Tetgen

In[46]:=

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.

In[47]:=

Creating the input file for the mesher program Tetgen

In[54]:=

The program Tetgen is going to be used to create the mesh. You can download Tetgen from http://tetgen.berlios.de/

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.

In[55]:=

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

In[59]:=

Run the program Tetgen

In[63]:=

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 http://tetgen.berlios.de/ . 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.

In[64]:=

In[66]:=

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

In[67]:=

Out[67]=

Reading Tetgen's output and inserting boundary conditions

In[68]:=

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:

In[69]:=

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

In[77]:=

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

In[78]:=

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.

In[82]:=

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):

In[83]:=

Assembly of the global stiffness matrix

In[84]:=

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

In[85]:=

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

In[86]:=

Here we create all the element stiffness matrices

In[94]:=

Out[95]=

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

In[96]:=

Out[96]=

Boundary conditions in the global stiffness matrix

In[97]:=

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

In[98]:=

Solve the linear system

In[104]:=

The variable solution holds the values for every node:

In[105]:=

Out[105]=

Plot of the solution

In[106]:=

First we load the required package

In[107]:=

Here is a 3D plot of the solution:

In[108]:=

Derivative of

In[110]:=

Now we need to calculate the partial derivative of with respect to x. One way to do it is using shape functions for derivatives: ∇=, this leads to the system dΩ =∇dΩ which we will solve with the imsNFEMConvection and imsNFEMTransientMatrix operators.

In[111]:=

Here we create all the element stiffness matrices

In[112]:=

Out[112]=

The element stiffness matrices are assembled into the global stifness matrix

In[113]:=

Out[114]=

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

In[115]:=

Graph the derivative

In[116]:=

We graph the derivative of with respect to

In[117]:=

Integral of the x-component of the flux

In[118]:=

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

=k(x,y,z)dxdydz

First we load the required packages.

In[119]:=

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.

In[121]:=

Out[121]=

Final answer

In[122]:=

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

In[123]:=

Out[123]=

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.

In[124]:=

References

[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) |