Homat

$\renewcommand{\vec}[1]{\boldsymbol{#1}}$

Introduction¶

Multi-scale analysis of components with a heterogeneous microstructure¶

Material properties characterize the functionality of almost any product. For many industrial applications the final material properties result from a number of process steps. Each of them allows the tailoring of the microstructure in order to achieve the desired functionality. Bridging the scales is an essential feature of the multi-scale analysis of components with locally different material properties ¹. In the case of multi-scale material modeling and in homogenization in general, one proceeds from the detailed description of all material heterogeneities at the lower scale upward in order to obtain effective material properties accurately. Alternatively, in process modeling it is important to be able to step down to the scale of the real, heterogeneous material microstructure. The corresponding technique is known as the localization method. This coupling between the scales corresponds to a serial, hand-shake type exchange of data (e.g. stresses, strains, etc.) passing up and down the scales.

In Fig. 1 the multi-scale analysis approach is outlined, using the example of the U-O forming production step of an X65 steel pipeline tube with a pearlite-ferrite microstructure. The corresponding effective thermal properties (thermal conductivity, capacity, density) and thermo-elastic properties (Hooke matrix, Young and shear moduli, Poisson ratios, thermal expansion coefficient, etc.) are evaluated by the homogenization tool HOMAT by applying the asymptotic homogenization method to heterogeneous materials. The theoretical background of this method is briefly summarized in Appendix Homogenization of Materials with Periodic Microstructure.

Figure 1: Multi-scale simulation of the U-O forming of pipeline steel with pearlite-ferrite microstructure.

The effective non-linear materials properties, like the elasto-plastic flow curves, are determined via uniaxial virtual testing by applying, at the moment, the commercial finite element software package ABAQUS.

Then, the simulation of the U-O forming process of the pipeline tube is performed at the macro-scale and critical areas are detected. Locally in such a critical area, the macro-strain states are extracted and applied as boundary condition to the ferrite-pearlite microstructure in order to determine the corresponding distribution of micro-strains and stresses on the Representative Volume Element (RVE) and the location of possible failure. The governing equations of the localization procedure within in the framework of the asymptotic homogenization method, are briefly outlined in Appendix Post-processing of macro-scale results: the localization step.

The homogenization and virtual testing tools are based on an accurate geometrical description of the complex microstructure (i.e. the topology of the individual phases and their spatial distribution) and the properties of the individual phases. The strength of these methods is their ability to consider transverse isotropic, cubic, orthotropic or generally anisotropic material properties. The phase properties may be either given by ab-initio calculations or provided by models, experiments and/or databases. The geometry of the microstructure can be either simulated directly by the multi-component phase-field program MICRESS or by the 3-D cellular automata SphaeroSim for semi-crystalline thermoplastics ^2,3,4, or specified experimentally via 3-D CT scans or SEM micrographs ⁵. As mentioned in Appendix Homogenization of Materials with Periodic Microstructure, not the entire microstructure is discretized at the micro-scale but only a RVE. For a material with periodic microstructure its size is defined by the period while for materials with random microstructure 1-D spectral analysis or 2-point probability functions are used ⁵ to specify a first guess of the RVE size. A more detailed description about the adopted stochastic homogenization procedure for materials with random microstructure is given in the bookchapter "Prediction of effective properties" ¹.

The homogenization method further requires the automatic generation of interface elements with zero thickness at the boundaries between phases, grains, spherulites or material domains. These interface elements allow the calculation of implicit forces/fluxes induced by the discontinuity of the material properties at these boundaries. To achieve this task a dedicated interface program, named Mesh2HOMAT, is supplied (see Fig. 2).

Figure 2: Overview of the prediction of effective material properties of a simulated MICRESS microstructure.

The pre-processor tools Mesh2HOMAT and Mesh2ABAQUS¶

The interface program Mesh2XXX has a generic name and it includes both specific interface tools: Mesh2HOMAT and Mesh2ABAQUS. The aim of this interface program is twofold. Firstly, it allows the conversion of a structured or unstructured mesh from one data format to another. Four different input file formats are currently available: the standard platform exchange format VTK, the universal format UNV of IDEAS, the Abaqus input format INP ⁶ and the MICRESS format, for which the possibility to read the compressed grain information directly is retained.

This preprocessor automatically detects the interfaces between adjacent material domains, grains or spherulites and generates new nodes and the corresponding interface elements. It evaluates the outward normal to the free surfaces and defines the pair of periodic surfaces as a function of the user specifications. Depending on the actual problem, either 1-D, 2-D or, by default, 3-D periodic boundary conditions can be specified on the unit cell. The microstructure solutions of MICRESS and SphaeroSim, available in the initial VTK format, are rewritten in a VTK file for the newly generated, unstructured finite element mesh of the RVE. Eventually, Mesh2Homat writes a geometry file PROBLEM_NAME.hgeo and a template of the HOMAT command file PROBLEM_NAME.hin.

By specifying Mesh2ABAQUS as program name, this interface program generates automatically the six (or three in 2-D) input files necessary to perform uniaxial virtual tests of the discretized RVE of the considered microstructure within ABAQUS. The program automatically generates the required periodic boundary conditions. It also reads the orientation of each grain of the simulated microstructure by MICRESS (PROBLEM_NAME.TabO file) and transforms these orientations into a format readable by ABAQUS. The detailed description of the common and specific commands of both interface program versions Mesh2HOMAT and Mesh2ABAQUS is given in a separate manual.

Overview and Features of the Homogenization Tool HOMAT¶

Survey of HOMAT¶

The homogenization tool HOMAT is based on the asymptotic homogenization method applied to materials with periodic microstructures as outlined in Appendix Homogenization of Materials with Periodic Microstructure. This method allows the determination of the following effective properties:

• thermal conductivity, heat capacity and density
• Darcy permeability
• elastic stiffness (Hooke) matrix, thermal expansion coefficients, Young and shear moduli, Poisson ratios and volumetric eigenstrains

HOMAT started with the homogenization of cooled multi-layer systems by the prediction of effective thermal properties of these systems ⁷. Only periodic boundary conditions are specified on the RVE, even for materials with random microstructure. For these materials the periodicity is applied in a weak manner: nodes without a projection point in the same phase on the periodic face are eliminated from the periodic boundary condition. This is mostly relevant for Stokes flow homogenization ⁸. The tool evaluates microscopic periodic displacements on the unit cell, which are used to express the unknown effective properties. Depending on the problem nature, up to 8 (thermo-elasticity with eigenstrains) or 3 (thermal and fluid flow) boundary value problems are discretized by finite elements and solved on the unit cell successively. For the mechanical homogenization, the rigid body modes are fixed automatically and the user can specify the kind of the desired homogenization analysis: 3-D, plane strain or plane stress. The right-hand sides of these problems correspond to the material discontinuities at the grain or material boundaries plus the variation of the property inside each grain/domain (e.g. with the temperature or/and orientation of the domain). The evaluated effective properties of the unit cell are written in an ASCII file PROBLEM_NAME.hres and in a material data file (.mat) of ABAQUS type. Hence, the effective material properties can be directly transferred to the macro-scale simulation programs (ABAQUS, DEFORM, SYSWELD, etc.). Moreover, for each right-hand side and component, the distributions of the implicit fluxes, forces or boost forces and of the corresponding microscopic displacements on the unit cell are written in files of VTK format. In presence of mixed finite elements (modelling incompressible materials) also the periodic pressure fields are provided in VTK format. Since version 6.0, it is possible in thermoelastic homogenizations to write also the microscopic strain distribution on the RVE. In the case of Stokes flow homogenization, the periodic velocity and pressure fields are written by default. In order to reduce the storage data, these files are, by default, binary ones if the microstructure is simulated by MICRESS or SphaeroSim. These microscopic results can then be visualized with Paraview.

In the present user guide all the keywords of the input file PROBLEM_NAME.hin and of the associated geometry file PROBLEM_NAME.hgeo are detailed.

The program HOMAT is launched easily by typing:

homat PROBLEM_NAME

without any postfix on the command line or by dragging the command file PROBLEM_NAME.hin on the icon/executable of HOMAT. Note that if you execute HOMAT without any model information on the command line, the program will display usage information:

Usage: homat_name <problem_name> <options>
<problem_name> : name of the homat input file without postfix
<options>
   -check : check of the consistency of the hin input file

The new option -check on the command line permits to check the consistency of the introduced input file PROBLEM_NAME.hin. If activated, HOMAT will check if the required geometry and declared material properties and the input files (Microsim and/or Macro) are present. Moreover, it prints the chosen homogenization or localization analysis, the type of applied periodicity and some mesh informations, like the total number of nodes, the number of pure geometrical nodes and the number of additional interface nodes.

A detailed description how to use HOMAT in an efficient way is outlined and different homogenization and localization examples are provided with this version. Each example has its own definition file, allways named Def_Example.txt. All the information required to run the considered example are outlined there. The specification of the command file PROBLEM_NAME.hin for a discretized RVE of the considered heterogeneous material can be divided in different parts, which all have to be specified for a homogenization run:

• The definition of the thermal and/or thermo-elastic material properties of each phase or domain of the microstructure have to be specified. A specific data file, named MATERIAL.hmat, has to be created for each phase or material. The material properties specific keywords will be described in Sec. Definition of the Materials Properties.
• The geometry of the RVE is defined in the PROBLEM_NAME.hgeo file, created by Mesh2HOMAT. This file contains a list of the nodes, the bulk elements, the interface elements defined between phases or domains and the element faces belonging to outer surfaces. The geometry and model related keywords are presented in Sec. Definition of the Microstructure Model.
• In Sec. Homat input files, the specific HOMAT keywords of the command file PROBLEM_NAME.hin are specified. This chapter involves also the description of the definition of localization runs in selected Gauss points of the preceding macro-simulation.
• Finally, in Sec. HOMAT Output, in function of the specified homogenization analysis, the different kinds of output files are detailed and examples provided. The different solution files of a localization analysis are also explained.

Features of the HOMAT Version 6.0¶

In this section the key features of the HOMAT Version 6.0 are briefly outlined. The corresponding commands and input data are detailed in Sec. Definition of HOMAT- specific Commands, dedicated to the specific HOMAT commands.

• HOMAT automatically evaluates the periodic displacement fields on the RVE (for example, 6 for a purely mechanical homogenization in 3-D) required by the calculation of the desired effective properties.
• Several homogenizations can be defined and performed in one HOMAT run. For example, for a given microstructure RVE effective properties can be evaluated at several specified temperatures in the same HOMAT simulation.
• In the presence of a heterogeneous multi-layer material, like the cooled combustion chamber wall composed of a superalloy, a bondcoat and a thermal barrier coating (see Fig. 25), the user can choose between multi-layer or mono-layer homogenization. In the case of mono-layer homogenization, the user has to specify the layer for which effective properties are determined.
• Different kind of homogenizations (thermal, thermo-elastic, ...) can be performed in the same HOMAT simulation. Note that the new Stokes flow homogenization cannot be combined with the two other kinds of homogenization.
• The periodic microscopic displacement fields $\chi^{\text{rs}}$, solutions of the linear problems [see Eq. \eqref{eq:solu_2}] and the implicit fluxes (thermal homogenization) and/or implicit forces (thermo-elastic homogenization) are saved in files written in VTK format. In presence of mixed finite elements, used to accurately handle material incompressibility, the predicted microscopic pressure fields are also written in VTK format. In Stokes flow homogenization, the microscopic velocity and pressure fields are provided. The different solution files of HOMAT can be visualized by Paraview on the initial, undeformed unit cell or RVE.
• The considered periodic fields on the RVE (temperature, displacement, pressure, etc) are discretized in 3-D by linear or quadratic finite elements of hexahedral, pentahedral and/or tetrahedral type: Since version 6.0, the initial element geometry can be also curved (using quadratic ansatz functions).
• In order to model accurately the quasi-incompressibility of phases like the melt or the amorphous phase of a semi-crystalline polymer, several mixed element types have been implemented in HOMAT 6.0: namely the equal order linear displacement and pressure elements (P1P1 and Q1Q1) stabilized by the double Gauss point integration, proposed by Li and He ^9,2, the mini element version, where the displacement field is augmented by a volume bubble. As a tetrahedron, this element discretization is automatically stabilized but not its hexahedral variant. Eventually, new stabilized elements with volume bubble and Li-He stabilization have been developed and tested successfully.
• The number of Degree of Freedoms (DOF) per node is now variable in the RVE model. If a quasi-incompressible phase is present in combination with compressible ones, a 3-D model will have nodes with 4 DOFs as well as other nodes with classical 3 DOFs.
• The elements can be either fully or reduced integrated ¹⁰. By default, the reduced integration scheme is adopted. The user can specify the number of Gauss points per element to define the integration type.
• The predicted effective properties are written in a PROBLEM_NAME.hres file. This file serves for further post-processing the homogenization results and also to print them in a convenient way. Since version 6.0, the user can write the effective properties also in a specified local axis system, e.g. In the grain axis system. Moreover, the effective material properties are written in a PROBLEM_NAME.mat file of ABAQUS type. Thus, these predicted material properties can be directly included in the finite element model of the specimen at the macro-scale.
• The mechanical phase or domain properties can be either isotropic, transverse isotropic, cubic, orthotropic, monoclinic or anisotropic. In thermoelasticity, the thermal expansion can be either isotropic or orthotropic; whereas the eigenstrain is assumed to be isotropic. The thermal phase and domain properties can be now isotropic, transverse isotropic, orthotropic or anisotropic.
• If the phases or domains of the considered RVE are cubic with cubic axes corresponding to the RVE axis system, the effective Young and shear moduli, the effective Poisson and Zener anisotropy coefficients of the heterogeneous materials are expressed and written in the PROBLEM_NAME.hres result file.
• Dependent on the nature of the homogenization problem, the user can choose in thermo-elasticity to perform either a 3-D homogenization (default case) or a 2-D one of plane strain or stress type. In the case of thermal or fluid flow homogenization, either 3-D or 2-D homogenization can be selected.
• On the RVE 3-D periodic boundary conditions (B.C.) are, by default, specified. Dependent on the nature of the heterogeneous materials of the multi-layer system, 2-D or 1-D periodic B.C. can be also specified.
• HOMAT automatically detects in the list of the periodic nodes, the double nodes on wedges, which built the intersection of periodic surfaces, and removes the rigid body modes. The user has no specific fixations to define! Starting from a node on one wedge, a special link procedure is developed for 2-D and 3-D periodicity, like shown in Fig. 3. As the finite element mesh on one surface of a periodic pair can be different from the mesh on the other surface, we obtain a chain of linked D.O.F, which must be open in order to ensure equal movement of the linked nodes. In presence of 2-D periodicity the corner nodes of the lowest face are fixed, like shown in Fig. 3 (left), whereas in the case of 3-D periodicity all 8 corner nodes are fixed and three chains of linked wedge nodes are build (Fig. 3, right).

Figure 3: Definition on the RVE of the open chains of linked nodes in the case of 2-D or 3-D periodicity respectively. The fixed nodes are there marked by a cross.

• At the interface between two material domains, tied contact is achieved. The slave nodes are condensed, only the master nodes are retained in the linear system. Thus, the system to be solved is reduced. This procedures leads to save computational time. After resolution, the results are duplicated on the condensed slave nodes.
• In the version 6.0 the concept of master-slave nodes is extended to the periodic surfaces- Independent, if periodic surfaces have directly corresponding nodes or not, the slave nodes are eliminated and replaced, via a projection algorithm, detailed in Appendix Periodic boundary conditions via a master-slave projection algorithm, either by the directly corresponding master node or by the master nodes of the projection element. This way, the total system to be solved is reduced further.
• A major difference between the version 6.0 and the previous versions 5.x is the new treatment of the periodic B.C. via a master-slave projection algorithm, which doesn't use any penaty factor. This way, no additional stiffness at periodic surface nodes is introduced leading to larger microscopic displacement variations over the RVE and, subsequently, to larger corrections of volume averaging material property (see Eqs. \eqref{eq:eff_hoke_mat} and \eqref{eq:eff_eigenstress}). Therefore, the effective material properties predicted with HOMAT version 6.0 differ from previous predictions (version 5.3).
• The CSR sparse matrix format is now used to store the global system matrix. In this format, flexible number of DOFs per node are more easily treated and the implementation, in near future, of direct solvers is facilitated.
• An improved iterative BiCGstab solver is used to solve the linear systems of the mechanical, thermal or flow problem. The solver is enhanced by an automatic restart option and a new symmetric pre-conditioning, based on a scaling of the diagonal terms. This scaling is mainly important for the calculation of the displacement field $\psi$ (calculation of the effective thermal expansion coefficients).
• The application of the periodicity B.C. on a heterogeneous material with random microstructure, like an open-cell metallic foam, is delicate. For each periodic pair, the nodes belonging on the same phase or material domain on both surfaces are linked classically. But, nodes not having a projection point/element in the same phase on the master periodic surface are eliminated from the periodic boundary condition. This way, weak periodic B.C. are specified. When evaluating of the Darcy permeability, a special box can be defined inside the RVE where the velocities are used in the averaging process ¹¹. Thus, the main affected boundary layers near the free periodic surfaces can be omitted in the effective permeability evaluation.

Homat input files¶

Common features of HOMAT input files¶

A HOMAT input file is a plain ASCII data file. A template for the parameters of the keywords of the command file PROBLEM_NAME.hin is created by the pre-processor Mesh2HOMAT. The user has to verify if the default options are adequate for the considered homogenization or localization run. He can modify and adapt the predefined parameters by using any text editor. The input files consist of a series of lines containing keywords, parameters and data. The following rules are applicable to all keywords and data lines of the command, geometry and material files read by HOMAT:

• The first non-blank character of each keyword line must be the slash symbol "/".
• The keyword must be followed by a comma "," if any parameter is specified.
• Parameters must be separated by commas.
• Blanks on any line are ignored.
• A line is always limited to 256 characters, including blanks.
• Any line beginning with the hash symbol "#" is a comment line. Such lines can be placed anywhere in the input files of HOMAT, they are ignored by the “read line interpreter”.
• Keywords and parameters are not case sensitive and be either written in upper or lower case letters, except file names.
• If a parameter has a value, the equal sign "=" is used. The value can be either an integer, a float number or a character string, depending on the parameter's meaning.
• Continuation of a keyword line is sometimes necessary: for example, in presence of a large number of parameters of the considered keyword. If the last character on a keyword line is a comma, the next line is interpreted automatically as a continuation line of the previous line.
• Character strings can be built of up to 80 characters and are not case sensitive.
• Labels such as group or surface names are case insensitive. Their unique name has max. 80 characters and must not contain blanks.
• The name of the keywords and their parameters have been chosen as similar as possible to the keywords of the ABAQUS input file .inp. It is clear that HOMAT has its own keywords but mainly the keywords describing the material and geometry features of the model are made compatible to the ABAQUS input file. An easy change of the slash "/" into "*" allows an easy transformation of a HOMAT model in an ABAQUS model.

Definition of the Materials Properties¶

In this chapter the thermal and thermo-elastic properties per phase or per domain are presented respectively. Per property class they are listed in alphabetical order of the keywords. All these temperature dependent material properties are collected, per material, in a special file, named MATERIAL_NAME.hmat. Note that properties of the mushy zone can depend on the fraction solid, $f_{S}$, as the temperature is in non-unique above $T_{sol}$ due to recalescence phenomena ³. This file can be included easily in the command file of a specific homogenization application via the /INCLUDE command explained more in detail in Sec. Definition of HOMAT- specific Commands. This inclusion happens after the definition of the material by the keyword /MATERIAL and its parameter. This way, the same phase or constituent material properties are specified once for several homogenization runs.

/MATERIAL, NAME = name, TYPE = type

This keyword is part of the model definition of the heterogeneous material at the micro-scale level. It is used to indicate the start of a material definition. A material data block can be defined directly in the command file PROBLEM_NAME.hin or imported via the /INCLUDE command by loading an existing MATERIAL_NAME.hmat file, involving all properties of the considered phase or material. A material data block is defined by all keywords between the /MATERIAL line and either another /MATERIAL line or a keyword line that does not specify material properties. If a property is specified more than once for the same material, always the last definition will be used. Example:

/MATERIAL, NAME=ferrite, TYPE=ISO

The /MATERIAL keyword has two required parameters, the NAME of the constituent material or phase, like ferrite, austenite, ... and the TYPE of the specified material model. The name of the material must be defined and is an alias for the defined material properties. The parameter NAME will be used in the /SOLID SECTION options in order to assign the material to a group of elements (see the definition of these keywords in Sec. Definition of the Microstructure Model). The second parameter TYPE is essential not only for thermo-elasticity but, since version 5.3, also for the thermal homogenization. It specifies the material model used for the considered phase or material domain. Indeed, if the user specifies TYPE=ISO, which means that the ferrite phase is considered in the homogenization run as an isotropic material, only the Young modulus and the Poisson coefficient are considered in this simulation even if the cubic elastic constants are also provided in the material properties file .hmat or in the command file .hin. To specify the bcc ferrite phase as a cubic symmetric phase, characterized by the elastic constants $H_{11}$, $H_{12}$ and $H_{44}$, the user has to specify the material model ORTHO. Considering the elastic cubic symmetric behavior of the ferrite phase will, of course, improve the homogenization results obtained with the isotropic model. Therefore, the correct choice of the appropriate model is essential and has a direct impact on accuracy of the homogenization results. It is important to note that, by specifying TYPE=ORTHO for the ferrite phase, the present Young and shear modules and Poisson coefficients in the material file .hmat are ignored because the material is assumed to be of cubic symmetry in the considered homogenization run. Following material models can be specified via the TYPE parameter:

• TYPE=ISO: isotropic material behavior, defined only by the Young modulus and the Poisson coefficient and a single value for the thermal conductivity;
• TYPE=TR_ISO: transverse isotropic material behavior, requiring the introduction of 5 independent material coefficients for the elastic behavior and two coefficients for thermal conductivity;
• TYPE=ORTHO: orthotropic material behavior, requiring usually the definition of 9 independent material quantities and three values (or two in 2-D homogenization) for the thermal conductivity;
• TYPE=ANISO: anisotropic material behavior, requiring the definition of 21 material coefficients and 6 values for the thermal conductivity tensor;
• TYPE=SPHERU: spherulite model for semi-crystalline thermoplastics. This model is detailed in the recent publication ³. Please contact ACCESS, if you want to use this special model in combination with SphaeroSim results.

Remark: The cubic symmetry is a subclass of the orthotropic material behavior and requires only the introduction of three material constants: User friendly introduction of this material behavior can be achieved via the /ELASTIC keyword (see Sec. Thermo-elastic properties). The specification of material models with mechanical anisotropy but thermal isotropy is delicate. In a simulation including thermal and thermoelastic homogenization, the material has to be specified either TR_ISO, ORTHO or ANISO. Subsequently the definition of the thermal conductivity has to be in agreement with the adopted material behavior and also be specified as TR_ISO, ORTHO or ANISO phase. On the other side, if only a pure thermal homogenization will be performed, this phase or material can be specified as isotropic and only one thermal conductivity value is specified. In this case, only isotropic mechanical are accepted in the material file MATERIAL_NAME.hmat. For example, the presence of an orthotropic or anisotropic Hooke matrix definition will lead to an error message. Therefore, it is recommended to specify all phase properties to the anisotropy level of its thermoelastic behavior.

Note: The thermal conductivity definition of some materials has been modified in some existing examples in order to respect the agreement with the specified material model. It concerns Nicrofer wires in the hollow_grid example, which are assumed to be transverse isotropic, the Yttrium stabilized Zirconium thermal barrier of the S30-F example, which is orthotropic and the sapphire fiber, which is anisotropic.

Remark: The calculation of the effective Darcy permeability tensor uses only unitary viscosities and densities per fluid domain see ⁸, the Stokes flow homogenization in HOMAT requires for each specified /SOLID SECTION keyword the name of a defined material. Of course, the real mechanical or thermal anisotropy of the specified material is neglected in this analysis.

Thermal properties¶

• /CONDUCTIVITY, TYPE = type

  This keyword defines the heat conductivity of the considered phase or material. Since the HOMAT version 5.3, the mandatory parameter type specifies the nature of the thermal conductivity: isotropic, transverse isotropic, orthotropic or anisotropic. As already mentioned, the definition of the thermal conductivity has to be specified in agreement with the adopted material model (see keyword /MATERIAL). All values of the thermal conductivity may be temperature or fraction solid dependent. As HOMAT is units independent, the user can specify any type of coherent unit system (SI, CGS). As millimeter is defined as default length, classically the unit for thermal conductivity is W/(mm K). Depending on the phase nature, following different thermal conductivity definitions are available via the TYPE parameter:

  • TYPE=ISO: isotropic thermal behavior
    After the 1^st line with the keyword /CONDUCTIVITY and the parameter TYPE, follow at the next lines the isotropic thermal conductivity with the corresponding temperature or fraction solid. The temperatures or fraction solids must be introduced in an ascending manner. Below the 1^st temperature, here 0° C, the thermal conductivity is set to the 1^st value, here 0.009 W/(mm.K). Above T=1000°C (last value), it is set to 0.2205 W/(mm.K). Between two introduced temperatures, HOMAT interpolates linearly.

    Note: The above detailed definition of a temperature dependent material property is general and will not be repeated by each property except in special cases like the elastic constants.

    #
    /CONDUCTIVITY, TYPE=ISO
    0.009000, 0.0
    0.009120, 20.0
    0.010180, 100.0
    0.011500, 200.0
    ............
    0.022050, 1000.0
    #
    # lambda in W/(mm * K)
    #
    /DENSITY
    0.0087, 0.0
    #
    # temperature in C
    # Rho in g/mm**3
    #
    /HEAT CAPACITY, TYPE=ISO
    0.003370, 20.0
    0.003500, 200.0
    0.003880, 400.0
    ................
    0.007130, 1000.0
    #
    # rhocp = density * specific heat in J/(mm**3 *K)
    # Data provided by John Fernihough (1997)
    #
    

    Example 1: Temperature dependent thermal properties of the single crystal CMSX-4 superalloy.
    - TYPE=TR_ISO: transverse isotropic thermal behavior
    The transverse isotropic case is characterized by a plane of isotropy. Only in the direction perpendicular to this plane, the material presents a different behavior than in the plane. Classical examples are unidirectional fibers. In order to specify a transverse isotropic thermal behavior, the user has only to introduce two independent thermal conductivities: $\lambda_{\text{t}}$, $\lambda_{\text{p}}$, and the temperature $T$ or fraction solid $f_{\text{S}}$.
    $$\begin{equation*} \lambda_{\text{t}}\,, \lambda_{\text{p}}\,, T \end{equation*}$$

    Note that the subscript t stands for “transverse” and subscript p for “plane”. Example 2 shows such a definition of transverse isotropic thermal conductivity for a Nicrofer wire.

    /CONDUCTIVITY, TYPE=TR_ISO
    0.0107, 0.0113, 20.
    0.0117, 0.0127, 100.
    0.0133, 0.0144, 200.
    0.0149, 0.0160, 300.
    0.0164, 0.0176, 400.
    0.0180, 0.0192, 500.
    0.0194, 0.0206, 600.
    0.0209, 0.0222, 700.
    #
    # unit W/(mm*K)
    #
    

    Example 2: Transverse isotropic thermal behavior of a Nicrofer wire.
    - TYPE=ORTHO: orthotropic thermal behavior
    A material presents an orthotropic behavior if we can identify three orthogonal symmetry planes in a local axis system ($e_{1}$, $e_{2}$, $e_{3}$). In this local axis system, the thermal conductivity tensor of such a phase has only diagonal terms: $\lambda_{11}$, $\lambda_{22}$, $\lambda_{33}$ and the off-diagonal terms are zero. To describe such a thermal material behavior, the user introduces:
    $$\begin{equation*} \lambda_{11}\,, \lambda_{22}\,, \lambda_{33}\,, T \end{equation*}$$

    The Example 3 describes the definition of the orthotropic thermal conductivity of an yttrium stabilized zirconium oxide layer ⁵.

    /CONDUCTIVITY, TYPE=ORTHO
    0.001436, 0.001436, 0.001436, 20.
    0.001397, 0.001397, 0.001397, 200.
    0.001379, 0.001379, 0.001379, 600.
    0.001359, 0.001359, 0.001359, 800.
    0.001385, 0.001385, 0.001385, 1000.
    #
    # unit W/(mm*K)
    #
    

    Example 3: Orthotropic thermal behavior of a yttrium stabilized zirconia oxide.
    - TYPE=ANISO: anisotropic thermal behavior
    In the full anisotropic case, the material has no specific symmetry plane, whereas in the monoclinic case, it has only one symmetry plane. The thermal conductivity tensor $\lambda_{ij}$ of an anisotropic material in 3-D is usually a full $3\times3$ matrix. Due to its symmetry, this tensor has only six independent conductivity terms:
    $$\begin{equation} \label{eq:lambda_ten} \lambda = \ \begin{bmatrix} \lambda_{11} & \lambda_{12} & \lambda_{13} \\ \bullet & \lambda_{22} & \lambda_{23} \\ \bullet & \bullet & \lambda_{33} \\ \end{bmatrix} \end{equation}$$

    To introduce such an anisotropic thermal behavior of the material in a local axis system ($e_{1}$, $e_{2}$, $e_{3}$), the user has to specify the six conductivity terms plus eventually the temperature in following manner:

    $$\begin{equation*} \lambda_{11}\,, \lambda_{22}\,, \lambda_{33}\,, \lambda_{12}\,, \lambda_{13}\,, \lambda_{23}\,, T \end{equation*}$$

    As example, the effective thermal behavior of pearlite, which is an eutectoid phase mixture and composed of alternate ferrite and cementite layers, is outlined for the 18CrNiMo7-6 steel. These effective anisotropic thermal conductivities are results of a thermal homogenization run at the nanoscale, where a bi-lamella of ferrite and cementite is homogenized. Please look to the example ferr_pearl and the reference ¹² for more details.

    /CONDUCTIVITY, TYPE=ANISO
    3.29137E-02, 5.54883E-02, 4.96613E-02, -4.8956E-04, -2.4484E-04, 3.88469E-03, 27.
    ##
    #
    #unit W/(mm*K)
    #
    

    Example 4: Effective anisotropic thermal behavior of the eutectoid pearlite phase of the 18CrNiMo7-6 steel grade.

    Note that this definition of anisotropic thermal conductivity differs from the Abaqus definition.

  Remarks: If the thermal behavior of the material is not isotropic, then the keyword /ORIENTATION has to be specified in order to define the directions of material orthotropy or anisotropy. Alternative keyword definitions that HOMAT interprets as /CONDUCTIVITY are /THERMAL CONDUCTIVITY and /HEAT CONDUCTIVITY.
  

• /DENSITY

  This keyword allows the definition of the mass density $\rho$ of the considered material. Usually the mass density is specified in g/mm³. This thermo-physical property can be defined as temperature dependent:

  $$\begin{equation*} \rho\,, T \end{equation*}$$

  In Ex. 1 of the superalloy CMSX-4, a constant density (value = 0.0087 g/mm³) is defined for all temperatures.
  

• /HEAT CAPACITY

  The user has two definition possibilities: either he specifies the heat capacity or the specific heat per unit mass of the material. Both physical quantities are linked by the simple relation:

  $$\begin{equation} \text{heat capacity (rcp)} = \text{density (rho)} * \text{specific heat (cp)} \end{equation}$$

  In Ex. 1 of the superalloy CMSX-4, the temperature dependence of the heat capacity is specified and defined in the reference unit system J/(mm³ K). In presence of phase changes it takes the corresponding latent heat into account and a strong variation of its value is observed. The user specifies

  $$\begin{equation*} \text{rcp}\,, T \end{equation*}$$

• /SPECIFIC HEAT

  The specific heat allows an alternate definition of the heat capacity of the material. The specific heat per mass unit of the material is defined as:

  $$\begin{equation} \label{eq:spec_heat} c_{p}\left( T \right) = \frac{\text{dU}(T)}{\text{dT}}\,, \end{equation}$$

  where $U(T)$ is the internal energy density of the material per unit mass due to thermal effects only.

  Remarks: The user must choose: either he specifies the heat capacity of the material or the specific heat but not both physical quantities for the same phase or constituent material at the same time. HOMAT evaluates the effective density by volume averaging on the RVE, heat capacity and specific heat of the heterogeneous material. The emissivity of the phase or material does not need to be specified. HOMAT neglects, the radiation exchange inside the discretized unit cell. The user specifies

  $$\begin{equation*} c_{p}\,, T \end{equation*}$$

Thermo-elastic properties¶

• /ELASTIC, TYPE = type

  This option is used to define the elastic properties of a phase or a constituent material. Its linear thermo-elastic behaviour is defined by the generalized Hooke law, linking the stress tensor $\sigma_{ij}$ to the total elastic strain tensor $\varepsilon_{kl}$:

  $$\begin{equation} \label{eq:stress} \sigma_{ij} = H_{ijkl}(T)\left[ \varepsilon_{kl}(\vec{u}) - \alpha_{kl}(T)\Delta T - \kappa_{kl} \right]\,. \end{equation}$$

  The Hooke matrix $H_{ijkl}(T)$ is a fourth order elasticity tensor and the thermal expansion coefficient $\alpha_{kl}(T)$ are both temperature or fraction solid dependent quantities. The volumetric eigenstrains $\kappa_{kl}$ are isotropic, per phase, and generally present in polycrystalline materials. The keyword /ELASTIC outlines the definition of the Hooke matrix. Depending on the number of symmetry planes for the elastic properties, a phase or constituent material can be classified as isotropic (infinite number of symmetry planes) or anisotropic (no symmetry planes). Some materials have a restricted number of symmetry planes. For example, orthotropic materials have three orthogonal symmetry planes for the elastic properties and monoclinic materials only one symmetry plane. The number of independent components of the Hooke matrix depends on such symmetry properties. The level of anisotropy and the mode of defining the elastic properties are chosen by the TYPE parameter, which must be always specified on the keyword line. If the material is not isotropic, then the /ORIENTATION option has to be specified in order to define the directions of anisotropy.

  • TYPE=ISO: isotropic elasticity
    The elastic properties are completely defined by specifying only the Young's modulus, E, and the Poisson coefficient $\nu$. For each temperature²⁶ $T$ the user provides one line of data

    $$\begin{align*} E, \nu, T \end{align*}$$

    Indeed, the shear modulus G of an isotropic elastic material is given by

    $$\begin{equation} \label{eq:shear_mod_iso} G = \frac{E}{2(1 + \nu)}\,. \end{equation}$$

    These two parameters can be specified as function of the temperature (see Example 5 of the ferrite phase).

    Remark: If the elastic properties of a constituent material are specified isotropic then the thermal expansion coefficient must be also defined as isotropic.

    /ELASTIC, TYPE=ISO
    211670.0, 0.27, 20.
    208220.0, 0.27, 95.
    ..........
    137210.0, 0.27, 760.
    95300.0, 0.27, 911.
    #
    # unit MPa
    #
    /EXPANSION, TYPE=ISO
    1.28925E-5, 20.
    1.31538E-5, 100.
    ................
    1.72930E-5, 850.
    1.74279E-5, 911.
    #
    

    Example 5: Elastic properties and thermal expansion coefficient of an assumed isotropic ferrite phase.
    - TYPE=TR_ISO: transverse isotropic case
    A special class of orthotropic material behaviour is the transverse isotropic case, which is characterized by a plane of isotropy at every point of the material. Assuming the plane 1-2 is the plane of isotropy, transverse isotropy requires $E_{11} = E_{22} = E_{\text{p}}$ and $\nu_{12} =\nu_{\text{p}}$ and $G_{\text{p}} = \frac{E_{\text{p}}}{2(1 + \nu_{\text{p}})}$, where the subscript p stands for “in-plane”; but also $G_{13} = G_{23} = G_{\text{t}}$ and $\nu_{31} = \nu_{32} = \nu_{\text{tp}}$ and $\nu_{13} = \nu_{23} = \nu_{pt}$ and $G_{13} = G_{23} = G_{\text{t}}$, where subscript t stand for “transverse” to the plane of isotropy. Thus, while $\nu_{\text{tp}}$ is the Poisson ratio that characterizes the strain in the plane of isotropy induces by a stress normal to it, $\nu_{pt}$ characterizes the strain in the normal direction to the plane of isotropy resulting from a in-plane stress. In general these quantities are not equal. They are related by:
    $$\begin{equation} \label{eq:nu_triso} \frac{\nu_{\text{tp}}}{E_{\text{t}}} = \ \frac{\nu_{\text{pt}}}{E_{\text{p}}} \end{equation}$$

    The corresponding compliance matrix $\mathbb{S}$, inverse of the Hooke matrix $\mathbb{H}$, has the following expression in the Voigt notation

    $$\begin{equation} \mathbb{S} = \mathbb{H}^{-1} = \begin{bmatrix} \frac{1}{E_{\text{p}}} & - \frac{\nu_{\text{p}}}{E_{\text{p}}} & - \frac{\nu_{\text{tp}}}{E_{\text{t}}} & 0 & 0 & 0 \\ - \frac{\nu_{\text{p}}}{E_{\text{p}}} & \frac{1}{E_{\text{p}}} & - \frac{\nu_{\text{tp}}}{E_{\text{t}}} & 0 & 0 & 0 \\ - \frac{\nu_{\text{tp}}}{E_{\text{t}}} & - \frac{\nu_{\text{tp}}}{E_{\text{t}}} & \frac{1}{E_{\text{t}}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{\text{p}}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{\text{t}}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G_{\text{t}}} \\ \end{bmatrix}\,. \end{equation}$$

    For each temperature²⁶ the user introduces only five independent material properties

    $$\begin{align*} E_{\text{t}}, E_{\text{p}}, \nu_{\text{tp}}, \nu_{\text{p}}, G_{\text{t}}, T \end{align*}$$

    to specify the transverse isotropic material behavior, like shown in Example 6 for a unidirectional reinforced glass-epoxy composite.

    /ELASTIC, TYPE=TR_ISO
    43565.0, 12800., 0.29, 0.307, 4143.85, 20.
    #
    # unit MPa
    #
    /EXPANSION, TYPE=ORTHO
    4.51477E-6, 1.02361D-05, 1.02361D-05, 20.
    #
    

    Example 6: Tranverse isotropic elastic behaviour of a UD reinforced glass-epoxy composite.
    - TYPE=ENG: orthotropic elasticity specified by engineering constants
    Linear elasticity of an orthotropic material is most easily defined by giving its engineering constants: three Young moduli $E_{1}$, $E_{2}$, $E_{3}$; the Poisson coefficients $\nu_{12}$, $\nu_{13}$, $\nu_{23}$ and the three shear moduli $G_{12}$, $G_{13}$ and $G_{23}$ associated to the material's or phase's principal material directions. This set of 9 independent moduli defines following elastic compliance matrix
    $$\begin{equation} \mathbb{S} = \ \mathbb{H}^{\mathbf{- 1}} = \begin{bmatrix} \frac{1}{E_{1}} & - \frac{\nu_{21}}{E_{2}} & - \frac{\nu_{31}}{E_{3}} & 0 & 0 & 0 \\ - \frac{\nu_{12}}{E_{1}} & \frac{1}{E_{2}} & - \frac{\nu_{32}}{E_{3}} & 0 & 0 & 0 \\ - \frac{\nu_{13}}{E_{1}} & - \frac{\nu_{23}}{E_{2}} & \frac{1}{E_{3}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{12}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{13}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G_{23}} \\ \end{bmatrix}\,. \end{equation}$$

    The Poisson coefficient $\nu_{ij}$ measures the transverse strain in $j$ direction, when the material is stressed in $i$ direction. The symmetry of the Hooke and compliance matrices implies

    $$\begin{equation} \label{eq:poisson_rat} \frac{\nu_{\text{ij}}}{E_{i}} = \frac{\nu_{\text{ji}}}{E_{j}}\,. \end{equation}$$

    Of course, the engineering constants can also be given as functions of the temperature, if necessary. For each temperature²⁶, the user provides two lines

    $$\begin{align*} &E_{1}, E_{2}, E_{3}, \nu_{12}, \nu_{13}, \nu_{23}, G_{12}, G_{13},\\ &G_{23}, T \end{align*}$$

    as shown in Example 7 for an yttrium stabilized zirconium oxide layer ⁵.

    #
    /ELASTIC, TYPE=ENG
    5.096191E+04, 3.124513E+04, 6.37023875E+04, 2.705990E-01, 3.442350E-01, 3.434799E-01, 1.386069E+04, 1.386069E+04,
    1.386069E+04, 700.
    5.077967E+04, 3.121674E+04, 6.34745875E+04, 2.698298E-01, 3.424140E-01, 3.416010E-01, 1.374312E+04, 1.374312E+04,
    1.374312E+04, 750.
    ..............
    5.011533E+04, 3.099163E+04, 6.26441625E+04, 2.682793E-01, 3.369597E-01, 3.351691E-01, 1.351601E+04, 1.351601E+04,
    1.351601E+04, 900.
    4.975417E+04, 3.087421E+04, 6.21927125e+04, 2.673659E-01, 3.347734E-01, 3.331901E-01, 1.349206E+04, 1.349206E+04,
    1.349206E+04, 950.
    #
    /EXPANSION, TYPE=ORTHO
    9.579838E-06, 8.633990E-06, 8.215011E-08, 700
    9.278161E-06, 9.050553E-06, 5.160571E-07, 750
    ..................................
    9.449730E-06, 9.183234E-06, 3.867535E-07, 900
    9.380581E-06, 9.157542E-06, 5.460644E-07, 950
    #
    

    Example 7: The orthotropic elastic behaviour of yttrium stabilized zirconia oxide specified by engineering constants.

    As the elastic behaviour is orthotropic, the expansion coefficient of the material must be also specified as orthotropic.

    Note: An alternative and totally equivalent definition of the elastic orthotropic material behaviour can be achieved via the combined definition of orthotropic Young and shear moduli and of orthotropic Poisson coefficients (please see the corresponding keywords and Ex. 12).

  • TYPE=ORTHO: orthotropic elasticity via terms of the elastic Hooke matrix
    Linear elasticity of an orthotropic material can also be defined by giving directly the nine independent elastic stiffness terms, as functions of the temperature, if desired. In this case, the Hooke matrix, written in Voigt notation, has following expression

    $$\begin{equation} \mathbb{H} = \begin{bmatrix} H_{11} & H_{12} & H_{13} & 0 & 0 & 0 \\ H_{21} & H_{22} & H_{23} & 0 & 0 & 0 \\ H_{31} & H_{32} & H_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & H_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & H_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & H_{66} \\ \end{bmatrix}\,. \end{equation}$$

    The symmetry of the Hooke matrix $\mathbb{H}$ implies $H_{ij} = H_{ji}$. In text books about composite and anisotropic materials ^13,14, the reader can find not only the expressions linking the stiffness terms of expression (8) to the engineering constants ($E_{i}$, $G_{ij}$ and $\nu_{ij}$) but also the stability criteria which these stiffness terms must fulfill in order to represent a material with positive deformation energy. The user has to specify two lines for each temperature²⁶:

    $$\begin{align*} &H_{11}, H_{12}, H_{22}, H_{13}, H_{23}, H_{33}, H_{44}, H_{55},\\ &H_{66}, T \end{align*}$$

    The orthotropic elastic properties of the cementite phase of X65 pipeline steel are given as Ex. 8.

    /ELASTIC, TYPE=ORTHO
    388000., 156000., 345000., 164000., 162000., 322000., 15000., 134000.,
    134000., 23.
    #
    # Reference: Jiang et al. Jnl. of applied Physics, vol. 103, 043502, 2008.
    #
    /EXPANSION, TYPE=ORTHO
    0.000006, 0.000006, 0.000006, 23.
    0.000007, 0.000007, 0.000007, 200.
    ........
    0.0000190, 0.0000190, 0.000190, 1000.
    #
    

    Example 8: Definition of the orthotropic elastic properties of the cementite phase.

  • TYPE=CUBIC: cubic symmetric elasticity
    Metal phases, like the bcc ferrite or the fcc austenite in steels, present often with cubic symmetry, which is a special subclass of orthotropic materials. In this case, the diagonal stiffness terms $H_{ii}$ $(i=1,3)$ are identical; the shear stiffness terms $H_{ii}$, $i=4,6$ are also identical as well as the off-diagonal terms $H_{ij}$, $i\neq j$, $i,j=1,3$. In fact, the material is defined only by 3 independent elastic constants: $H_{11}$, $H_{12}$ and $H_{44}$. The user has only to specify:

    $$\begin{equation*} H_{11}, H_{12}, H_{44}, T \end{equation*}$$

    /ELASTIC, TYPE=CUBIC
    233100., 135440., 117830., 27.
    228425., 131744., 115467., 100.
    ..............
    127525., 73550., 64463., 900.
    #
    # Ref.: Ghosh & Olson, Acta Mat. 50 (2002) 2655-75 plus experimental data from
    # Rayne & Chandrasekhar, Phys. rev, 122 (1961) 1714-16.
    #
    

    Example 9: Definition of the elastic coefficients of the ferrite phase of a FeCMnSi alloy.
    - TYPE=ANISO: fully anisotropic elasticity
    For a fully anisotropic elastic material, 21 independent elastic stiffness parameters have to be specified. The symmetric Hooke matrix, in Voigt notation, of an anisotropic material is given by:
    $$\begin{equation} \label{eq:H_aniso} \mathbb{H} = \begin{bmatrix} H_{11} & H_{12} & H_{13} & H_{14} & H_{15} & H_{16} \\ & H_{22} & H_{23} & H_{24} & H_{25} & H_{26} \\ & & H_{33} & H_{34} & H_{35} & H_{36} \\ & S & & H_{44} & H_{45} & H_{46} \\ & & Y & & H_{55} & H_{56} \\ & & & M & & H_{66} \\ \end{bmatrix} \end{equation}$$

    In order to define this Hooke matrix, the user has to specify the 21 stiffness terms for each temperature²⁶ in following manner:

    $$\begin{align*} &H_{11}, H_{12}, H_{22}, H_{13}, H_{23}, H_{33}, H_{14}, H_{24},\\ &H_{34}, H_{44}, H_{15}, H_{25}, H_{35}, H_{45}, H_{55}, H_{16},\\ &H_{26}, H_{36}, H_{46}, H_{56}, H_{66}, T \end{align*}$$

    Note: HOMAT makes internally a further differentiation: it distinguishes between fully anisotropic material behavior with 21 independent Hooke matrix terms and the monoclinic case, with a symmetry plane either parallel to the $Z=0$ plane or to the X or Y planes. In these special cases, only 13 independent material coefficients have to be specified. In the thermo-elastic case, the orthotropic thermal expansion coefficients $\alpha_{i}\,,\text{ for } i=1, 3$ must be specified. A full anisotropic thermal expansion tensor $\alpha_{ij}$ cannot be introduced at the moment. As example, the monoclinic elastic behaviour of a mono-crystalline sapphire fibre is outlined. The present isotropic Young modulus and Poisson coefficient are ignored when the parameter TYPE of the /MATERIAL option is set to ANISO; whereas the /ELASTIC keyword is ignored if TYPE is set to ISO.

    /YOUNG MODULUS, TYPE=ISO
    414000.0, 20.0
    407000.0, 300.0
    -----------
    390000.0 900.0
    385000.0, 1200.0
    #
    # unit 1 N/mm2 = 1 Mpa
    #
    /POISSON RATIO, TYPE=ISO
    0.27, 0.0
    #
    /ELASTIC, TYPE=ANISO
    456900., 149400., 456900., 106100., 106100., 462200., 0.0, 0.0,
    0.0, 153700., 0.0, 0.0, 0.0, 20700., 135400., 20700.0,
    -20700., 0.0, 0.0, 0.0, 135400., 1100.
    #
    # thermal expansion coefficient 1/(Grad Celsius) function of temperature
    /EXPANSION, TYPE=ORTHO
    0.0000069, 0.0000033, 0.0000033, 100.0
    0.0000083, 0.0000077, 0.0000077, 500.0
    0.0000090, 0.0000083, 0.0000083, 1000.0
    0.0000095, 0.0000088, 0.0000088, 1500.0
    #
    

    Example 10: Definition of the monoclinic elastic behaviour of the mono-crystalline sapphire fibre.

• /EXPANSION, TYPE = type

  This keyword defines the thermal expansion coefficients of a phase or constituent material of the RVE. They are interpreted as total expansion coefficients with respect to a reference temperature $T_{ref}$, i.e. the thermal strain $\epsilon_{th}$ of a material at final temperature²⁶ $T$ is given by

  $$\begin{equation} \label{eq:th_expasion} \varepsilon_{\text{th}} = \mathbf{\alpha}(T)(T - T_{\text{ref}})\,, \end{equation}$$

  where $\alpha(T)$ is the thermal expansion tensor at temperature $T$. As the objective of HOMAT is to determine the effective thermal expansion of the heterogeneous material, the temperature difference is always set to 1°C without loss of generality. Therefore, we have: $\varepsilon_{\text{th}} = \alpha(T)$ The optional parameter TYPE on the keyword line specifies the adopted material model: ISO, TR_ISO or ORTHO. By default, TYPE=ISO. TYPE=TR_ISO has to be specifies for the transverse isotropic elastic behaviour, whereas TYPE=ORTHO has to be defined for an orthotropic, cubic symmetric, fully anisotropic or monoclinic elastic material behaviour. Depending on the TYPE value, the user specifies for each temperature²⁶:

  • TYPE=ISO
    

    $$\begin{equation*} \alpha, T \end{equation*}$$

  • TYPE =TR_ISO
    

    $$\begin{equation*} \alpha_{\text{t}}, \alpha_{\text{p}}, T \end{equation*}$$

  • TYPE=ORTHO
    

    $$\begin{equation*} \alpha_{11}, \alpha_{22}, \alpha_{33}, T \end{equation*}$$

  where $\alpha_{11}$, $\alpha_{22}$ and $\alpha_{33}$ are the thermal expansion in the directions 1, 2 or 3 of the local material axis system of the material domain (specified by the /ORIENTATION option) respectively; for the transverse isotropic case, $\alpha_{\text{t}}$ is the thermal expansion coefficient in the transverse direction and $\alpha_{\text{p}}$ the thermal expansion coefficient in the isotropic plane. The unit of the thermal expansion coeffcient is usually °C^-1. Several definition examples have been already provided (see Ex. 5-8 and 10).

  Note: Alternatively to /EXPANSION the user can introduces the keyword /THERMAL EXPANSION. Only the diagonal terms of the thermal expansion tensor can be introduced at the moment. Off-diagonal terms are neglected.
  

• /MOLAR VOLUME

  This keyword defines the volume occupied by a mole of the considered phase/material. As mm is the default geometry dimension, its unit is mm³/mol. In presence of different phases, this physical quantity must be specified for each phase in order to evaluate the corresponding isotropic, volumetric eigenstrain. For example, it occurs during a solid or liquid phase transformation, like the austenite to ferrite one in heat treatment processes. This quantity is directly linked to the phase/material density and to the molar weight and is, of course, temperature²⁶ dependent. No additional parameter has to be specified. The user enters simply:

  $$\begin{equation*} V_{mol}, T \end{equation*}$$

  The temperature dependent molar volume of the ferrite phase is given in Ex. 11.

  #
  /MOLAR VOLUME
  7099.52, 20.
  7121.51, 100.
  ............
  7368.35, 850.
  7391.90, 911.
  # expressed via the density and the molar weight of the alloy: 55.8309
  # in mm^3/mol
  #
  

  Example 11: Definition of the temperature dependent molar volume of the ferrite phase of a FeCMnSi alloy.

• /POISSON RATIO, TYPE = type

  This keyword is used in the case of isotropic, transverse isotropic or orthotropic material behaviour, defined by engineering constants, to introduce the Poisson coefficients of the considered phase or constituent material. Again in accord with the specified material model on the /MATERIAL keyword, the parameter TYPE have to be defined on the keyword line. Classically, the user specifies:

  • TYPE=ISO

    $$\begin{equation*} \nu, T \end{equation*}$$

  • TYPE=TR_ISO

    $$\begin{equation*} \nu_{\text{tp}}, \nu_{\text{p}}, T \end{equation*}$$

  • TYPE=ORTHO

    $$\begin{equation*} \nu_{12}, \nu_{13}, \nu_{23}, T \end{equation*}$$

  It is important to note that, for a correct isotropic material behaviour definition, the option /POISSON RATIO must be accompanied by the keyword /YOUNG MODULUS, whereas in the case of transverse isotropic or orthotropic material behaviour, this option must be specified in conjunction with the keywords /SHEAR MODULUS and /YOUNG MODULUS. Note that for the isotropic transverse material behaviour, the definition of the transverse and in-plane Poisson coefficients, $\nu_{\text{tp}}$ and $\nu_{\text{p}}$ respectively, are identical as for the keyword /ELASTIC. Ex. 10 illustrates the definition of an “isotropic” sapphire fibre. The orthotropic case is outlined in Ex. 12 for the yttrium stabilized zirconia oxide used as thermal barrier coating.

  # Young modulus (T) MPa
  #
  /YOUNG MODULUS, TYPE=ORTHO
  5.096191E+04, 3.124513E+04, 6.37023875E+04, 700.
  5.077967E+04, 3.121674E+04, 6.34745875E+04, 750.
  .........................
  5.011533E+04, 3.099163E+04, 6.26441625E+04, 900.
  4.975417E+04, 3.087421E+04, 6.21927125e+04, 950.
  #
  /POISSON RATIO, TYPE=ORTHO
  2.705990E-01, 3.442350E-01, 3.434799E-01, 700
  2.698298E-01, 3.424140E-01, 3.416010E-01, 750
  ..........................
  2.682793E-01, 3.369597E-01, 3.351691E-01, 900
  2.673659E-01, 3.347734E-01, 3.331901E-01, 950
  #
  # shear modulus (T) MPa
  #
  /SHEAR MODULUS, TYPE=ORTHO
  1.386069E+04, 1.386069E+04, 1.386069E+04, 700.
  1.374312E+04, 1.374312E+04, 1.374312E+04, 750.
  .......................
  1.351601E+04, 1.351601E+04, 1.351601E+04, 900.
  1.349206E+04, 1.349206E+04, 1.349206E+04, 950.
  #
  

  Example 12: Specification of the orthotropic material behaviour of the yttrium stabilized zirconia oxide by the Young, shear modules and the Poisson ratio's. This definition is equivalent to the specification via the /ELASTIC keyword (see Example 7).

• /SHEAR MODULUS, TYPE = type

  The user needs to specify this option only in the case of transverse isotropic or orthotropic material behavior, defined by a combination of Young and shear moduli and the Poisson coefficients. The TYPE parameter equals either TR_ISO or ORTHO. The shear modules are usually specified in MPa and can depend on the temperature²⁶. The user introduces these modules as follows:

  • TYPE=TR_ISO

    $$\begin{equation*} G_{\text{t}}, 0.0, T \end{equation*}$$

    Note that a second value must be specified for technical reasons. It is best set to 0.0.

  • TYPE=ORTHO

    $$\begin{equation*} G_{12}, G_{13}, G_{23}, f_{S} \end{equation*}$$

    Note: The definition of the shear planes, noted here 1-2, 1-3 and 2-3, are achieved via the /ORIENRATION option. This option must be specified for each orthotropic or anisotropic material/phase.

  Ex. 12 illustrates also the definition of this keyword.

• /YOUNG MODULUS, TYPE = type

  The Young modulus describes the basic engineering concept of elasticity. The keyword /YOUNG MODULUS can be used to describe either an isotropic, transverse isotropic or orthotropic material behaviour, defined by engineering constants. The additional parameter TYPE defines the nature of the material model. Three values are permitted: TYPE=ISO (default), TYPE=TR_ISO and TYPE=ORTHO. Only if its value corresponds to the value specified in the /MATERIAL option, then it will be used to specify the Hooke matrix of the material, otherwise it will be ignored. Note that this option has to be used in combination with the keyword /POISSON RATIO for an isotropic material and with /POISSON RATIO and /SHEAR MODULUS for a transverse isotropic or orthotropic material behaviour. The Young moduli are usually specified in MPa and can depend on the temperature²⁶. The user introduces either:

  • TYPE=ISO

    $$\begin{equation*} E, T \end{equation*}$$

  • TYPE=TR_ISO

    $$\begin{equation*} E_{\text{t}}, E_{\text{p}}, T \end{equation*}$$

  • TYPE=ORTHO

    $$\begin{equation*} E_{1}, E_{2}, E_{3}, T \end{equation*}$$

    where 1,2,3 specify the directions of orthotropic symmetry; $E_{\text{t}}$ defines the transverse Young modulus and $E_{\text{p}}$ the in-plane Young modulus (see also /ELASTIC, TYPE=TR_ISO).

  Example 8 illustrates the definition of this keyword.

Definition of the Microstructure Model¶

In this chapter all the keywords describing the geometrical model of the RVE on the micro-scale are outlined in an alphabetical order. The geometrical model has been pre-treated by Mesh2HOMAT in order to generate the needed interface elements at phase/material boundaries and the periodic boundary conditions. All the geometrical information about the microstructure is summarized in the PROBLEM_NAME.hgeo file, which is read by HOMAT. This file can be included easily in the command file of a specific homogenization application via the /INCLUDE command explained more in detail in Sec. Definition of HOMAT- specific Commands. This inclusion happens directly after the header, the problem, output and unit definitions, as shown in Ex. 13. After the description of the geometrical information provided by Mesh2HOMAT, the specific HOMAT keywords related to the microstructure model, like the /SOLID SECTION, /INTERFACE or /PERIODICITY options, are outlined in detail.

/HEADING
Date: 25.04. Time: 10:24:20
Total number of Elements: 296
Total number of Nodes: 732
Total number of new Nodes: 76
Total number of Materials: 2
#
/PROBLEM, TYPE=THELAS, HOM_MOD=0, PER_MOD=2
#
/OUTPUT, FORMAT=VTK, TYPE=ASCII, IOP=2, VOLAVRG=YES, ENG=NO, VISLAV=YES, LOC_AXIS=NO
U
#
/UNITS, GEOMETRY=1.0, LENGTH UNIT=MM, TEMP OFFSET=273.15, ANGLE=GRAD
#
/INCLUDE, NAME=ehz3D_iso.hgeo.gz
#
/MATERIAL, NAME=VERRE, TYPE=ISO
/INCLUDE, NAME=verre.hmat
#
/MATERIAL, NAME=EPOXY
/INCLUDE, NAME=epoxy.hmat
/SOLID SECTION, NUMBER=1, ELSET=FIBER, NAME=fibre, MATERIAL=VERRE, TEMP=NO, MECH=CALC, FLOW=NO
#
/SOLID SECTION, NUMBER=2, ELSET=RESIN, NAME=matrice, MATERIAL=EPOXY,
TEMP=NO, MECH=HOM, FLOW=NO
#
/INTERFACE, NAME=FIBER:RESIN, NUMBER=1,
SURFACE1=INTERFACE_FIBER:RESIN,
SURFACE2=INTERFACE_RESIN:FIBER
#
/PERIODICITY, NAME=X-X+, TYPE=TRANS, NUMBER=1, START=X-, TARGET=X+
3, 0.0, 0.0
#
/PERIODICITY, NAME=X+X-, TYPE=TRANS, NUMBER=2, START=X+, TARGET=X-
-3, 0.0, 0.0
#
/PERIODICITY, NAME=Y-Y+, TYPE=TRANS, NUMBER=3, START=Y-, TARGET=Y+
0.0, 3, 0.0
#
/PERIODICITY, NAME=Y+Y-, TYPE=TRANS, NUMBER=4, START=Y+, TARGET=Y-
0.0, -3, 0.0

Geometrical model information provided by Mesh2HOMAT¶

• /ELEMENT, TYPE = type, ELSET = elset

  This keyword defines the volume finite elements (FE) in a specific RVE domain of the heterogeneous material. Two additional parameters must be specified: TYPE and ELSET. The parameter TYPE specifies the kind of finite element used for the considered domain. The element geometries available in HOMAT are summarized in Tab. 1.

  Table 1: Element types available in HOMAT.

  type	description
  C3D4	4-node linear tetrahedral element
  C3D6	6-node linear pentahedral element
  C3D8	8-node linear hexahedral element
  C3D10	10-node quadratic tetrahedral element
  C3D15	15-node quadratic pentahedral element
  C3D20	20-node quadratic hexahedral element

  The user familiar with ABAQUS ⁶ has directly recognized that the same element nomenclature is used in HOMAT as in ABAQUS.

  Note: The use of quadratic elements is recommended to special applications where the accuracy of the prediction of effective properties is mandatory or in the localization analysis, where accurate local stresses near material interfaces are needed. Their use will severely increase of the number of D.O.F and, subsequently, of the CPU time and the stored data. In order to model accurately the mechanical and fluid flow behaviour of quasi-incompressible materials and fluids, like melt in solidification processes, amorphous phase of semi-crystalline polymers, mixed linear volume finite elements have been implemented in HOMAT 6.0. For each volume element type (tetra-, penta- and hexahedron), following different mixed variants have been introduced, detailed here for tetrahedral elements:

  • C3D4E 4-node equal order tetrahedral element with linear displacement and pressure fields. This element version is not stabilized. Therefore, its use is not recommended and restricted to comparison purpose in scientific and benchmark applications.
  • C3D4L 4-node mixed stabilized equal order tetrahedral element with bilinear displacement and pressure fields. Following Li & He ⁹, this element version is stabilized by a double Gauss point integration of the pressure term, avoiding pressure oscillations.
  • C3D4M 4-node tetrahedral element where the linear displacement field is enriched by a 3^d order volume bubble term. Its pressure field is always linear. As Arnold et al. ¹⁵ has proven, this tetrahedral element, named mini element, satisfies the inf-sup conditions and is automatically stable. However, the extension of the mini element to hexahedral C3D8M or pentahedral C3D6M elements is not straightforward, as, for example for the hexahedron, three internal nodes instead one for the tetrahedron has to be added in order to satisfy the inf-sup condition.
  • C3D4B 4-node mixed stabilized tetrahedral element with bilinear displacement field enriched by a volume bubble term and a bilinear pressure field. This “bubble” element version is further stabilized by the Li-He method, by applying a double Gauss point integration of the pressure term.

  For elements with pentahedral geometry, we can specify following mixed element variants: C3D6E, C3D6L, C3D6M and C3DB. For the hexahedrons, we can define the same variants, namely: C3D8E, C3D8L, C3D8M and C3D8B. Note that for penta- and hexahedral geometries, only the LiHe and the bubble variants are stable. In Mesh2HOMAT the new option discre allows the user to specify the desired finite element discretization for each domain/phase. By default, the classical bilinear displacement elements C3D4, C3D6 and C3D8 are used. The numerical benchmark tests (homogenization of a silicate sphere in an epoxy matrix, permeability prediction of flow around a cylinder, ..) have shown that the stabilized bubLi variant (C3D4B, C3D6B and C3D8B) with enriched displacement field, is the most accurate element version. It is recommended to use this element version in mechanical applications with quasi-incompressible materials and in Stokes flow homogenizations.

  The parameter ELSET specifies the name of the group of elements, which will be introduced under the same keyword /ELEMENT. They define a specific element set, which can be easily addressed by its unique name. In Ex. 14 it defines, for example, the glass fibre. The keyword /ELEMENT can be repeated several times to describe for example a complex domain composed of different types of elements. Independent of the element type, an element of the considered set is defined by its identification number ID_ele followed by the list of the element nodes. For a linear hexahedral element, we have:

  ID_el, N1, N2, N3, N4, N5, N6, N7, N8
  

  In Ex. 14 a selection of the geometry file generated by Mesh2HOMAT of the simple RVE of a reinforced UD glass-epoxy composite (see Fig. 4) is presented. The /ELEMENT keyword has been used twice: once to define the fibre and once to specify the matrix around the fibre.

  

  Figure 4: Finite element mesh of the 3-D RVE of the UD reinforced glass-epoxy composite.

  #
  /NODE, NSET=ALLNODES
  1, -1.190000e+00 ,+0.000000e+00 ,+2.500000e+01
  2, +1.190000e+00 ,+0.000000e+00 ,+2.500000e+01
  3, -1.173767e+00 ,+1.958856e-01 ,+2.500000e+01
  ...................
  731, +1.125508e+00 ,-3.864357e-01 ,+2.500000e+01
  732, +1.125523e+00 ,-3.863900e-01 ,-2.500000e+01
  /ELEMENT, TYPE=C3D8, ELSET=FIBER
  1, 39, 21, 1, 3, 131, 124, 112, 111
  2, 20, 2, 32, 73, 94, 93, 113, 165
  3, 49, 39, 3, 4, 141, 131, 111, 110
  ...........................................
  223, 421, 424, 423, 420, 499, 502, 501, 498
  /ELEMENT, TYPE=C3D8, ELSET=RESIN
  76, 209, 254, 189, 185, 267, 335, 289, 266
  77, 249, 218, 186, 195, 330, 299, 282, 283
  ............................................
  295, 586, 582, 583, 587, 647, 643, 644, 648
  296, 594, 592, 593, 595, 655, 653, 654, 656
  #
  /ELSET, NAME=ALLELEMENTS
  FIBER
  RESIN
  #
  /SURFACE, NAME=INTERFACE_FIBRE:RESIN
  1 , 6
  149 , 6
  2 , 5
  150 , 5
  3 , 6
  .........................
  163 , 3
  164 , 5
  165 , 5
  166 , 5
  167 , 5
  #
  /SURFACE, NAME=INTERFACE_RESIN:FIBRE
  81 , 5
  229 , 5
  80 , 6
  ......................
  269 , 4
  270 , 6
  272 , 6
  273 , 6
  274 , 6
  /SURFACE, NAME=X-
  78, 5
  79, 5
  .................
  253, 5
  254, 5
  /SURFACE, NAME=X+
  76, 6
  77, 6
  ................
  234, 6
  235, 6
  /SURFACE, NAME=Y-
  76, 4
  78, 4
  ...............
  99, 4
  100, 4
  /SURFACE, NAME=Y+
  224, 3
  226, 3
  ...............
  247, 3
  248, 3
  

  Example 14: Geometrical model file UD_compos.hgeo of the UD reinforced glass-epoxy composite generated by Mesh2HOMAT.

• /ELSET, NAME = name, [ GENERATE ]

  This option is used to assign elements to a specific element set or to combine existing element sets to a new one. The parameter NAME must be specified on the keyword line. In the Ex. 14, the NAME is ALLELEMENTS. An element set with the same name must not already exist. The set name is case insensitive. The user has three ways to specify an ELSET:

  1. List of element identification numbers, maximum 16 entries per line:

     /ELSET, NAME=CRIT_AREA
     45, 56, 89, 95, 85, 206, 218
     

  2. In order to facilitate the explicit definition of the element numbers, belonging to a specific set, an optional parameter GENERATE can be specified on the keyword line. If present, element ranges can be expressed by their initial element number, their final number and an increment.

     /ELSET, NAME=E2, GENERATE
     20, 60, 5
     

     In this case, the elements 20, 25, 30, 35, 40, 45, 50, 55 and 60 have been assigned to element set E2.

  3. As shown in Ex. 14, an ELSET may be specified as the union of previously defined element sets. Each subset appears on a separate line. In the special case of Ex. 14, the new element set ALLELEMENTS is the union of the element sets FIBER and MATRIX.

• /NODE, [ SCAL ]

  This keyword specifies the nodes with their 3-D coordinates of a part of the RVE or of the whole RVE of a heterogeneous material. The coordinates are given in the structural axis system of the RVE. At the moment, no local axis system can be defined for the introduction of nodal coordinates. On the keyword line the optional parameter SCAL can be introduced. This parameter allows the multiplication of the coordinates by the value of this parameter. Note that this scale factor overrides the scale value introduced via parameter GEOMETRY of the more general keyword /UNITS (see Sec. Options introduced before the RVE description). The users specifies for example:

  /NODE, SCAL=0.001
  1, 0., 0., .0.
  2, 10., 0., 0.
  3, 0., 10., 0.
  

  The coordinates of nodes 1-3 are multiply by the factor 10^-3 in order to pass from a micrometre to a millimetre definition of the RVE geometry.

• /NSET, NAME = name, [ GENERATE ]

  Analogous to the definition of the element set /ELSET, keyword /NSET is used to assign nodes to a specific node set. The parameter NAME, specifying the name of the considered node set, is required. A node set with the same name must not already exist. The set name is also case insensitive. The user has three ways to specify a specific NSET:

  1. List of node numbers, maximum 16 entries per line

     /NSET, NAME=AROUND_HOLE
        45, 56, 89, 95, 85, 206, 218
     

  2. In order to facilitate the generation of the list of node numbers an optional parameter GENERATE can be specified on the keyword line. If present, nodes ranges can be expressed by their initial node number, their final number and an increment.

     /NSET, NAME=FREE_SURF, GENERATE
     12, 24, 2
     

     In this case, the nodes 12, 14, 16, 18, 20, 22 and 24 have been assigned to the node set FREE_SURF.

  3. List of previously defined nodes sets - one per line - to specify a new node set. This way it is easy to join already defined node sets.

• /SURFACE, NAME = name

  This keyword is used to assign element faces to a specific surface group. The parameter NAME, specifying the name of the considered surface, is required. Surfaces are used to specify the interfaces between different phases, materials, where jumps of the thermo-physical properties occur. They are also used to define the periodic boundary conditions at the free surfaces of the RVE and contact heat transfer between phases, constituents inside the RVE for thermal homogenization. The numbering of the faces, belonging to the surface, depends on the element type (see Figs. 5 and 6 and Tabs. 2 - 4).

  

  Figure 5: Linear hexahedron (left), linear pentahedron (center), linear tetrahedron (right).

  

  Figure 6: Quadratic tetrahedron (left), quadratic pentahedron (center), quadratic hexahedron (right).

  Table 2: Faces and nodes of tetrahedral elements.

  face	nodes (linear)	nodes (quadratic)
  1	1, 2, 3	1, 2, 3, 5, 6, 7
  2	1, 4, 2	1, 4, 2, 8, 9, 5
  3	2, 4, 3	2, 4, 3, 9, 10, 6
  4	3, 4, 1	3, 4, 1, 10, 8, 7

  Table 3: Faces and nodes of hexahedral elements.

  face	nodes (linear)	nodes (quadratic)
  1	1, 2, 3, 4	1, 2, 3, 4, 9, 10, 11, 12
  2	5, 8, 7, 6	5, 8, 7, 6, 18, 15, 14, 12
  3	1, 5, 6, 2	1, 5, 6, 2, 17, 12, 18, 9
  4	2, 6, 7, 3	2, 6, 7, 3, 18, 14, 19, 10
  5	3, 7, 8, 4	3, 7, 8, 4, 19, 15, 20, 11
  6	4, 8, 5, 1	1, 5, 8, 4, 17, 16, 20, 12

  Table 4: Faces and nodes of pentahedral elements.

  face	nodes (linear)	nodes (quadratic)
  1	1, 2, 3	1, 2, 3, 7, 8, 9
  2	4, 6, 5	4, 6, 5, 12, 11, 10
  3	1, 4, 5, 2	1, 4, 5, 2, 13, 10, 14, 7
  4	2, 5, 6, 3	2, 5, 6, 3, 14, 11, 15, 8
  5	1, 3, 6, 4	1, 3, 6, 4, 9, 15, 12, 13

  Independent of the element type, the faces are numbered in such a way that the normal to this face is always oriented inside the volume of the element. This face definition is completely identical to the notation used in ABAQUS. The user specifies on the keyword line the parameter NAME and its value. Then, on the following line the element number and the number of face is specified. This line is repeated until all faces belonging to the considered surface are introduced.

  /SURFACE, NAME=FreeSurf_1
  1, 6
  2, 1
  

  This entry assigns the nodes of face 6 of element 1 and the nodes belonging to face 1 of element 2 to the surface with the name FreeSurf_1. For the RVE of the UD reinforced composite of Fig. 4, the surfaces X-, X+, Y- and Y+ are specified on the free surfaces of the RVE, which are normal to the X or Y direction respectively. These surface definitions are required by the /PERIODICITY keyword (see Sec. Model information defined in the command file). At the interface between the fibre and the matrix two surfaces named INTEFACE_FIBRE:RESIN and INTERFACE_RESIN:FIBRE respectively, are automatically defined by Mesh2HOMAT. The nomenclature of the surface names of an interface between two grains, phases, or material regions follows always following rules:

  • 1^st surface: INTERFACE_{1st material region}:{2nd material region}
  • 2^nd surface: INTERFACE_{2nd material region}:{1st material region}

  where the definition of the 1^st and 2^nd material regions depends of their order in the list of ELSET's.

• /XTUPLE, *=MAX ==*/=nbr=/, *=NUM ==*/=nodes=/

  This new keyword provides for each geometrical master node, where several (> 2) materials are connected, the list of nodes with the same coordinates. The mandatory second parameter MAX specifies the maximum number of nodes having the same coordinates; the second parameter NUM specifies the number of XTUPLE specified automatically by Mesh2HOMAT.

  For the hollow_grid example, Mesh2HOMAT writes in the grid_mono_th.hgeo file following /XTUPLE information:

  #
  /XTUPLE, MAX=3, NUM=648
  5, 32154, 32156
  7, 32174, 32176
  17, 32242, 32245
  24, 32292, 32295
  .....
  2176, 34294, 34295
  2177, 34211, 34218
  2178, 34258, 34262
  

  Example 15: Part of geometrical model file grid_mono_th.hgeo of the monolayer of hollow wire frame generated by Mesh2HOMAT.

Model information defined in the command file¶

• /HEAT TRANSFER, TYPE = IDEAL | FCT, [ HVAL = hval ]

  This keyword is only read in thermal homogenization simulations. It is a sub-keyword related to a specific interface and must be introduced directly behind the concerned /INTERFACE. As millimeter is the default unit of length, the heat transfer coefficient is usually specified in W/(mm² K). Dependent on the nature of the two regions in contact, the heat transfer coefficient can be assumed to be ideal or has a temperature dependent value, obtained from experiments. The user must define on the keyword line the nature of the heat transfer coefficient: either IDEAL or FCT (= function of the temperature).

  • Ideal heat transfer

    /HEAT TRANSFER, TYPE=IDEAL, HVAL=0.008
    

    In the case of ideal heat transfer between two contact regions, the user can defined a specific value corresponding to this ideal heat transfer. By default, HVAL=0.01 W/(mm² K). In the present example, it is defined a little lower: 0.008 W/(mm²K). If the parameter HVAL is specified, its value overrides the parameter HTR_IDEAL, which is defined in the /HOM PARA keyword (see Sec. Definition of HOMAT- specific Commands) for the whole RVE.

  • Function of the temperature

    /HEAT TRANSFER, TYPE=FCT
       0.0005, 500.0
       0.0008, 700.0
       0.001. 900.0
    

    Linear interpolation with the temperature is achieved between these values. If the temperatures of both regions in contact differ, HOMAT takes, per convention, the temperature of the 1^st region in the region list as the reference temperature at which the interpolation happens.

• /INTERFACE, NAME = name , NUMBER = number , SURFACE1 = surface1 , SURFACE2 = surface2

  This keyword is automatically generated by Mesh2HOMAT and allows the definition of an interface between two grains, phases or materials at the micro-scale. Both surfaces belonging to an interface with zero thickness and having double coincident nodes (one per surface in contact) must have a coherent surface mesh: the number of nodes and triangular or quadrangular surface elements must be exactly the same on both interface sides. The keyword line can be followed by a /HEAT TRANSFER option in the case of thermal homogenization. For the unique interface between fibre and matrix of Ex. 13, Mesh2HOMAT generates the following interface definition:

  /INTERFACE, NAME=FIBER:RESIN, NUMBER=1,
  SURFACE1=INTERFACE_FIBER:RESIN,
  SURFACE2=INTERFACE_RESIN:FIBER
  

  On the keyword line, the parameter NAME is required and specifies the name of the interface. The name is always a combination of the names of the two domains in contact, here FIBER and RESIN. Further a second parameter NUMBER is introduced there. It specifies the identification number of the interface. In this simple example, this number equals 1. On the second and third line, the names of both surface groups, previously defined in the PROBLEM_NAME.hgeo file, are introduced respectively. Parameter SURFACE1 specifies always the surface of the first specified region in the geometrical model description and SURFACE2 the 2^nd one. In Ex. 13 the fiber domain is the 1^st region and the resin the 2^nd region. Mesh2HOMAT generates such a specific interface for each combination of different regions, having their own properties. In the case of grains of a polycrystal, belonging to the same phase, e.g. ferrite or austenite, the material properties are the same, of course, but not their orientation. In this case also a specific interface is generated between adjacent grains. For complex polycrystal microstructures (> 100 grains) several interfaces (> 1000) have to be specified. In this case manual definition is hardly an alternative to automatic generation via Mesh2HOMAT.

• /ORIENTATION, NAME = name, DEFINITION = definition, [ ROTX = rotx ], [ OUT_ORIE ]

  This option is used to define a local coordinate system for a region, grain in order to specify its transverse isotropic, orthotropic or anisotropic material properties or in order to specify the local axis system in which the calculated effective properties are output. On the keyword line the required parameter NAME assigns a name to each orientation definition. This name is referenced in the /SOLID SECTION keyword (see below). The second parameter, DEFINITION, specifies the nature of the introduced orientation definition. At the moment, HOMAT recognizes following orientation definitions:

  • DEFINITION=QUATERNION
    A quaternion has four components and is composed of a vector plus a scalar: $\mathbf{q} \triangleq \lbrack\mathbf{v}^{T},\ q_{4}\rbrack$, whose norm is one: $\left| v \right|^{2} + \ q_{4}^{2} = \ \mathbf{q}^{T}\mathbf{q} = 1$. A quaternion representing without the rotation by an angle α around a unit vector $\vec{a}$ in 3-D has the following expression

    $$\begin{equation} \label{eq:quaternion} [a_{1}\sin(\alpha/2), a_{2}\sin(\alpha/2), a_{3}\sin(\alpha/2), \cos(\alpha/2)]\,. \end{equation}$$

    The user introduces for example:

    /ORIENTATION, NAME=grain_12, DEFINITION=QUATERNION
    0.880735, 0.2870, -0.33212, 0.177857
    

  • DEFINITION=EULER
    In 3-D any rotation of a structural axis system to a local axis system and back, can be described by a sequence of three Euler rotations about coordinate axes, where two successive rotations may not be done around the same axis. Several ways to define the Euler angles exist in the literature. HOMAT uses the convention shown in Fig. 7, where the first rotation is about the Z structural axis (angle $\phi$, rotation matrix $\vec{D}$); followed by a rotation about the new X' axis (angle $\theta \in [0, \pi]$, rotation matrix $\vec{C}$) and finally a rotation about Z' is achieved (angle $\psi$, rotation matrix $\vec{B}$). Any kind of rotation matrix $\vec{A}$ from a structural to a local axis system can then be expressed in 3-D by following successive matrix products:

    $$\begin{equation} \label{eq:euler_ang} \vec{A} = \vec{B} \vec{C} \vec{D} \end{equation}$$

    

    Figure 7: Definition of the Euler angles $\phi$, $\theta$, $\psi$ for a transformation from the structural to local axis system.

    For an Euler angle orientation definition in degrees, the user introduces following sequence:

    /ORIENTATION, NAME=grain_8, DEFINITION=EULER
    90, 14.0, 0.
    

    where $\phi=90$°, $\theta=14$° and $\psi=0$°.

  • DEFINITION=ANGLE_AXIS
    An orientation in the 3-D space can also be defined by a unit axis, defined by its components in the structural system and by a rotation angle $\alpha$. In this case the user provides simply at first the angle in degrees followed by the three components of the rotation axis in the 3-D space.

    /ORIENTATION, NAME=grain_8, DEFINITION=ANGLE_AXIS
    90.0, 0.0, 1.0, 0.0
    

  • DEFINITION=ROTMA
    The user can also directly specify the unitary rotation matrix $\vec{A}$ (eq. 12), which is composed of the scalar products of the unit vectors of structural ($\vec{E}_{j}$) and the new local $\vec{e}_{i}$ axis system:

    $$\begin{equation*} a_{ij} = \vec{E}_{i} \cdot \vec{e}_{j} \end{equation*}$$

    The explicit expression of the rotation matrix $\vec{A}$ is given by

    $$\begin{equation} \label{eq:rot_mat} \vec{A} = \left( \begin{matrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{matrix} \right)\,. \end{equation}$$

    The definition of a local orientation axis system via the rotation matrix $\vec{A}$ is achieved by introducing its components. This way, a direct link to the output of EBSD and metallographic measurements is achieved. Note that, HOMAT performs the inverse transformation from the local axis system, in which the material properties are defined, to the structural one, in which the homogenization analysis is performed.

    /ORIENTATION, NAME=grain_18, DEFINITION=ROTMA
    a11, a12, a13, a21, a22, a23, a31, a32, a33
    

    or , for a particular example,

    /ORIENTATION, NAME=grain_18, DEFINITION=ROTMA
    0.29237, 0.0, -0.9563, 0.0, 1.0, 0.0, 0.9563, 0.0, 0.29237
    

  • DEFINITION=COORDINATES
    Analogous to the orientation definition in ABAQUS, the user can also specify a local, rectangular coordinate axis system by defining the origin of this system $\vec{c}$ and two points $\vec{a}$ and $\vec{b}$, as shown in Fig. 8 Point $\vec{a}$ must be on the local axis X' and point $\vec{b}$ must be in the plane X'-Y', not necessary located on axis Y'. Intuitively, the user selects point $\vec{b}$ on or near the local axis Y' (see Fig. 8). If definition COORDINATES is specified, points $\vec{a}$ and $\vec{b}$ are defined directly by introducing their coordinates. If the user wants to define an origin of the local axis different from the structural one, he must also specify the coordinates of a third point $\vec{c}$. Otherwise, the coordinates of the two nodes $\vec{a}$ and $\vec{b}$ suffice. In the general case with different origin, the user introduces following orientation definition sequence:

    /ORIENTATION, NAME=some_name, DEFINITION=COORDINATES
    a1, a2, a3, b1, b2, b3, c1, c2, c3
    

    or , for a particular example,

    /ORIENTATION, NAME=wire_hori, DEFINITION=COORDINATES
    2.0, 0.0, 0.0, 1.0, 3.0, 0.0, 1.0, 0.0, 0.0
    

    

    Figure 8: Material orientation system defined by points $\vec{a}$, $\vec{b}$ and, if desired, pointn $\vec{c}$, the origin of the local axis system.

  • DEFINITION=NODES
    Similar to DEFINITION=COORDINATES, but the points $\vec{a}$, $\vec{b}$ and ($\vec{c}$) are located at two (three) nodes of the finite element mesh of the analysed RVE. If orientation definition NODES is selected, the user specifies the nodes corresponding to points $\vec{a}$, $\vec{b}$ and ($\vec{c}$) by their global node numbers in the mesh. The input sequence might look like this:

    /ORIENTATION, NAME=wire_vert, DEFINITION=NODES
    12, 34, 54
    

  • DEFINITION=OFFSET TO NODES
    In this special definition of a local axis system, points $\vec{a}$, $\vec{b}$ and ($\vec{c}$) are defined by specifying the local node numbers, relative to an element. By default, the origin corresponds to the first node number of each element belonging to the considered element set. For example, local node number 4 is given as the location of point $\vec{a}$ and local number 5 to define point $\vec{b}$. Then, the local X' axis is defined to be parallel to the (1, 4) side of the hexahedral element and local axis Y' is in the plane formed by the nodes 1, 4 and 5 of this element. In this case, the user introduces simply:

    /ORIENTATION, NAME=fabric, DEFINITION=OFFSET TO NODES
    4, 5
    

    Of course, if the origin of the local axis systems differs from the first element node, the user has to specify also point $\vec{c}$ by the corresponding local origin node. This orientation definition is very flexible and powerful. Indeed, in presence of an orthotropic curved layer, this orientation definition allows the user to specify easily a different local axis system in each element of the considered element set. The curvature of the layer implies such a local orientation definition per element. Note that, as the local axis system varies from element to element, also the Hooke matrix, expressed in the structural axis system, varies in the considered material domain, leading to the generation of implicit body forces, which are taken into account in the homogenization analysis (see Eq. \eqref{eq:solu_2} of Appendix Homogenization of Materials with Periodic Microstructure).

  • ROTX=90 for 2-D MICRESS microstructures
    The additional optional parameter ROTX must be specified on the keyword line in the case of 2-D homogenization of microstructure predicted by MICRESS and no specific orientation file (.TabO file) is transferred to HOMAT. Indeed, in 2-D phase field simulation MICRESS defines the discretization domain in the plane X-Z. Structural axis Y is perpendicular to this plane. Whereas HOMAT specifies always a 2-D plane strain or stress model to homogenize in the classical X-Y plane. Therefore, an additional rotation of $\alpha = 90$° about the common X axis must be achieved to turn the MICRESS model with its grain orientation specified in the $\vec{E}_{X}$, $\vec{E}_{Z}$, $\vec{E}_{Y}$ axis system (see Fig. 9), in the $e_{X}$, $e_{Y}$, $e_{Z}$ axis system of HOMAT. In 2-D, the user must specify this additional rotation via the ROTX parameter:

    /ORIENTATION, NAME=grain_1, DEFINITION=QUATERNION, ROTX=90
    0.880735, 0.287, -0.3321, 0.177857
    

    

    Figure 9: MICRESS 2-D model ($\vec{E}_{X}$, $\vec{E}_{Z}$, $\vec{E}_{Y}$ axis system) rotated around the common axis $\vec{E}_{X} = e_{X}$ in the $e_{X}$-$e_{Y}$ plane of the HOMAT 2-D analysis.

  • OUT_ORIE
    This additional optional parameter permits to define the local axis system in which the calculated effective properties are printed in the result files (PROBLEM_NAME.hres and PROBLEM_NAME.mat) in addition to the structural axis system of the RVE. In fact, it permits to save the name of the actual orientation definition for this additional output.

    /ORIENTATION, NAME=or_grain3, DEFINITION=EULER, OUT_ORIE
    45, 14.0, 60.
    

    The orientation or_grain3, defined by the Euler angles $\phi=45$°, $\theta=14$° and $\psi=60$°, specifies in this example the local axis system for the output of the effective orthotropic or anisotropic material properties.

    Note: The OUT_ORIE parameter can be specified only once for a specific orientation of the whole homogenization model. The locally variable orientation definition OFFSET TO NODES cannot be used to specify the output local axis system.

• /PERIODICITY, NAME = name, TYPE = TRANS | ROTA, NUMBER = number, START = start, TARGET = target

  This important boundary condition must at least be defined in one direction, for example structural direction X, on the microstructure RVE. Classically the RVE is periodic in all three directions. As further outlined in Appendix Homogenization of Materials with Periodic Microstructure, the periodicity B.C. introduces a constraint between displacements belonging to two opposite surfaces separated by a distance called period. Mesh2HOMAT automatically identifies, for each structural direction the appropriate free surfaces. A special nomenclature is introduced: the surface normal to the X, Y or Z directions, whose X, Y or Z coordinates are the smallest, are named X-, Y- or Z-; whereas the surfaces with the largest X, Y or Z coordinates are named the X+, Y+ or Z+ surfaces. Depending on the adopted periodicity assumption, e. g. 3-D, 2-D or 1-D, Mesh2HOMAT prepares the definition of the linked surface pairs. In Ex. 13 of the UD reinforced composite, 2-D periodicity is specified, as the unit cell can be repeated in X and Y direction by a period equal to its extension, but not in the fibre direction, as the fibre is assumed to be endless. In this case, Mesh2HOMAT generates following:

  /PERIODICITY, NAME=X-X+, TYPE=TRANS, NUMBER=1, START=X-, TARGET=X+
  -3.0, 0.0, 0.
  /PERIODICITY, NAME=X+X-, TYPE=TRANS, NUMBER=2, START=X+, TARGET=X-
  3.0, 0.0, 0.
  /PERIODICITY, NAME=Y-Y+, TYPE=TRANS, NUMBER=3, START=Y-, TARGET=Y+
  0.0, -3.0, 0.
  /PERIODICITY, NAME=Y+Y-, TYPE=TRANS, NUMBER=4, START=Y+, TARGET=Y-
  0.0, 3.0, 0.
  

  After the keyword /PERIODICITY;, the parameter NAME specifies the name attributed to the considered link between two periodic surfaces. This name is a composition of the names of the two linked surfaces. The first surface name corresponds here to the surface, which is projected on the other surface, specified by the second name. Then, the parameter TYPE defines how the surfaces are linked. HOMAT allows the introduction of two different transfer modes: translation or rotation.

  • If TYPE=TRANS is specified, a translation vector, defining the period, is introduced on the second line. Its components specify, per direction, the increments need to go from the slave face (1^st surface) to the master surface on the opposite RVE side (2^nd surface). The required parameter START specifies the name of the starting (slave) surface and the parameter TARGET defines the master surface, on which the slave nodes are projected.

  • If TYPE=ROTA is specified, the user provides in the X-Z plane a rotation angle $\theta_{y}$ (in degrees) around the common structural axis Y. This special feature allows to define periodic B.C. on curved cylindrical RVE's, is detailed in reference ¹⁶.

    /PERIODICITY, NAME=X-X+, TYPE=ROTA, NUMBER=1, START=X-, TARGET=X+
    0, -15.7., 0.
    

    Note: at the moment, only a rotation $\theta_{y}$ about the common axis Y is implemented in HOMAT. For each linked pair of periodic surfaces, for example X- and X+, the program Mesh2HOMAT specifies always two /PERIODICITY commands: the 1^st one to define the translation/rotation in one direction and the other one to specify the translation/rotation in the opposite direction. Indeed, it is important to specify each surface one time as slave face and the other time as master face. As 2-D periodicity is defined on the glass-epoxy RVE of Fig. 4, two pairs of linked surfaces, in X and Y directions respectively, are generated with four /PERIODICITY commands in Ex. 13.

• /SOLID SECTION, NUMBER = number, ELSET = elset, NAME = name , MATERIAL = material , [ ORIENTATION = orientation, ] TEMP = NO | HOM | CALC, MECH = NO | HOM | CALC, FLOW = NO | HOM | CALC

  This keyword is used to define physical regions with common properties, defined previously by the /ELSET option and corresponding, for example, to grains of a polycrystal or to regions of the RVE of a heterogeneous material. Some parameters are always required and the ORIENTATION parameter must only be specified if the considered section is transverse isotropic, orthotropic or anisotropic (see for example keyword /ELASTIC). On the command line the user specifies successively the following parameters:

  • ELSET: this parameter is set equal to the name of the element set containing all elements for which the material behaviour is being defined via this option. Note that alternatively, the parameter GROUP can be specified. It has exactly the same meaning.
  • NUMBER: defines the actual position in the list of all specified solid sections. Its count starts from 1 and it is incremented at each new solid section definition.
  • NAME specifies the name of the considered solid section.
  • MATERIAL set this parameter equal to the name of the previously specified material to be applied to this solid section.

  • ORIENTATION for the transverse isotropic, orthotropic and anisotropic materials the user has to specify the name of the previously defined local axis system via keyword /ORIENTATION. Of course, this parameter can be neglected for isotropic materials.

    Remark: In order to reduce drastically the definition effort in the case of MICRESS RVE's with a lot of grains, having each another grain orientation, a direct import of the orientations from the TabO file can be specified via the parameter ORIE of option /MICROSIM (see Sec. Definition of HOMAT- specific Commands). In this special case the parameter ORIENTATION can be omitted in the /SOLID SECTION definition.

  In addition to these parameters common with ABAQUS, the following special HOMAT parameters have to be specified. Their definition is essential as they define if the considered solid section participates in the actual thermal and/or thermoelastic or Stokes flow homogenization run.

  • TEMP: If TEMP=CALC or HOM, then the considered solid section is part of the planed thermal homogenization analysis. If TEMP=NO, this solid section does not participate to the thermal homogenization analysis. This way, the simulation model can be easily reduced to solid sections of interest and omit non interesting regions of the RVE.
  • MECH: this parameter can have also the values: CALC, HOM or NO. In the case CALC or HOM, the considered solid section is part of the planed thermo-elastic homogenization and, of course, in the case NO, it does not participate to the analysis.
  • FLOW: this parameter can have also the values: CALC, HOM or NO. In the case of CALC or HOM, the considered "solid section" is a fluid domain and participates to the Stokes flow homogenization. Please specify for all the solid domains, FLOW=NO, as in these regions no flow happens.

  Note: The new parameter problem of Mesh2HOMAT specifies not only the desired homogenization analysis but also provides, for each homogenization analysis type, a default definition of the parameters TEMP, MECH and FLOW. Of course, the user can modify this definition to suit his needs. In Ex. 13 of the UD reinforced composite a thermo-elastic homogenization analysis is defined. Therefore for both constituents, the fibre and matrix regions, a /SOLID SECTION keyword is specified with TEMP=NO and MECH=CALC or HOM:

  /SOLID SECTION, NUMBER=1, ELSET=FIBER, NAME=fibre, MATERIAL=verre, TEMP=NO, MECH=CALC
  #
  /SOLID SECTION, NUMBER=2, GROUP=RESIN, NAME=matrix, MATERIAL=epoxy,TEMP=NO, MECH=HOM
  

Definition of HOMAT- specific Commands¶

In the following section the HOMAT-specific keywords are listed not alphabetically but in the order they are classically written by Mesh2HOMAT in the HOMAT command file PROBLEM_NAME.hin.

Options introduced before the RVE description¶

• /HEADING

  The command file PROBLEM_NAME.hin always starts with the /HEADING keyword. It has no specific parameters and permits the user to specify the intended homogenization simulation by providing a problem description. When generated by Mesh2HOMAT, the header section of the command file contents 6 lines, where the second line is left blank for the user to add a short problem description. The first line contains the date and time of the Mesh2HOMAT call. The 3^th until 6^th line outlines the total number of elements, nodes, new generated nodes and materials respectively.

  Note: It is important that the total number of nodes and new nodes indicating in the heading corresponds to the number of nodes specified in the PROBLEM_NAME.hgeo file. HOMAT checks this coincidence. Ex. 13 of the thermo-elastic homogenization of a UD reinforced composite starts with following header:

  /HEADING
  Date: 20.11.19 Time: 12:45:24
  Total number of Elements: 296
  Total number of Nodes: 732
  Total number of new Nodes: 76
  Total number of Materials: 2
  

• /PROBLEM, TYPE = type, HOM_MOD = hom_mod, PER_MOD = per_mod

  This option defines the desired homogenization analysis and is the keyword with the largest impact on the predicted effective properties. The nature of the homogenization analysis is specified here via the parameter TYPE. The user has also to choose one of following three homogenization models: 3-D (default), 2-D plane strain or 2-D plane stress model, via the parameter HOM_MOD. The distinction between plane stress and strain is only relevant for thermos-elastic homogenization. If a 2-D thermal or fluid flow homogenization is desired, HOM_MOD = 1 is always set. The user must specify the nature of the imposed periodicity boundary conditions: 1-D, 2-D or 3-D. Since HOMAT 6.0, 3-D periodic BCs are assumed by default. The TYPE specifier can take the following values:

  • TYPE=THERMAL
    A purely thermal homogenization analysis is performed and only the effective thermal conductivities, capacities and densities are evaluated.
  • TYPE=THERMEC
    A full thermo-mechanical homogenization analysis is performed. In this case, a thermal homogenization and a thermo-elastic one are realized successively in the same run. Effective thermal conductivities, capacities and densities are determined as well as the effective Hooke matrix, the effective thermal expansion coefficient and the generalized eigenstrain, if present.
  • TYPE=MECHAN
    A purely mechanical homogenization run is performed. In this case, only the Hooke matrix and its derivate effective properties (Young and shear moduli, Poisson coefficient,...) are evaluated but not the effective thermal expansion coefficient and the volumetric eigenstrain.
  • TYPE=THELAS
    A thermo-elastic homogenization analysis is performed. In this case all effective thermo-elastic properties (Hooke matrix and its derivate, effective thermal expansion coefficient, ..) are determined, but no thermal properties.
  • TYPE=STOKES
    A Stokes flow homogenization analysis is performed. The effective Darcy permeability tensor is evaluated, by applying unitary boost volume forces on the fluid domains.
  • TYPE=LOCALI
    A mechanical localization analysis is performed including a previous mechanical homogenization step in order to calculate the required, microscopic, periodic displacement fields $x^{rs}$ on the RVE but the effective mechanical properties are not evaluated. The micro-stresses and -strains are evaluated on the RVE for a given macro-strain state.
  • TYPE=LOCATH
    A thermo-elastic localization analysis is performed including a necessary thermo-elastic homogenization step. Total micro-strains, including thermal expansion and volumetric eigenstrains, and micro-stresses are evaluated on the RVE for a given macro-strain state.

  The HOM_MOD parameter can take the following values:

  • HOM_MOD = 0
    3-D homogenization model. In this case, for the thermal analysis, three linear problems are solved. They correspond to the microscopic displacement fields $\zeta^{X}$, $\zeta^{Y}$ and $\zeta^{Z}$ respectively. In the case of mechanical homogenization, HOMAT establishes and solves six linear systems, corresponding to the microscopic displacement fields $x^{rs}$ $(x^{XX},\, x^{YY},\, x^{ZZ},\,x^{XY},\, x^{XZ},\, x^{YZ})$ respectively. These fields are induced by the corresponding unitary macro-strain tensor $\vec{E}_{rs}$. In the case of flow homogenization, three special Stokes flow problems are solved providing the velocity fields $\omega^{x}$, $\omega^{y}$ and $\omega^{z}$ and the pressure fields $p^{x}$, $p^{y}$, $p^{z}$ on the RVE respectively. The 3-D homogenization model is specified by default.
  • HOM_MOD = 1
    2-D plane strain homogenization model. In the case of thermal or flow homogenization, two linear problems are solved in the X-Y plane. They provide either the planar microscopic displacement fields $\zeta^{X}$, $\zeta^{Y}$ in the thermal analysis or the planar velocity and pressure fields ($\omega^{x}$, $\omega^{y}$ and $\omega^{z}$ and $p^{x}$, $p^{y}$). 2-D plane strain mechanical behaviour ($\epsilon_{ZZ}=0$ but $\sigma_{ZZ} \neq 0$, corresponding to a section through an infinite thick component) is defined in the mechanical homogenization leading to solve four linear problems, corresponding to the displacement fields $x^{XX}$, $x^{YY}$, $x^{ZZ}$,$x^{XY}$.
  • HOM_MOD = 2
    2-D plane stress homogenization model. In the case of thermal and flow homogenization, this model equals the previous one. No difference between plane strain and plane stress assumptions is made there. In the case of mechanical homogenization, the plane stress assumption ($\sigma_{ZZ}=0$ and $\epsilon_{ZZ} =0$, thin membrane behaviour) leads to solve only three linear systems, providing the displacement fields $x^{XX}$, $x^{YY}$, $x^{XY}$.

  The PER_MOD parameter determines the nature of the applied periodic boundary conditions:

  • PER_MOD = 1
    1-D periodic boundary condition is specified. In this case, two opposite free surfaces are linked. Classically, they are located normal to the X direction, but they can also be defined in other direction. The keyword /PERIODICITY defines the pair of linked surfaces.
  • PER_MOD = 2
    2-D periodic boundary conditions are defined. In this case two pairs of linked surfaces must be specified for an RVE of cubic or rhombic shape via the /PERIODICITY option. Classically, the free surfaces in X direction are linked together as well as both free surfaces in Y direction. But, it is also possible to specify other linked pair of opposite surfaces, e.g. in X and Z directions or in Y and Z directions.

  • PER_MOD = 3
    3-D periodic boundary conditions are applied on the RVE. In this case, three pairs of linked surfaces must be specified via the /PERIODICITY option.

    Note: At the moment, Mesh2HOMAT automatically generates the periodic boundary conditions only for RVEs of cubic or rhombic shape but not for unit cells of hexagonal or octahedral shape (see Fig. 10 and 11).

    

    Figure 10: Unidirectional reinforced composite with hexagonal unit cell.

    

    Figure 11: 3-D finite element mesh of hexagonal RVE with three pairs of periodic faces and 2-D periodicity.

  The SOLVER option specifies the adopted solver for the discretized linear global system.

  • SOLVER=ITER
    An improved iterative BiCGstab solver is used, by default, to solve the linear systems of the mechanical, thermal or flow problems. This solver is coupled with a symmetric pre-conditioner and uses, since version 6.0, a sparse matrix storage format for the system matrix.

    Note: At the moment, the iterative BiCGstab solver is the only available solver for the thermal, mechanical and fluid flow homogenization problems.

    In Ex. 16, the 2^nd part of the command file for the thermo-mechanical homogenization of a UD reinforced NiAl composite with Al₂O₃ (sapphire) fibres is outlined. It is the part of the command file with all specific HOMAT keywords, which are introduced after the keywords describing the RVE model. The /PROBLEM keyword indicates there that successive thermal and thermo-elastic homogenizations will be performed on a RVE with square shape on which 2-D periodic boundary conditions are applied.

  #
  /PROBLEM, TYPE=THERMEC, HOM_MOD=0, PER_MOD=2
  #
  /DEF LAYER, TYPE=MULTI, NUMBER=1
  LAYER=1, NAME=RVE, NB_DOM=ALL
  #
  /TEMP DOMAIN, NUMBER=2, NTEMP=2
  DOMAIN=1, 800.0, 1100.0
  DOMAIN=2, 800.0, 1100.0
  #
  # Solver settings
  #
  /NUMERICS, ITMAX=2000, PRECISION=5.0D-05, PENAL=1.0D+05,
  ITCRIT=1, OUTFREQ=20, EPSCOORD=1.0D-02, TINY=1.0D-16
  #
  # Parameter definition for the homogenization module
  #
  /HOM PARA, MINTEMP=0.0, MAXTEMP=1400.0, CUBVAL=0.1, TREF=20.0
  #
  /OUTPUT, FORMAT=VTK, TYPE=ASCII, IOP=2, DBG=YES, VOLAVRG=YES, VISLAV=YES
  T, U
  #
  /GAUSS INTEGRATION,* TYPE=HEX, N_X=2, N_Y=2, N_Z=2
  #
  /GAUSS INTEGRATION*, TYPE=PENT, N_T=3, N_Z=2
  /GAUSS INTEGRATION*, TYPE=TET, NUMBER=4
  #
  /GAUSS INTEGRATION,* TYPE=QUAD, NUMBER=4
  /GAUSS INTEGRATION,* TYPE=TRIA, N_T=3
  #
  

  Example 16: HOMAT specific keywords defined for the thermal and thermo-elastic homogenization of a UD reinforced Al₂O₃-NiAl composite.

• /OUTPUT

  This keyword permits the user to define the output of the considered homogenization or localization analysis. Sec. HOMAT Output details the different output possibilities for homogenization and localization analyses. The parameters of the keyword line are the following:

  • FORMAT: this parameter is set by default to VTK. In fact, the VTK format is used as standard exchange protocol between HOMAT and other programs, working on the RVE level, such as MICRESS or SphaeroSim. Moreover, HOMAT relies on the open-source, multi-platform visualization program Paraview to visualize the implicit fluxes/forces and periodic displacement, micro-strain, micro-stress, velocity and pressure fields.

  • TYPE: the VTK output files can be either ASCII or BINARY. The BINARY option allows to reduce the size of the result files and is the default setting if the HOMAT analysis is performed after a microstructure simulation with MICRESS or SphaeroSim. For other homogenization and localization runs, the ASCII option is adopted as the default TYPE (see Ex. 16).

  • DBG: this parameter triggers additional control output either on the screen or in a log file, named PROBLEM_NAME.hout. By default, this parameter is set to YES. The activation of this output option is recommended as it does not produce much overhead, but allows to follow the individual program steps in more detail.

  • IOP: this option controls the level of printed intermediate results, like the Hooke matrix per grain/material, the evaluated implicit forces, the element stiffness matrix and so on either on the screen or in the PROBLEM_NAME.hout file. This output control parameter varies between 0 and 5 and a higher value produces more output. If IOP=0, the screen result output is strictly limited. This value corresponds to the default option for large scale applications. Note that IOP=2 is a very frequently used choice and recommended if the user want to control the main intermediate output results. Values above IOP=2 are only needed to debug the program.

  • ENG: this parameter allows, if ENG=YES, the direct output in the .mat file of the effective orthotropic Young and shear moduli and of the Poisson coefficients. This output is not classical ABAQUS one, but can be used by other macro-scale simulation tools. By default, ENG=NO.

  • VOLAVRG: if VOLAVRG=YES, then in the effective properties result file, PROBLEM_NAME.hres, not only the effective thermal and/or thermo-elastic properties are displayed; but also the corresponding mean volume averaged properties of the layers (monolayer homogenization) or of the whole RVE (multilayer homogenization). This way, the user can better assess the accuracy increase, induced by the correction term of the first order homogenization [cf. Eq. \eqref{eq:solu_3}]. Mesh2HOMAT writes in the template of the command file PROBLEM_NAME.hin the default value VOLAVRG=YES.

  • POLE: this parameter allows to draw the pole figure in the plane X-Y of the lamella orientations perpendicular to the radial growth direction in each spherulite of semi-crystalline thermoplastic microstructure. The pole figure are realized by stereographic projection of the orientations on the plane X-Y (see Fig. 12).

  • VISLAV this parameter permits the accurate visualization of forces and fluxes within Paraview. It allows to regenerate the original values of the forces and fluxes in the master nodes, which in the solved reduced system content not only its value but also the values of corresponding linked slave nodes. By default, VISLAV=YES.

  • LOC_AXIS: if LOC_AXIS=YES, then the effective material properties are not only written in the structural RVE axis system but also in the specified local axis system in both result files: PROBLEM_NAME.hres and PROBLEM_NAME.mat. Of course, the local axis system has to be specified before via the /ORIENTATION or the /MICROSIM keyword. By default, LOC_AXIS=NO.

  • SHEAR_TENS this parameter governs the output of the shear-tension coefficients, which can be derived from a fully anisotropic effective Hooke matrix. If SHEAR_TENS=YES, these shear-tension coefficients are calculated and printed in the PROBLEM_NAME.hres file. In order to reduce the output of fully anisotropic Hooke matrices in this file, SHEAR_TENS=NO by default.

  In the data line of keyword /OUTPUT the user has to specify the nature of results fields on the RVE, which will be saved by HOMAT in the VTK format. A detailed description of the results files and their nomenclature is given in Sec. HOMAT Output. The following options are available:

  • ALL: This option is self-explanatory. All the evaluated microscopic fields are saved for possible visualization. By specifying this option, the user does not need to select the microscopic displacement and flux fields individually, results of the thermal homogenization via parameter T or/and the displacement, micro-strain, pressure and implicit force fields, results of the thermo-elastic homogenization, via parameters U, P and E. Both are saved in the case of a combined thermal and a thermo-elastic homogenization. For a localization run, the microscopic strain and stress fields are both saved in separate files. For a Stokes flow homogenization, the boost forces, the velocity and pressure fields are saved in separated files.

  • T: In the case of thermal homogenization, the microscopic displacement and implicit flux fields are saved in separated files for visualization purpose.

  • U: In the case of thermo-elastic or mechanical homogenization, the microscopic displacement fields and the implicit forces are saved in separated files for possible visualization.
  • E: In localization, mechanical as well as thermo-elastic homogenization analyses, the variation of the microscopic engineering strain tensor over the considered layer or over the whole RVE are saved for visualization and/or post-treatment. In the case of localization, the equivalent strain is also saved for possible visualization.

  • S: The components of the microscopic stress tensor and the equivalent stress is saved per layer or per RVE for visualization purpose in a localization run.

  • P: In presence of mixed finite elements (thermo-elastic homogenization involving quasi-incompressible materials or Stokes flow homogenization), the parameter P indicates that the local pressure fields in these elements are saved for visualization purpose.

  • V: The components of the velocity fields, for Stokes flow problems, are saved per layer or RVE for visualization.

  If the parameter POLE=YES is set, then the user has to specify via the additional parameter NB_SPHERU on the third line (second data line) the number of spherulites for which the pole figured will be build. On the following line, the user specifies the list of spherulites for which a specific file containing the stereographic projection points of the lamella orientations, perpendicular to the radial growth direction, on the X-Y plane is built. The identification number of a spherulite is the NUMBER value specified in the /SOLID SECTION definition. In the spherulite example, mscal_PP, we have introduced following output definition:

  /OUTPUT, FORMAT=VTK, TYPE=BINARY, DBG=YES, IOP=2, VOLAVRG=YES, ENG=NO, POLE=YES
  U
  NB_SPHERU=3
  80, 103, 104
  

  In Figure 13, the pole figure for the spherulite 103, shown in Fig. 12, is presented.

  

  Figure 12: Voxel mesh of spherulite 103.

  

  Figure 13: The pole figure of lamella orientations in the X-Y plane for spherulite 103.

• /UNITS, GEOMETRY = sfac, LENGTH UNIT = ME | CM | MM | MU | NM, TEMP OFFSET = offset, ANGLE = GRAD | RAD

  This keyword permits the user to specify the desired unit system for the homogenization run. Following optional parameters can be specified on the keyword line:

  • GEOMETRY
    This parameter allows the introduction of a scaling factor, which will be applied to all node coordinates. By default, no scaling is defined: GEOMETRY=1.
  • LENGTH UNIT
    This parameter, linked to the GEOMETRY parameter, can take either the value ME for meter, CM for centimeter, MM for millimeter, MU for micrometer and NM for nanometer. It allows HOMAT to print, for example, the solid section volumes using correct units. MM is the default value of this parameter.
  • TEMP OFFSET
    This parameter allows to change easily the temperature unit from grad Celsius to Kelvin. By default, the temperature is noted in °C. The offset is then: 273.15K.
  • ANGLE
    This parameter allows to specify the unit for angles: either as degree (ANGLE=GRAD) or as radians (ANGLE=RAD). By default, an angle is defined in degree.

  In Ex. 13, following unit definitions are introduced:

  /UNITS, GEOMETRY=1, TEMP OFFSET=273.15, ANGLE=GRAD
  

  The geometry of the RVE is not scaled. Therefore, it was not necessary to specify the length unit as the default one is used. The angles are specified here in degree. If the user wants to use in the simulation microns as length unit and the original mesh is in millimeters, the following /UNITS parameters are specified:

  /UNITS, GEOMETRY=1.0E+3, LENGTH UNIT=MU, TEMP OFFSET=273.15, ANGLE=GRAD
  

• /INCLUDE, NAME = name

  The /INCLUDE keyword can be used at any time in the command file in order to include an external file that contains a part of the HOMAT command file, like the PROBLEM_NAME.hgeo file, which describes the geometrical model of the heterogeneous material, whose effective properties will be evaluated. Note that, in order to reduce the size of the geometry file, Mesh2HOMAT writes, by default, a gzipped geometrical model. When a reference to an external file is encountered, HOMAT will unzip the file and process the data within the specified external file. At the end of this external file, HOMAT will automatically continue to interpret the calling file, which is generally the command file PROBLEM_NAME.hin. The maximum number of nested /INCLUDE operations is 4. This keyword has only one mandatory parameter NAME. This parameter specifies the name of the external file. In Ex. 13, the RVE geometry and the thermo-physical material properties of the fibre and matrix are specified externally and included via the option /INCLUDE, as follows:

  /INCLUDE, NAME=ehz50_3D.hgeo.gz
  #
  /MATERIAL, NAME=VERRE, TYPE=ANISO
  /INCLUDE, NAME=verre_ani.hmat
  #
  /MATERIAL, NAME=EPOXY, TYPE=ANISO
  /INCLUDE, NAME=epoxy_ani.hmat
  

• /MICROSIM, NAME = name, TIME = time, [ ORIE = YES | NO ], [ OUT_ORIE = grain ]

  This keyword tells HOMAT that the microstructure to homogenize comes from a previous simulation with MICRESS. This option allows an easy and direct definition of the grain orientations of a polycrystal. In fact, the MICRESS problem name is known by HOMAT via the mandatory parameter NAME. The user has to provide in the working directory the ASCII file MICRESS_NAME.TabO containing the mean orientation per grain. As this text file contents orientation information for the whole microstructure evolution, it is essential to specify to HOMAT the MICRESS simulation time for which effective properties will be evaluated. This information is given via the mandatory parameter TIME. Setting parameter ORIE=YES, indicates to HOMAT that, by default, the grains/regions of the microstructure receive their orientation directly via the TabO file. Thus, the orientation definition via the keyword /ORIENTATION and the parameter ORIENTATION of each solid section can be omitted. However, if this keyword is specified for a specific solid section, than its orientation definition overrides the default setting via the TabO file.

  Remark: In the case of localization analysis with different structural axis system at micro- and macro-scale, the /MICROSIM keyword must be specified before the keyword /ORIENTATION in the .hin command file, if the user wants to specify the grain orientations easily via the TabO file.

  In the presence of a 2-D MICRESS microstructure, the parameter ROTX=90 must be specified. It defines the rotation angle between the MICRESS axis system $\vec{E}_{X}$, $\vec{E}_{Z}$, $\vec{E}_{Y}$ to the HOMAT structural axis system $e_{X}$, $e_{Y}$, $e_{Z}$ (see Fig. 8). In the 2-D case, HOMAT rotates, first, each vector or tensor quantity from the local crystal axis system to the MICRESS axis one and then from this system to its own reference system. This parameter has no relevance in 3-D.

  The new optional parameter OUT_ORIE permits to specify the name of the grain orientation which specifies the local axis system for the output of the effective material properties in this local axis system. Only one grain can be selected for this additional output and a definition via the /ORIENTATION keyword is obsolete.

  For the austenite to ferrite phase transformation of FeCMnSi steel, discretized in 3-D and detailed further in the example manual, the user introduces the following parameters:

  /MICROSIM, NAME=dim1.5, TIME=42.0, ORIE=YES, OUT_ORIE=grain2
  

  The MICRESS name corresponds to the driving file name. Here, the microstructure after 42sec, containing 49.46% ferrite grains, is homogenized and all grain orientations are directly provided via the TabO file. Moreover, the effective thermo-elastic properties are additionally written in the local axis system of grain2.

Options introduced after the RVE model description¶

• /DEF LAYER, TYPE = MONO | MULTI, NUMBER = number, LAYER = , NAME = name, NB_DOM = dom

  In presence of a RVE build of different material layers, like the cooled multilayer build of a substrate of superalloy CMSX-4, a MCrAlY bondcoat and a yttrium stabilized thermal barrier coating (see Fig. 14), the keyword /DEF LAYER allows to select either multilayer or monolayer homogenization via the keyword parameter TYPE. Multilayer homogenization predicts the effective properties of the whole RVE; whereas in the case of monolayer homogenization, the user can specify individually the layers, which are homogenized in an analysis. Of course, in presence of only one material layer, both homogenization types are identical. Note that the multi- or monolayer homogenization option can be specified in combination with any kind of homogenization analysis.

  • TYPE=MULTI: multilayer homogenization
    By default, TYPE is set to MULTI. In this case, the other keyword parameters are easily pre-defined by Mesh2HOMAT, like shown as following:

    /DEF LAYER, TYPE=MULTI, NUMBER=1, LAYER=1, NAME=RVE, NB_DOM=ALL
    

    Only one layer is defined (NUMBER=1 or 0), the number of the considered LAYER is always 1 and the layer NAME is set to RVE and the parameter NB_DOM, specifying the number of considered domains or grains is set to ALL.

  • TYPE=MONO: monolayer homogenization
    The number of layers considered by the homogenization run is introduced via parameter NUMBER. If you desire the effective monolayer properties of each layer, the user introduces the list of the concerned material domains, as shown below for the cooled multilayer plate of Fig. 14:

    /DEF LAYER, TYPE=MONO, NUMBER=3
    LAYER=1, NAME=TBC_HO, NB_DOM=6
    1, 2, 3, 4, 5, 16
    LAYER=2, NAME=BOND_HO, NB_DOM=6
    6, 7, 8, 9, 10, 17
    LAYER=3, NAME=SUBSTRATE, NB_DOM=6
    11, 12, 13, 14, 15, 18
    

    The number of materials domains per layer is specified via the parameter NB_DOM. Each layer of the cooled plate has 6 domains, the solid and 5 cooling holes. On the next line, the user has to specify, via their identification number, the list of the domains which belong to the layer. The identification number of a material domain is the NUMBER value specified in the /SOLID SECTION definition.

  Remark: At the moment, HOMAT v6.0 allows to perform monolayer homogenization on layers of a multilayer RVE with 2D or 1D periodic boundary conditions. 3D periodic BCs on individual layers are not allowed.

  

  Figure 14: RVE of a cooled multilayer plate composed of a thermal barrier coating, a bondcoat and a superalloy substrate layer. The cooling channels are inclined by a 30° angle.

• /TEMP DOMAIN, NUMBER = number, [ NTEMP = ntemp ]

  This keyword allows the definition of the domain temperature: either this temperature can vary inside each domain or the temperature is constant per domain. In the 1^st case a temperature field is provided, coming from a previous analysis via an input file (see option /INPUT FILES). In the 2^nd case, the user has to specify, for each domain, the temperatures at which homogenization will be performed in one run. This keyword requires the definition of the parameter NUMBER, which specifies the nature of the temperature field.

  1. In the first case, the user specifies only following keyword line:

     /TEMP DOMAIN, NUMBER=0
     

     Parameter NUMBER equals zero is directly interpreted by HOMAT that all domains have no constant temperature and the required local temperature field is imported via an input file. Only one homogenization simulation is performed with this set temperature profile on the RVE.

  2. In the second case, the user has, at first, to define via parameter NUMBER the number of material domains or grains affected by the present temperature definition. Per domain, a special entry line has to be introduced, as shown below. Further, the user must introduce the number of different temperatures, specified for each concerned material domain or grain, via the keyword parameter NTEMP. If NTEMP differs from one, HOMAT will perform automatically NTEMP homogenizations in one simulation run. This option is very interesting in applications where the temperature dependence of the effective properties is investigated. This feature allows to perform NTEMP homogenizations in a single run. At the moment, HOMAT limits NTEMP to 9 simultaneous homogenizations. An example is given below:

     /TEMP DOMAIN, NUMBER=2, NTEMP=3
     DOMAIN=1, 27., 100., 200.
     DOMAIN=4, 27., 120., 220.
     

     The keyword line states that two domains with three different temperature sets are specified. Per domain, the three set temperatures are introduced on each data line: The domains are identified by the parameter NUMBER of the /SOLID SECTION keyword. Of course, the temperatures can be different, per domain/grain. HOMAT will perform a first homogenization at 27°C for both domains, then a 2^nd homogenization with domain 1 at 100°C and domain 4 at 120°C. Finally, it performs a third homogenization with domain 1 at 200°C and domain 4 at 220°C. This simple example shows the great flexibility of defining uniform temperatures per domain for different homogenization analyses. In order to condense the required input, mainly in presence of a polycrystal with several grains set to constant temperatures, the parameter NUMBER can be set to ALL. In this case, the definition for each domain separately can be replaced by one common set of temperatures:

     /TEMP DOMAIN, NUMBER=ALL, NTEMP=2
     27., 100.
     

  In this example, two homogenization runs are specified with all domains set to 27.0°C and 100.0°C respectively. For a MICRESS microstructure, Mesh2HOMAT sets by default, NUMBER=ALL.

• /INPUT FILES, INTEMP = intemp | MACRO = macro , [ ORIENTATION = orientation ]

  This keyword permits to specify the name of different input files coming from linked, previous simulations, like a multi-phase field analysis of metallic alloys with MICRESS. At the moment, only two parameters can be specified: INTEMP and MACRO. They are briefly described below:

  /INPUT FILES, INTEMP=ReX_orien.vtk
  

  The user introduces with this command line the file name of a temperature distribution over the RVE as initial condition for the homogenization run. This set temperature field may vary locally inside each domain/grain. The temperatures are there nodal values, specified via the POINT data format (see VTK documentation).

  /INPUT FILES, MACRO=/your/path/to/U+O_forming.txt
  

  The parameter MACRO defines the path to the ASCII file containing the macro-strain results of the linked thermo-mechanical macro-simulation. A python script, named aba2hom, extracts these macro strains from the ABAQUS results database (OBD file ⁶) in critical areas of the analyzed structure. HOMAT uses these strain states to perform localization runs. The different contents of this file are outlined in Sec. Transfer of macro-scale results from ABAQUS to HOMAT, dedicated specially to this topic.

  Note: The user can introduce either the local file name with extension in the working directory or the global name including the total path name, as outlined above. If the structural axis system in HOMAT differs from the used structural axis in the macro-scale analysis, the user must specify the optional parameter ORIENTATION on the keyword line:

  /INPUT FILES, MACRO=Critical_orth_lin.txt, ORIENTATION=rot_macro
  

  In this case, the orientation definition rot_macro allows HOMAT to rotate the macro-strain results of ABAQUS in the structural axis system of the microstructure in order to perform the localization analysis.

• /NUMERICS

  This keyword specifies numerical parameters for the iterative BiCGStab solver ¹⁷. Mesh2HOMAT has endowed these parameters with a default value, which the user may adapt to his specific problem via the present command. Following parameters can be specified via this keyword:

  • ITMAX
    The number of maximum allowed iterations for the iterative BiCGStab solver to converge to the unknown periodic displacement fields $\zeta^{r}$ or $x^{rs}$ or velocity fields $\omega^{r}$. If the number of maximum iterations is reached without fulfilling the convergence criterion, HOMAT takes the displacement field corresponding to the smallest residuum as the result. By increasing the maximum iteration number, the solver accuracy is increased but, in the same time, the CPU time increases. By default, ITMAX is set to 2000.

  • PRECISION
    This parameter specifies the value, which the convergence criterion of the iterative BiCGStab must undershoot in order to converge. By default, in conjunction with the relative convergence criterion (ITCRIT=1) this convergence parameter is set to PRECISION=5.10^{-5}.

  • ITCRIT
    This option specifies the adopted convergence criterion of the BiCGstab solver. By default, the relative convergence criterion ITCRIT=1 is specified

    $$\begin{equation} \label{eq:relconv} \mathbf{\text{crit}}\left( \mathbf{x}^{i} \right) = \frac{\left| \mathbf{A}\mathbf{x}^{i} - \mathbf{b} \right|^{2}}{\left| \mathbf{A}\mathbf{x}^{0} - \mathbf{b} \right|^{2}} < \mathtt{PRECISION}\,, \end{equation}$$

    where $\vec{A}$ is the symmetric system matrix and $\vec{b}$ the right hand side vector of the linear system $\vec{A}\cdot\vec{x} = \vec{b}$ and $x^{0}$ and $x^{i}$ are the initial guess of the displacement field and its expression at iteration i respectively. If ITCRIT=2, following absolute convergence criterion is adopted

    \begin{equation} \label{eq:absconv} \mathbf{\text{crit}}\left( \mathbf{x}^{i} \right) = \ \left| \mathbf{A}\mathbf{x}^{i} - \mathbf{b} \right|^{2} < \mathtt{PRECISION} \,. \end{equation} - OUTFREQ
    This option defines the frequency of printing intermediate solver results. The default value OUTFREQ=20 means that after 20 iterations, HOMAT prints the value of the actual convergence criterion, the norm of the actual residuum: $\vec{r}^{i} = \vec{A}\cdot\vec{x}^{i} - \vec{b}$, the mean displacement or velocity and the best iteration. - EPSCOORD
    This numerical parameter is used in the implementation of periodic BCs to determine if a node has a direct correspondence on the master surface. This numerical parameter specifies a relative precision, which is multiplied by the mean element size of the slave face in order to define the local absolute precision. HOMAT ensures that the absolute precision is always smaller than the smallest element size on the considered free surface. By default, EPSCOORD=1.D-02.

  • IPRE
    This parameter specifies the adopted pre-conditioning scheme. IPRE=1 means that a diagonal matrix $\vec{D}$ built of the square roots of the diagonal terms of system matrix $\vec{A}$ is applied in symmetric way to the system matrix $\vec{A}$:

    $$\begin{equation} \label{eq:pre1} \vec{D}\cdot\vec{A}\cdot\vec{D}\cdot\vec{x} = \vec{D}\cdot\vec{b}\cdot\vec{D}\,,\text{ with } D_{ii} = \frac{1}{\sqrt{\left|A_{ii}\right|}} \end{equation}$$

    In the case IPRE=2 the precondition matrix $\vec{D}$ is built from the absolute value of each diagonal term of matrix $\vec{A}$ and only a left-hand side pre-conditioning is performed:

    $$\begin{equation} \label{eq:pre2} \vec{D}\cdot\vec{A}\cdot\vec{x} = \vec{D}\cdot\vec{b}\,,\text{ with } D_{ii} = \frac{1}{\left|A_{ii}\right|} \end{equation}$$

    Numerical investigation has shown that the new symmetric pre-conditioning outperforms the left-side one, used in HOMAT v5.x versions. Therefore, IPRE=1 is the default option.

  • TINY
    This parameter defines a small number, which is used throughout HOMAT as numerical equal zero. By default, TINY=1.0E-16.

  • ADMIS
    This parameter permits to specify the max. admissible increase of the mean solution (displacement or velocity) between two successive iterations: If the actual mean solution exceeds its admissible increase, a warning message is written in the log file PROBLEM_NAME.hout. Mainly in the first iterations such a warning appears. If such message appears more than five times during the solution of a specific linear system, a detailed analyzed of the problem and its definition is recommended. By default, ADMIS=200.

• /HOM PARA

  This keyword allows the specification of different global parameters, which are used throughout the homogenization analyses.

  • MINTEMP defines the lowest temperature (in °C) used in the definition of all temperature dependent thermo-physical properties. By default, MINTEMP=0.

  • MAXTEMP defines the maximum temperature (in °C) used in the definition of all temperature dependent thermo-physical properties. By default, MAXTEMP=1400.

  • CUBVAL defines the precision value used in the criterion for detecting cubic symmetry in the RVE coordinate system. If

    $$\begin{equation*} \frac{\left| H_{\text{ii}} - H_{\text{mean}} \right|}{H_{\text{mean}}} < \mathtt{CUBVAL}\,, \end{equation*}$$

    the equivalent material is considered cubic²⁷. Here $H_{\text{mean}} = \frac{H_{\text{ii}}}{3}$, $i = 1,3$ is the mean value of the diagonal Hooke matrix terms $H_{ii}$. In this case, the effective Young's modulus E, the elastic constants $C_{11}$, $C_{12}$ and $C_{44}$ and the anisotropic Zener coefficient $A^{Z}$ [see Eq. \eqref{eq:eff_zener}] are provided in addition to its effective orthotropic Young's and shear moduli $E_{i}$, $G_{ij}$ and Poisson coefficients $\nu_{ij}$ in the PROBLEM_NAME.hres output file. If

    $$\begin{equation*} \frac{\left| H_{ii} - H_{\text{mean}} \right|}{H_{\text{mean}}} > \mathtt{CUBVAL}\,, \end{equation*}$$

    the equivalent material is orthotropic. Only its effective Young moduli, $E_{i}$, per structural direction i, the shear moduli $G_{ij}$ per plane i-j and the Poisson coefficients $\nu_{ij}$ are given in the .hres result file. By default, CUBVAL is set to 0.01.

  • HTR_IDEAL specifies the value of an ideal heat transfer between two materials. By default, this value is set to 0.01 W/(mm²K). This value is used for all interfaces with no specified ideal heat transfer coefficient.

  • TREF
    This parameter specifies the reference temperature in °C for the homogenization and localization runs. This reference temperature plays an important role in the definition of eigenstrains via molar volume of the different phases/materials, as it defines the initial state without any eigenstrain in the heterogeneous material. In the thermo-elastic localization analysis this temperature specifies also the reference strainless state. By default, TREF=20.

  • FOAMBOX
    If FOAMBOX=YES, then the user has the possibility to define a box inside the RVE where the velocities inside the porous medium is used to express the Darcy permeability. Velocity values outside this box are ignored in the permeability evaluation. This way, the impact of the periodic boundary condition on a unit cell of a material with random microstructure and of delicate pore cutting, like for open-cell metallic foams, is reduced in the permeability prediction ¹¹. By default, this parameter is set to NO. If FOAMBOX=YES, the user has to specify on the next line the extension of the foam box in the adopted unit system:

    /HOM PARA; MINTEMP=0, MAXTEMP=1400., FOAMBOX=YES
    0.1, 0.2, 0.2
    

  • RANDOM
    This parameter specifies the seed used by the random number generator. Random numbers are defined, for example, by the spherulite model to generate, for each spherulite, the initial orientation of the bilamella of crystalline and amorphous phases in the plane perpendicular to the radial growing direction. If no random seed is specified by the user, HOMAT sets it to 256.

• /GAUSS INTEGRATION

  The finite element mesh of the RVE is composed of hexahedral, pentahedral and/or tetrahedral elements. The variation of the field (temperature, ..) inside an element is approximated either via a linear or quadratic shape function. The evaluation of quantities over the element volume, like the element stiffness matrix, is performed numerically via the Gauss-Legendre method, which requires, per element type, the definition of integration points, named Gauss points. HOMAT offers the possibility to choose between a full or a reduced integration schemes ¹⁰. For each element type, the default definition of the Gauss integration scheme set by Mesh2HOMAT is as follows:

  • linear and quadratic hexahedron: $2 \times 2 \times 2$ Gauss points, which means that 2 Gauss points are defined per structural direction, as shown in Fig. 15(a).
  • linear and quadratic pentahedron: $3 \times 2$, means 3 Gauss points in planes parallel to both triangular faces and 2 Gauss points in the element thickness, defining the Z position of these two planes [see Fig. 15(b)].
  • linear and quadratic tetrahedron: 4 Gauss points, located as shown in Fig. 15©.
  • linear and quadratic quadrangular element: $2 \times 2$ Gauss points, located as shown in Fig. 15(d).
  • linear and quadratic triangular element: 3 Gauss points, located as shown in figure Fig. 15(e).

  As reduced integration can lead to the development of hourglass modes (known as zero energy modes) ¹⁰ and no hourglass control is implemented in HOMAT, reduced integration, like using only one Gauss point for a linear hexahedral element, is extremely delicate to adopt. Therefore, the modification of the number of default integration points is reserved to advanced users, having profound knowledge about the finite element method. Per element type, the user has following possibilities for integration:

  • linear hexahedron element: $1 \times 1 \times 1$ or $2 \times 2 \times 2$ Gauss points;
  • quadratic hexahedron element: $2 \times 2 \times 2$ or $3 \times 3 \times 3$ Gauss points;
  • linear pentahedron element: N_T = 1, 3 or 4 and N_Z = 1 or 2, where N_T are the number of GP defined per triangular face and N_Z the number in the thickness direction;
  • quadratic pentahedron element: N_T = 3, 4, 6 or 7 and N_Z = 2 or 3 Gauss points;
  • linear tetrahedral element: 1 or 4 Gauss points;
  • quadratic tetrahedral element: 4 or 5 Gauss points;
  • linear quadrangular element: $1 \times 1$ or $2 \times 2$ Gauss points
  • quadratic quadrangular element: $2 \times 2$ or $3 \times 3$ Gauss points
  • linear triangular element: 1 or 3 Gauss points
  • quadratic triangular element: 3 or 4 Gauss points

  Note: Independent of the element type, the new stabilized mixed elements (equal order, LiHe, mini or bubble) have the same Gauss point definition as the corresponding linear element version.

  

  Figure 15: Per element type, location of the Gauss points: a) linear hexahedral element b) linear wedge element c) linear tetrahedral element d) linear quadrangular element and e) linear triangular element.

  To specify the integration scheme in HOMAT, the user must define on the keyword line the concerned element type: TYPE=HEX, TYPE=PENT or TYPE=TET for a hexahedral, pentahedral or tetrahedral volume element, respectively, or TYPE=QUAD or TRIA for a quadrangular or triangular surface element. Then, he has the choice to specify the number of integration points directly via the parameter NUMBER or by specifying the number of integration points per direction. For example, he can specify: N_X=2, N_Y=2 and N_Z=2. Note that the intrinsic directions of each element type are defined by the element node numbering (see Figs. 5 and 6). Examples 17 and 18 illustrate different equivalent Gauss point definitions for linear elements:

  #
  /GAUSS INTEGRATION, TYPE=HEX, N_X=2, N_Y=2, N_Z=2
  #
  /GAUSS INTEGRATION, TYPE=PENT, N_T=3, N_Z=2
  /GAUSS INTEGRATION, TYPE=TET, NUMBER=4
  #
  /GAUSS INTEGRATION, TYPE=QUAD, NUMBER=4
  /GAUSS INTEGRATION, TYPE=TRIA, N_T=3
  #
  

  Example 17: Gauss point definition for linear volume and surface elements.

  #
  /GAUSS INTEGRATION, TYPE=HEX, NUMBER=8
  #
  /GAUSS INTEGRATION, TYPE=PENT, NUMBER=6
  /GAUSS INTEGRATION, TYPE=TET, N_T=4
  #
  /GAUSS INTEGRATION, TYPE=QUAD, N_X=2, N_Y=2
  /GAUSS INTEGRATION, TYPE=TRIA, NUMBER=3
  #
  

  Example 18: Alternative, equivalent Gauss point definition for linear and surface elements.

• /EIGENSTRAIN, NUMBER = number

  This keyword must be specified in thermo-elastic homogenization, if the user wants to take isotropic volumetric eigenstrains into account. Eigenstrains can be introduced in HOMAT in two different ways:

  • In presence of a temperature dependent molar volume of each involved material or phase, the user has only to specify on the keyword line via the parameter MATRIX, the material or phase name of the matrix. All other materials or phases are then considered as inclusions, producing a volumetric eigenstrain. The reference temperature, TREF, specified by the /HOM PARA keyword, defines the temperature $T_{\text{ref}}$ at which it is assumed that no eigenstrains exist in the heterogeneous material. Thus, it defines the starting temperature for eigenstrain development in all inclusions via expression:

    $$\begin{equation} \label{eq:eigenstrain} \epsilon_{\text{eigen}} = \frac{\sqrt[3]{M^{\text{vol}}_{\text{incl}}(T)} - \sqrt[3]{M^{\text{vol}}_{\text{matrix}}(T)}}{\sqrt[3]{M^{\text{vol}}_{\text{matrix}}(T)}} + \frac{\sqrt[3]{M^{\text{vol}}_{\text{matrix}}(T)} - \sqrt[3]{M^{\text{vol}}_{\text{matrix}}(T_{\text{ref}})}}{\sqrt[3]{M^{\text{vol}}_{\text{matrix}}(T_{\text{ref}})}} \end{equation}$$

    This easy definition is recommended to be used in combination with a MICRESS simulation, for example, by solid state transformation of a steel microstructure, as the user has only to specify the matrix material. - constant eigenstrain can be defined for each material or phase by introducing the number of inclusion materials/phases presenting a constant eigenstrain via parameter NUMBER on the keyword line. Then on the next lines, the user has to specify, per inclusion, the name of the corresponding phase via the parameter PHASE and the constant value of the isotropic eigenstrain. Ex. 19 illustrates the definition of constant eigenstrain per inclusion for a steel alloy where the austenite phase is the matrix and ferrite and pearlite are considered as inclusions.

    #
    /EIGENSTRAIN, NUMBER=2
    PHASE=FERRITE, 4.3135E-03
    PHASE=PEARLITE, 3.2137E-03
    #
    

    Example 19: Definition of constant eigenstrain per inclusion phase for a steel alloy.

    Note: This definition neglects, at the moment, the anisotropic character of the eigenstrain of pearlite composed of lamellae of ferrite and cementite.

• /LOCA, NUMBER = number

  This option defines the number of critical locations in the sample at the macro-scale, for which a localization run will be performed on the RVE. On the keyword line the user has to specify the number of desired localization analysis via parameter NUMBER. Then, on the next lines and, for each localization analysis, the user specifies the element number of the macro-model via the parameter ELE and the Gauss point number via parameter GP. If GP=0, then all Gauss points of the element are considered. Then, depending of the element type at the macro-scale, defined in the transfer file (described in Sec. Transfer of macro-scale results from ABAQUS to HOMAT), NGP localization analyses will be performed. Moreover, if the parameter ELE=0, then all elements of the macro-simulation ELSET, defined in the transfer file via command (/INPUT FILES, MACRO=NAME_MACROSIM.txt) are selected. For example, ELE=0, GP=1, means that a localization run is performed for all elements of the transferred macro ELSET with the strain state taken always from the first Gauss point of each element.

  Ex. 20 illustrates the definition of 6 successive localization runs after O-forming of a pipeline tube (plane strain finite element analysis), which will be performed in the same HOMAT run.

  #
  /LOCA, NUMBER=3
  ELE=1297, GP=1
  ELE=1253, GP=0
  ELE=1248, GP=3
  #
  

  Example 20: Definition of elements and the Gauss points of the macro-model, where localization analyses will be performed.

  On the keyword line, parameter NUMBER gives the number of lines with specific localization input. For macro-elements 1297 (GP=1) and 1248 (GP=3) a localization run will be performed as well as for all Gauss points of macro-element 1253. As a plane strain forming simulation has been achieved at the macro-scale with linear quadrangular elements having 4 Gauss points, the definition of Ex. 20 will result in $1+4+1=6$ localization runs. If in the macro-simulation an ELSET=critical_elements is defined, which involves the elements 1297, 1253, 1248, then the following input leads to three localization runs: for each element of the ELSET the strain state in the first Gauss point is selected to be applied as boundary condition on the RVE of the steel microstructure:

  /LOCA, NUMBER=1
  ELE=0, GP=1
  

• /SPHERU CENTER, ANGSPH = angsph, NUMBER = number

  If the spherulite model (TYPE=SPHERU) is selected in the keyword /MATERIAL for several material domains, the user must then also specify the spherulite centres of the semi-crystalline microstructure. The microstructure evolution program SphaeroSim provides these centres to Mesh2HOMAT, which put them in a specific file, named PROBLEM_NAME.hsph, which can be read by HOMAT. On the keyword line, the mandatory parameter NUMBER has to be specified. It defines the total number of spherulite centres, specified in the RVE for the considered homogenization analysis. On the following lines and per spherulite, the user enters its identification number (ID), followed by the three coordinates of its centre. Note that the spherulite ID must corresponds to the number introduced in the /SOLID SECTION definition. A second, facultative parameter ANGSPH on the keyword line defines the rotation angle increment (in grad) between two successive axis locations in the plane perpendicular to the radial growth direction. This angle increment allows for a spherical spherulite to have an effective quasi-isotropic elastic behaviour. Numerical experiments on a single spherical spherulite have shown that this angle lies between 60 deg. and 90 deg. By default, ANGSPH=80 deg. More details about the definition of this angle increment can be found in ³.

  /SPHERU CENTER, ANGSPH=70, NUMBER=756
  1, 0.017, 0.009, 0.001
  ...
  736, 0.065,0.043, 0.099
  

  Example 21: The definition of spherulite centres for a solidified spherulite microstructure, extracted from a polypropylene plate.

  In this example, a RVE of extension 0.1 mm³ with 756 solidified spherulites is analysed. In the spherulite model an angle increment of 70° deg is adopted here. The coordinates of the spherulite centres are specified here in mm.

HOMAT Output¶

In the following section the different HOMAT result files and their nomenclature are introduced. Sec. Effective properties results concerns the description of the files containing the predicted effective properties of the considered heterogeneous material. In Sec. Field Results the nature and nomenclature of the VTK results files of a homogenization analysis are described in detail. These files contain field variables, which vary in the RVE of the heterogeneous material. In Sec. Localization Results the nature and nomenclature of the field output files of a localization analysis are outlined. Eventually, in Sec. Log File Output a brief overview of the principal control outputs in the log file, named either PROBLEM_NAME.hout for homogenization or PROBLEM_NAME.lout for localization is given.

Effective properties results¶

In the case of a monolayer or multilayer homogenization analysis, HOMAT provides always two results files: the first one is named PROBLEM_NAME.hres. It is an ASCII file containing the evaluated effective properties, presented in a documented and human-readable way. The second file, named PROBLEM_NAME.mat, contains the effective properties, written in the ABAQUS format for subsequent structural analysis on the macroscale.

Thermal homogenization results¶

For each selected layer or for the whole multilayer structure, HOMAT provides the calculated effective, orthotropic thermal conductivity tensor, expressed in the structural axis system and usually in W/(mm K). Moreover, the effective heat capacity [in J/(mm³ K)], density [in g/(mm³)] and specific heat [in J/(g K)], obtained by volume averaging (or homogenization of order zero), are printed in the result file .hres. Additionally, if desired by the user, HOMAT writes also the effective thermal conductivity tensor in the adopted local axis system. Ex. 22 shows the thermal homogenization result file .hres for a layer of a unidirectional (UD) composite with rectangular fibre shape, defined by Dvorak as a benchmark test ¹⁸ for his homogenization method. The rectangular fibre has a conductivity of 10^-2 W/(mm K); while the matrix has a larger conductivity of 1.0 W/(mm K). With his numerical approach, Dvorak predicts the following orthotropic effective thermal conductivities in the cross-section: $\lambda_{1} = 0.671$ W/(mm K) and $\lambda_{2} = 0.425$ W/(mm K). HOMAT provides with a coarse mesh of $12 \times 12 \times 2$ linear hexahedral elements the effective properties given in Ex. 22. The result file .hres starts always by giving the name of the present homogenization analysis, here hexas_scale, followed by the number of specified layers (here one) and the number of homogenizations, here also one. Then, the effective orthotropic thermal conductivity tensor is written in the structural axis system.

|                                                                          |
|     HOMOGENISATION RESULTS   -   EFFECTIVE  PROPERTIES                   |
|                                                                          |
|        SIMULATION:       hexas_scale_1                                   |
|                                                                          |
|         number of layers :          1                                    |
|         number of homogenisations :  1                                   |
|                                                                          |
|--------------------------------------------------------------------------|
|                       |                                                  |
|  heat conductivity    |    k_i1        |    k_i2        |    k_i3        |
|                       |                |                |                |
|          k_1j         |   6.776877E-01 |   1.803914E-12 |                |
|          k_2j         |   1.109493E-11 |   4.272344E-01 |                |
|                       |                |                |                |
|  mean value           |     Dk_x       |     Dk_y       |     Dk_z       |
|                       |                |                |                |
|      5.524610E-01     |      22.67     |     -22.67     |       0.00     |
|--------------------------------------------------------------------------|

|                       |                |                |                |
|  heat capacity        |    rhocp       |                |                |
|                       |                |                |                |
|                       |   6.500015E-02 |                |                |
|                       |                |                |                |
|  Density              |    dens        |                |                |
|                       |                |                |                |
|                       |   2.062501E-03 |                |                |
|                       |                |                |                |
|  specific heat        |    cp          |                |                |
|                       |                |                |                |
|                       |   2.750264E+01 |                |                |
|                       |                |                |                |
|--------------------------------------------------------------------------|

|--------------------------------------------------------------------------|
|  heat cond. vol. aver.|    kav_i1      |    kav_i2      |    kav_i3      |
|                       |                |                |                |
|          kav_1j       |   7.502492E-01 |   0.000000E+00 |                |
|          kav_2j       |   0.000000E+00 |   7.502492E-01 |                |

Figure 16: Cross-section of a UD reinforced composite with rectangular fibres (left) and Dvorak's triangular discretization (right).

Due to the rectangular shape of the fibres, the equivalent layer presents a pronounced orthotropic thermal behaviour with $\lambda_{1} = 0.67769$ W/(mm K) and $\lambda_{2} = 0.42723$ W/(mm K). In spite of the coarse mesh, these values differ from Dvorak's predictions only by about 1% for λ₁ and by 0.62% for $\lambda_{2}$. The next line provides the mean value and the relative deviation for each component from this mean value. A deviation of ± 22.67 % from the mean, isotropic conductivity is computed, which emphasizes the significant orthotropic character of the equivalent material. Effective heat capacity, density and specific heat are also provided. As the output parameter VOLAVRG=YES is set, the volume averaged thermal conductivity is written in the output file. As the volume averaging procedure takes only the volume fraction of the fibre and not its rectangular shape into account, it predicts an erroneous isotropic equivalent behaviour! This simple example shows the advantage of a first order asymptotic homogenization scheme in comparison to the volume averaging method. Ex. 23 shows the contents of the .mat file for a further thermal analysis at the macroscale with ABAQUS. As no name for the investigated composite is provided, HOMAT denotes the equivalent layer simply as a RVE material. The user can change this default name, by editing the .mat file.

,\*,\*------------------------------
,\*,\*
,\*,\* NUMBER OF SPECIFIED MATERIALS:      1
,\*,\*
,\*MATERIAL, name = RVE
,\*,\*
,\*DENSITY
  2.062501E-03
,\*CONDUCTIVITY, TYPE=ORTHO
  6.776877E-01,  4.272344E-01,  5.524610E-01
,\*,\*
,\*SPECIFIC HEAT
  2.750264E+01

Thermal and thermoelastic homogenization results¶

In order to illustrate the versatility of HOMAT, we predict in the same homogenization run effective thermal and mechanical properties of a metal matrix composite at different selected temperatures. The example of a FG75 NiAl matrix reinforced by UD sapphire fibres of cylindrical shape is considered here. It is assumed that the fibres have a volume content of 30%. The ASCII output result file .hres starts, after the definition of the problem name, with the print of the thermal homogenization results, followed by the thermo-elastic ones. Example 24 shows the thermal part of the .hres result file.

|                                                                          |
|     HOMOGENISATION RESULTS   -   EFFECTIVE  PROPERTIES                   |
|                                                                          |
|        SIMULATION:       fas30_thmec_ani                                 |
|                                                                          |
|                       SIMPLE PROBLEM DEFINITION                          |
|                                                                          |
|         number of layers :          1                                    |
|         number of homogenisations :  2                                   |
|                                                                          |
|--------------------------------------------------------------------------|
|                       |                                                  |
|  heat conductivity    |    k_i1        |    k_i2        |    k_i3        |
|                       |                |                |                |
|          800.00                                                          |
|          k_1j         |   3.434308E-02 |  -6.504486E-10 |   0.000000E+00 |
|          k_2j         |   1.710552E-11 |   3.434307E-02 |   0.000000E+00 |
|          k_3j         |  -6.310036E-11 |  -4.336384E-11 |   3.937139E-02 |
|                       |                |                |                |
|         1100.00                                                          |
|          k_1j         |   3.133230E-02 |  -8.107663E-10 |   0.000000E+00 |
|          k_2j         |  -5.581284E-11 |   3.133228E-02 |   0.000000E+00 |
|          k_3j         |  -7.202623E-11 |  -5.052411E-11 |   3.769298E-02 |
|                       |                |                |                |
|                       |                |                |                |
|  mean value           |     Dk_x       |     Dk_y       |     Dk_z       |
|                       |                |                |                |
|      3.601918E-02     |      -4.65     |      -4.65     |       9.31     |
|      3.345252E-02     |      -6.34     |      -6.34     |      12.68     |
|--------------------------------------------------------------------------|

|                       |                |                |                |
|  heat capacity        |    rhocp       |                |                |
|                       |                |                |                |
|          800.00       |   3.076047E-03 |                |                |
|         1100.00       |   3.532198E-03 |                |                |
|                       |                |                |                |
|  Density              |    dens        |                |                |
|                       |                |                |                |
|          800.00       |   4.742357E-03 |                |                |
|         1100.00       |   4.727331E-03 |                |                |
|                       |                |                |                |
|  specific heat        |    cp          |                |                |
|                       |                |                |                |
|          800.00       |   6.559209E-01 |                |                |
|         1100.00       |   7.480665E-01 |                |                |
|                       |                |                |                |

For each thermophysical quantity, effective orthotropic or isotropic properties are denoted at both defined temperatures, here at 800°C and 1100°C respectively. The cylindrical shape (see Fig. 4) induces an orthotropic effective behaviour. Indeed, deviations from the mean thermal conductivity of -4.65% ($\lambda_{11}$ and $\lambda_{22}$) and of 9.31% for $\lambda_{33}$ in fiber direction are predicted respectively at 800°C. Slightly higher deviations are predicted at 1100°C. The mechanical result part (see Ex. 25) starts by providing the calculated effective Hooke matrix of the UD reinforced composite, before its symmetrization. This way, the user can easily check how far the calculated Hooke matrix is from the expected symmetric one and obtains an error measure for the impact of the numerical discretization. As the sapphire fiber is monoclinic with only one symmetry plane perpendicular to the fibre, the effective Hooke matrix of the UD reinforced composite is fully anisotropic with extension-shear coupling terms. In order to derive effective engineering properties, the calculated Hooke matrix is symmetrized and inverted. The deduced effective shear and Young moduli as well as the Poisson coefficients of the composite layer are given in Example 26.

|   Hooke matrix 6*6    |    H_i1        |    H_i2        |    H_i3        |
|          800.00                                                          |
|          H_1j         |   3.001139E+05 |   1.328843E+05 |   1.229083E+05 |
|          H_2j         |   1.328979E+05 |   3.000802E+05 |   1.228837E+05 |
|          H_3j         |   1.195747E+05 |   1.195287E+05 |   3.190173E+05 |
|          H_4j         |  -3.502496E-03 |   2.038315E-03 |   1.421406E-04 |
|          H_5j         |  -2.420831E-04 |   2.026955E-04 |  -2.636610E-08 |
|          H_6j         |   2.955201E+03 |  -2.814527E+03 |   4.557293E-02 |
|                       |                |                |                |
|          part II      |    H_i4        |    H_i5        |    H_i6        |
|                       |                |                |                |
|          H_1j         |  -3.836078E-03 |  -7.183098E-04 |   2.301368E+03 |
|          H_2j         |   1.579003E-03 |   3.476790E-04 |  -2.301262E+03 |
|          H_3j         |   1.188207E-04 |  -1.084709E-05 |   1.058137E-01 |
|          H_4j         |   8.040503E+04 |   1.258861E+03 |  -5.246791E-04 |
|          H_5j         |   2.525469E+03 |   8.446847E+04 |  -5.621908E-04 |
|          H_6j         |  -1.232085E-04 |  -3.716822E-04 |   8.550840E+04 |

|  Shear moduli         |     G_xy       |     G_xz       |    G_yz        |
|                       |                |                |                |
|          800.00       |   8.036740E+04 |   8.442893E+04 |   8.542898E+04 |
|         1100.00       |   7.131272E+04 |   7.869988E+04 |   7.938645E+04 |
|                       |                |                |                |
|  mean value           |    DG_xy       |    DG_xz       |    DG_yz       |
|                       |                |                |                |
|      8.340844E+04     |       -3.65    |        1.22    |        2.42    |
|      7.646635E+04     |       -6.74    |        2.92    |        3.82    |
|--------------------------------------------------------------------------|
|                       |                |                |                |
|  Young moduli         |     E_x        |     E_y        |    E_z         |
|                       |                |                |                |
|          800.00       |   2.242213E+05 |   2.242187E+05 |   2.511521E+05 |
|         1100.00       |   2.018505E+05 |   2.018355E+05 |   2.332937E+05 |
|                       |                |                |                |
|  mean value           |    DE_x        |    DE_y        |    DE_z        |
|                       |                |                |                |
|      2.331973E+05     |       -3.85    |       -3.85    |        7.70    |
|      2.123265E+05     |       -4.93    |       -4.94    |        9.87    |
|--------------------------------------------------------------------------|
|                       |                |                |                |
|  Poisson coefficients |    nu_xy       |    nu_xz       |    nu_yz       |
|                       |                |                |                |
|          800.00       |   3.422293E-01 |   2.800512E-01 |   2.798895E-01 |
|         1100.00       |   3.459106E-01 |   2.791244E-01 |   2.791171E-01 |
|                       |                |                |                |
|  mean value           |    Dnu_xy      |    Dnu_xz      |    Dnu_yz      |
|                       |                |                |                |
|      3.007234E-01     |       13.80    |       -6.87    |       -6.93    |
|      3.013840E-01     |       14.77    |       -7.39    |       -7.39    |
|--------------------------------------------------------------------------|
|                       |                |                |                |
|  thermal expansion    |    alpha_xx    |    alpha_yy    |    alpha_zz    |
|                       |                |                |                |
|          800.00       |   1.265114E-05 |   1.240363E-05 |   1.100857E-05 |
|         1100.00       |   1.430439E-05 |   1.403192E-05 |   1.200152E-05 |
|                       |                |                |                |
|  mean value           |    Dal_xx      |    Dal_yy      |    Dal_zz      |
|                       |                |                |                |
|      1.202111E-05     |        5.24    |        3.18    |       -8.42    |
|      1.344594E-05     |        6.38    |        4.36    |      -10.74    |
|--------------------------------------------------------------------------|
|                       |                |                |                |
|  shear-extension coeff|    eta_i,12    |    eta_i,13    |    eta_i,23    |
|                       |                |                |                |
|          800.00                                                          |
|          eta_1j       |   5.340722E-08 |   5.363152E-09 |  -4.097442E-02 |
|          eta_2j       |  -3.741822E-08 |  -4.137568E-09 |   4.043315E-02 |
|          eta_3j       |  -8.072154E-09 |  -4.344513E-10 |   2.375421E-04 |
|         1100.00                                                          |
|          eta_1j       |   6.356865E-08 |   9.550613E-09 |  -4.018362E-02 |
|          eta_2j       |  -4.509562E-08 |   8.369636E-10 |   3.952047E-02 |
|          eta_3j       |  -7.773691E-09 |  -4.573540E-08 |   2.822128E-04 |

As expected, the equivalent composite material presents, in a first approximation, a transverse isotropic behaviour with a pronounced stiffness increase in fibre direction. Otherwise, the effective mechanical properties decrease with the temperature. As a thermoelastic homogenization analysis has been specified, HOMAT provides the effective thermal expansion coefficients in the structural directions. Due to the full anisotropy of the Hooke matrix, effective shear-extension coefficients ¹⁴, named $\eta_{ij}$ with $i \neq j$, are given in the .hres result file.

• In order to underline that the effective mechanical properties are expressed in the structural RVE axis system (x, y, z), each engineering property (e.g. Young, shear, Poisson, thermal expansion and eigenstrain) is doted by subscripts of type xx, xy, yy, etc.
• If an isotropic eigenstrain is specified per phase, like in solid phase transformation simulations with MICRESS, the generalized effective volumetric eigenstrain tensor $\kappa_{ij}$ is provided in addition to the thermal expansion coefficients. In the present analysis VOLAVRG is set to NO, otherwise the volume averaged expressions of all the thermal and mechanical properties would have been provided. The corresponding PROBLEM_NAME.mat file for a further thermal and/or thermoelastic analysis at the macroscale with Abaqus is given in Ex. 27. The symmetrised Hooke matrix is provided via the material keyword *ELASTIC. Please note also the large analogy between the present Abaqus definition and the HOMAT material definition in Sec. Definition of the Materials Properties. As the special output option ENG=YES is activated, HOMAT writes also the effective orthotropic Young and shear moduli and Poisson coefficients for the homogenized material. The anisotropic Hooke matrix suffices ABAQUS. Thus, the effective engineering constants are not read by ABAQUS but might be useful for later post-processing.

,\*,\*------------------------------
,\*,\*
,\*,\* NUMBER OF SPECIFIED MATERIALS:      1
,\*,\*
,\*MATERIAL, name = RVE
,\*,\*
,\*DENSITY
  4.742357E-03,     800.00
  4.727331E-03,    1100.00
,\*CONDUCTIVITY, TYPE=ORTHO
  3.434308E-02,  3.434307E-02,  3.937139E-02,     800.00
  3.133230E-02,  3.133228E-02,  3.769298E-02,    1100.00
,\*,\*
,\*SPECIFIC HEAT
  6.559209E-01,     800.00
  7.480665E-01,    1100.00
,\*ELASTIC, TYPE=ANISOTROPIC
  3.001139E+05,  1.328911E+05,  3.000802E+05,  1.212415E+05,  1.212062E+05,  3.190173E+05, -3.669287E-03,  1.808659E-03,
  1.304807E-04,  8.040503E+04, -4.801965E-04,  2.751873E-04, -5.436726E-06,  1.892165E+03,  8.446847E+04,  2.628284E+03,
 -2.557895E+03,  7.569334E-02, -3.239438E-04, -4.669365E-04,  8.550840E+04,     800.00
  2.694108E+05,  1.192916E+05,  2.693842E+05,  1.084949E+05,  1.084875E+05,  2.938578E+05, -3.875977E-03,  1.909514E-03,
  1.194593E-04,  7.134687E+04,  5.284290E-04,  1.202110E-03,  4.097078E-03,  1.782636E+03,  7.873757E+04,  2.393042E+03,
 -2.312515E+03,  9.735707E-02, -1.753401E-04, -2.586395E-04,  7.945887E+04,    1100.00
,\*,\*
,\*,\* 2. orthotropic variant neglecting shear-extension coupling
,\*,\*
,\*ELASTIC, TYPE=ENGINEERING CONSTANTS
  2.242213E+05,  2.242187E+05,  2.511521E+05,  3.422293E-01,  2.800512E-01,  2.798895E-01,  8.036740E+04,  8.442893E+04,
  8.542898E+04,     800.00
  2.018505E+05,  2.018355E+05,  2.332937E+05,  3.459106E-01,  2.791244E-01,  2.791171E-01,  7.131272E+04,  7.869988E+04,
  7.938645E+04,    1100.00
,\*,\*
,\*YOUNG MODULUS, TYPE=ORTHO
  2.242213E+05,  2.242187E+05,  2.511521E+05,     800.00
  2.018505E+05,  2.018355E+05,  2.332937E+05,    1100.00
,\*,\*
,\*SHEAR MODULUS, TYPE=ORTHO
  8.036740E+04,  8.442893E+04,  8.542898E+04,     800.00
  7.131272E+04,  7.869988E+04,  7.938645E+04,    1100.00
,\*,\*
,\*POISSON RATIO, TYPE=ORTHO
  3.422293E-01,  2.800512E-01,  2.798895E-01,     800.00
  3.459106E-01,  2.791244E-01,  2.791171E-01,    1100.00
,\*,\*
,\*EXPANSION, TYPE=ORTHO
  1.265114E-05,  1.240363E-05,  1.100857E-05,     800.00
  1.430439E-05,  1.403192E-05,  1.200152E-05,    1100.00
,\*,\*

2-D homogenization results¶

As the effective material can be either isotropic, orthotropic, cubic or generally anisotropic and the homogenization analysis can be of 3-D or 2-D nature, there exists several different output formats of the ASCII .hres file. In the present manual, only two 2-D variants: results of plane strain or plane stress homogenization, are briefly outlined.

• Plane strain homogenization

  The plane strain assumption ($\epsilon_{33}$ =0) leads to following strain and stress tensors:

  $$\begin{equation} \label{eq:str_planestrain} \epsilon^{\text{t}} = \left( \epsilon_{\text{xx}},\epsilon_{\text{yy}},\ 0,\ \gamma_{\text{xy}} \right)\,,\ \ \ \ \ \ \sigma^{\text{t}} = (\sigma_{\text{xx}},\ \sigma_{\text{yy}},\ \sigma_{\text{zz}},\ \tau_{\text{xy}})\,. \end{equation}$$

  The corresponding symmetric Hooke matrix, in $4\times 4$ notation, has the following expression:

  $$\begin{equation} \label{eq:H__planestrain} \mathbb{H} = \begin{bmatrix} H_{11} & H_{12} & H_{13} & H_{14} \\ \bullet & H_{22} & H_{23} & H_{24} \\ \bullet & \bullet & H_{33} & H_{34} \\ \bullet & \bullet & \bullet & H_{44} \\ \end{bmatrix}\,. \end{equation}$$

  The output format in the .hres file takes the reduced size $(4\times 4)$ of the Hooke matrix into account, as illustrated by Ex. 28 of the UD reinforced glass-epoxy composite, whose properties are given in Sec. Definition of the Materials Properties. As the out-plane shears $\tau_{13}$ and $\tau_{23}$ are both zero, only the shear modulus in the cross-section plane 1-2 (or x-y) is printed in the .hres file.

  ----------------------------------------------------------------------------
  |                                                                          |
  |     HOMOGENISATION RESULTS   -   EFFECTIVE  PROPERTIES                   |
  |                                                                          |
  |        SIMULATION:       per3D_scal05_evz                                |
  |                                                                          |
  |         number of layers :          1                                    |
  |         number of homogenisations :  1                                   |
  |                                                                          |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |   Hooke matrix 4*4    |    H_i1        |    H_i2        |    H_i3        |
  |                       |                |                |                |
  |          H_1j         |   1.443622E+04 |   4.898715E+03 |   5.274505E+03 |
  |          H_2j         |   4.898715E+03 |   1.443606E+04 |   5.274471E+03 |
  |          H_3j         |   5.274505E+03 |   5.274471E+03 |   4.626249E+04 |
  |          H_4j         |   1.344940E-08 |  -4.578791E-07 |  -5.370989E-08 |
  |                       |                |                |                |
  |          part II      |    H_i4        |                |                |
  |                       |                |                |                |
  |          H_14         |  -1.445253E-08 |                |                |
  |          H_24         |  -1.148775E-08 |                |                |
  |          H_34         |  -5.451318E-09 |                |                |
  |          H_44         |   3.065298E+03 |                |                |
  |                       |                |                |                |
  |                       |                |                |                |
  |  Shear moduli         |     G_xy       |                |                |
  |                       |                |                |                |
  |                       |   3.065298E+03 |                |                |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  Young moduli         |     E_x        |     E_y        |    E_z         |
  |                       |                |                |                |
  |                       |   1.250001E+04 |   1.249987E+04 |   4.338476E+04 |
  |                       |                |                |                |
  |  mean value           |    DE_x        |    DE_y        |    DE_z        |
  |                       |                |                |                |
  |      2.279488E+04     |      -45.16    |      -45.16    |       90.33    |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  Poisson coefficients |    nu_xy       |    nu_xz       |    nu_yz       |
  |                       |                |                |                |
  |                       |   3.106216E-01 |   2.727963E-01 |   2.727974E-01 |
  |                       |                |                |                |
  |  mean value           |    Dnu_xy      |    Dnu_xz      |    Dnu_yz      |
  |                       |                |                |                |
  

  Example 28: Plane strain effective elastic properties of an UD reinforced glass/epoxy composite.

• Plane stress homogenization

  In a similar way, the plane stress assumption leads to the following stress and strain tensors:

  $$\begin{equation} \label{eq:str_planestress} \epsilon^{\text{t}} = \left( \epsilon_{\text{xx}},\epsilon_{\text{yy}},\ \ \gamma_{\text{xy}} \right)\,;\qquad\sigma^{\text{t}} = (\sigma_{\text{xx}},\ \sigma_{\text{yy}},\ \ \tau_{\text{xy}}) \end{equation}$$

  The corresponding anisotropic Hooke matrix, expressed in Voigt notation, has following expression:

  $$\begin{equation} \label{eq:H_planestress} \mathbb{H} = \begin{bmatrix} H_{11} & H_{12} & H_{13} \\ \bullet & H_{22} & H_{23} \\ \bullet & \bullet & H_{33} \\ \end{bmatrix}\,. \end{equation}$$

  The output in the .hres file takes the reduced size of the Hooke matrix $(3\times 3)$ into account, as illustrated by Ex. 29 of the UD reinforced glass-epoxy composite.

  |  Hooke matrix 3*3     |    H_i1        |    H_i2        |    H_i3        |
  |                       |                |                |                |
  |          H_1j         |   1.027775E+04 |   6.249302E+02 |   0.000000E+00 |
  |          H_2j         |   6.248244E+02 |   1.027799E+04 |   0.000000E+00 |
  |          H_3j         |   0.000000E+00 |   0.000000E+00 |   3.932720E+03 |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  Young/shear moduli   |     E_x        |     E_y        |    G_xy        |
  |                       |                |                |                |
  |                       |   1.023976E+04 |   1.024000E+04 |   3.932720E+03 |
  |                       |                |                |                |
  |  mean value           |    DE_x        |    DE_y        |                |
  |                       |                |                |                |
  |      1.023988E+04     |       -0.00    |        0.00    |                |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  Elastic moduli       |     C_11       |     C_12       |    C_44        |
  |                       |                |                |                |
  |                       |   1.027787E+04 |   6.248773E+02 |   3.932720E+03 |
  |                       |                |                |                |
  |  Eng. moduli          |     E          |    nu          |    Zener coeff.|
  |                       |                |                |                |
  |                       |   1.023988E+04 |   6.079763E-02 |   8.148189E-01 |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  Poisson coefficients |    nu_xy       |    nu_yx       |                |
  |                       |                |                |                |
  |                       |   6.079763E-02 |   6.079903E-02 |                |
  

  Example 29: Plane stress effective elastic properties of the UD reinforced epoxy/glass composite.

  As in this analysis the parameter CUBVAL of the keyword /PARA is set to 0.1, the transverse quasi-isotropy of the UD reinforced composite leads HOMAT to provide also the effective elastic properties of the equivalent cubic material. Indeed, the three elastic constants $C_{11}$, $C_{12}$ and $C_{44}$ are given as well as the mean Young modulus, the mean Poisson coefficient and the Zener anisotropy coefficient, defined by

  $$\begin{equation*} Z = \frac{2C_{44}}{C_{11} - \ C_{12}}\,. \end{equation*}$$

  Moreover, in this run VOLAVRG=YES, therefore the volume averaged Hooke matrix is also printed. For the thermal expansion coefficient, not only the volume averaged expression is evaluated via the thermal stresses and printed but also the value obtained by the direct application of the mixture rule to the thermal expansion coefficient without using the thermal stress expression. The corresponding results are shown in Example 30 below.

  |  th. exp. vol. aver.  |    alpha_11    |    alpha_22    |                |
  |                       |                |                |                |
  |                       |   4.406312E-06 |   4.406312E-06 |                |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  th. expan. mixture   |    alpha_11    |    alpha_22    |                |
  |                       |                |                |                |
  |                       |   9.587308E-06 |   9.587308E-06 |                |
  

  Example 30: 2-D volume averaged and mixture thermal expansion coefficients of the UD reinforced epoxy/glass composite.

  Note: The different output of both 2-D homogenization analyses with some different parameters (CUBVAL, VOLAVRG) shows the great variety and the complexity of the output file .hres.

Results in a local axis system¶

Since version 6.0 the user has the possibility to specify a local axis system in which the calculated effective Hooke matrix is printed in the PROBLEM_NAME.hres and PROBLEM_NAME.mat files. To illustrate this new output a more industrial example of a hollow wireframe sandwich, built of two Nicrofer plates which surround woven firing and chaining wires (see Fig. 17), is adopted here. This hollow wire frame is placed on the stator of a steam turbine and it is cooled by steam at 425°C. The thermal and mechanical homogenizations are performed at 650°C. Due to the ECAP drawing process, the wires present elongated grains in the extrusion direction. Therefore, their mechanical properties are orthotropic. For each wire a local axis system in which the material properties are specified via the coordinates of two nodes (see Ex. 31), is used. In this example, the local axis system of the chain wire is chosen via parameter OUT_ORIE on the /ORIENTATION keyword for additional output of the effective anisotropic Hooke matrix of the cooled hollow multilayer. Note that this additional output requires the setting of LOC_AXIS=YES as parameter of the /OUTPUT keyword.

Figure 17: Mesh of the hollow sandwich built by two Nicrofer plates (only the bottom one is drawn) and the woven wireframe, where the chain wires are located along _y_ axis, whereas firing wires are along the _x_ axis.

#
# Output definition
#
/OUTPUT, FORMAT=VTK, TYPE=ASCII, IOP=2, VOLAVRG=YES, VISLAV=YES,
LOC_AXIS=YES, SHEAR_TENS=YES
U
/ORIENTATION, NAME=wire_X_axis, DEFINITION=COORDINATES, OUT_ORIE
0., 2.0, 0., 0., 0., 2.0
/ORIENTATION, NAME=wire_Y_axis, DEFINITION=COORDINATES
1.0, 0.0, 0.0, 0.0, 0.0, 1.0
#
# Definition of the solid sections
# per solid section, specify the desired analysis type
#
/SOLID SECTION, NUMBER=1, ELSET=SCHUSSDRAHT, NAME=firing_wire,
MATERIAL=HT-6025,
ORIENTATION=wire_Y_axis, TEMP=NO, MECH=HOM
#
/SOLID SECTION, NUMBER=2, ELSET=KETTDRAHT, NAME=chain_wire,
MATERIAL=HT-6025,
ORIENTATION=wire_X_axis, TEMP=NO, MECH=HOM

In fact, the local axis 1 is directed along structural axis Y; whereas local axis 2 corresponds to the structural direction Z and the 3^d orthogonal axis is along structural axis X. Therefore, the local Young, shear moduli, listed in Ex. 32 correspond to a simple permutation of the structural values. This permutation can easily be verified by the user.

• Output in the structural axis system:

  |  Shear moduli         |     G_xy       |     G_xz       |    G_yz        |
  |                       |                |                |                |
  |                       |   4.255631E+04 |   6.807693E+04 |   6.807691E+04 |
  |                       |                |                |                |
  |  mean value           |    DG_xy       |    DG_xz       |    DG_yz       |
  |                       |                |                |                |
  |      5.957005E+04     |      -28.56    |       14.28    |       14.28    |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  Young moduli         |     E_x        |     E_y        |    E_z         |
  |                       |                |                |                |
  |                       |   1.174384E+05 |   1.176799E+05 |   1.607221E+05 |
  |                       |                |                |                |
  |  mean value           |    DE_x        |    DE_y        |    DE_z        |
  |                       |                |                |                |
  |      1.319468E+05     |      -11.00    |      -10.81    |       21.81    |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  Poisson coefficients |    nu_xy       |    nu_xz       |    nu_yz       |
  |                       |                |                |                |
  |                       |   2.787595E-01 |   4.118723E-01 |   3.931044E-01 |
  |                       |                |                |                |
  |  mean value           |    Dnu_xy      |    Dnu_xz      |    Dnu_yz      |
  |                       |                |                |                |
  |      3.612454E-01     |      -22.83    |       14.01    |        8.82    |
  

  Example 32: Effective engineering moduli of the hollow multilayer expressed in the structural axis system.
  - Output in the local axis system

  |  local shear moduli   |     G_12       |     G_13       |    G_23        |
  |                       |                |                |                |
  |                       |   6.807691E+04 |   4.255631E+04 |   6.807693E+04 |
  |                       |                |                |                |
  |  mean value           |    DG_12       |    DG_13       |    DG_23       |
  |                       |                |                |                |
  |      5.957005E+04     |       14.28    |      -28.56    |       14.28    |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  local Young moduli   |     E_1        |     E_2        |    E_3         |
  |                       |                |                |                |
  |                       |   1.176068E+05 |   1.605981E+05 |   1.174333E+05 |
  |                       |                |                |                |
  |  mean value           |    DE_1        |    DE_2        |    DE_3        |
  |                       |                |                |                |
  |      1.318794E+05     |      -10.82    |       21.78    |      -10.95    |
  |--------------------------------------------------------------------------|
  |                       |                |                |                |
  |  local Poisson coeff. |    nu_12       |    nu_13       |    nu_23       |
  |                       |                |                |                |
  |                       |   2.879694E-01 |   2.789950E-01 |   4.116669E-01 |
  |                       |                |                |                |
  |  mean value           |    Dnu_12      |    Dnu_13      |    Dnu_23      |
  |                       |                |                |                |
  |      3.262104E-01     |      -11.72    |      -14.47    |       26.20    |
  

  Example 33: Effective engineering moduli of the hollow multilayer expressed in the local axis system.

Stokes flow homogenization input and results¶

For each specified layer or classically for the whole multilayer, HOMAT provides the calculated Darcy permeability tensor in the structural axis system of the RVE and usually in mm^-2. Additionally, the permeability tensor is also written in the local axis system, if such a local axis system has been specified via parameter LOC_AXIS on the /OUTPUT keyword and via parameter OUT_ORIE to select its orientation. As the definition of a Stokes flow homogenization problem is new since version 6.0, not only the permeability results are outlined for a benchmark test, but also the specific commands of the problem definition. The flow around and parallel to a cylindrical fiber of radius 2 mm in a cubic box of dimension $10 \times 10 \times 10$ mm³ is taken here as example. Ex. 34 outlines the specific commands for Stokes flow homogenization. Note that at least one domain has to be defined as fluid domain via following allowable solid section names: FLUID, AIR, GAS, MELT, LIQUID, POROUS or by following German names LUFT and SCHMELZE. Note: No other names are accepted.

#
/PROBLEM, TYPE=STOKES, HOM_MOD=0, PER_MOD=3
#
/OUTPUT, FORMAT=VTK, TYPE=ASCII, IOP=2, VOLAVRG=YES, ENG=NO, VISLAV=YES, LOC_AXIS=NO
V, P
#
/INCLUDE, NAME=FibBL_r2_bubLi.hgeo.gz
#
MATERIAL, NAME=fiber, TYPE=ISO
/INCLUDE, NAME=silica.hmat
#
/MATERIAL, NAME=AIR, TYPE=ISO
/INCLUDE, NAME=luft_mm.hmat
#
/SOLID SECTION, NUMBER=1, ELSET=Set-Matrix, NAME=fluid, MATERIAL=AIR,
TEMP=NO, MECH=NO, FLOW=CALC
#
/SOLID SECTION, NUMBER=2, ELSET=Set-Fiber, NAME=FIBER, MATERIAL=fiber
TEMP=NO, MECH=NO, FLOW=NO
#

In order to respect the general definition of the /SOLID SECTION keyword, also the material has to be specified, nevertheless the defined properties (density, viscosity, thermal conductivities, ...) are not used by HOMAT. HOMAT solves three (in 3-D) or two (in 2-D) special Stokes flow problems, where the density and the viscosity are both taken as unitary and a unitary volumetric boost force is applied in the structural directions X, Y and Z, respectively. Independent of the user specification FLOW=YES or NO for a solid domain, HOMAT exclude the solid domains from the reduced system, except the interface nodes, where no-slip boundary conditions are always applied. In Example 35 the predicted Darcy permeability tensor for the cylindrical fiber RVE is outlined. HOMAT 6.0 provides the principal permeabilities and the corresponding eigenvectors. Moreover, it computes the R anisotropy coefficient defined by the principal permeabilities $K_{1} > K_{2} > K_{3}$ via the expression ¹⁹

----------------------------------------------------------------------------
|                                                                          |
|     HOMOGENISATION RESULTS   -   EFFECTIVE  PROPERTIES                   |
|                                                                          |
|        SIMULATION:       FibBL_r2_bubLi                                  |
|        porosity:            87.45 %                                      |
|                                                                          |
|         number of layers :          1                                    |
|         number of homogenisations :  1                                   |
|                                                                          |
|--------------------------------------------------------------------------|
|                       |                |                |                |
|  permeability tensor  |    K_X         |    K_Y         |    K_Z         |
|                       |                |                |                |
|          comp. X      |   5.399124E+00 |   4.867793E-12 |  -1.735329E-15 |
|          comp. Y      |  -4.123367E-13 |   5.399125E+00 |  -6.511195E-16 |
|          comp. Z      |  -3.288843E-10 |   2.218531E-10 |   6.668920E+00 |
|                       |                |                |                |
|--------------------------------------------------------------------------|
|                       |                |                |                |
|  princ. perm. R_ani   |    K_1         |    K_2         |    K_3         |
|                       |                |                |                |
|    8.99775E-01        |   6.668920E+00 |   5.399125E+00 |   5.399124E+00 |
|--------------------------------------------------------------------------|
|                       |                |                |                |
|  unit vector - K_1    |    E1_1        |    E2_1        |    E3_1        |
|                       |                |                |                |
|                       |  -1.295035E-10 |   8.735858E-11 |   1.000000E+00 |
|                       |                |                |                |
|  unit vector - K_2    |    E1_2        |    E2_2        |    E3_2        |
|                       |                |                |                |
|                       |  -2.306980E-06 |  -1.000000E+00 |   8.735806E-11 |
|                       |                |                |                |
|  unit vector - K_3    |    E1_3        |    E2_3        |    E3_3        |
|                       |                |                |                |
|                       |   1.000000E+00 |  -2.306980E-06 |   1.295037E-10 |

Note that HOMAT also provides the porosity of the RVE as additional information. The predicted effective Darcy permeabilities are, as expected, transverse isotropic and the off-diagonal terms are nearly zero. The permeability parallel to the fibers, denoted $K_{\parallel}$, is larger than the transverse one due to less flow resistance parallel to the fibers, characterized by a lower tortuosity. These predictions are achieved with a coarse hexahedral mesh of 37120 elements and 41517 nodes but with a boundary layer around the fiber (see Figs. 20 and 21 of the pressure distribution around the central fiber). The permeability $K_{\parallel}$ is in excellent agreement with the analytical solution of Tamayol and Bahrami ²⁰ for the flow parallel to a square arrangement of cylindrical fibers: $K_{\parallel}^{a} = 6.6731$ mm^-2;; whereas the transverse permeability prediction $P_{\perp}= 5.3991$ mm^-2 overestimates the model value $K_{\perp}^{a}= 3.696$ mm² of Tamayol and Bahrami ²¹, leading to a less pronounced orthotropic permeability behaviour, in agreement with experimental results at high porosities (here $\epsilon = 0.874$).

Note: If the parameter LOC_AXIS=YES is specified, the permeability tensor is also printed in the local axis system (1, 2, 3). Compared to the output of Ex. 35, the descriptor becomes now local perm. tensor instead permeability tensor and the permeability components are noted now K_1, K_2 and K_3 instead K_X, K_Y and K_Z, indicating thus the output in the local axis system (1, 2, 3).

Field Results¶

In this section the previous example of the sapphire fibre in a metal matrix is reconsidered to detail the content of the different field output files in the case of a combined thermal and thermoelastic homogenization. Field results are characterized by nodal values, defined for each considered layer. They describe the local variation of the field on the layer or RVE and are written in the VTK format. They can be visualized via the open-source program Paraview.

Thermal field results¶

If the parameter T or ALL is specified in the /OUTPUT keyword, each thermal mono- or multilayer homogenization produces two or three different VTK files.

• The 1^st file contains all the generated implicit fluxes at the domain interfaces, where a jump of material properties exists. If the thermal conductivity varies inside each domain with the temperature, it contains also their volumetric contribution. The number of flux terms depends on the homogenization type: either 3 in the 3-D case or 2 in the 2-D case. The corresponding keyword in the VTK file is IMPL_FLUX. Since version 5.2, the implicit fluxes are saved as vector solutions in the VTK file. Its nomenclature is the following: to the problem name, fas30_fine, is added extension _flux. Then, two numbers are added: the 1^st one specifies the considered layer and the second one, the considered homogenization run. Each set temperature on the RVE (keyword /TEMP DOMAIN) leads to a separate homogenization analysis. In the present case, as two temperatures are specified ($T=800$°C and $T=1100$°C) on one layer, HOMAT writes following two implicit flux files for visualization purpose: fas30_fine_flux11.vtk and fas30_fine_flux12.vtk: the 1^st corresponds to the thermal homogenization at 800°C and the 2^d one to the homogenization at 1100°C.

• The second file contains the periodic microscopic displacements $x^{r}$ , with $r=1, 3$ in the 3-D case or $r=1, 2$ in the 2-D case on the undeformed RVE. Each displacement component is the solution of a linearized thermal system. The corresponding keyword in the VTK file for this vector solution is: DISPL_TH. Its nomenclature is similar to the one for implicit fluxes: HOMAT adds to the problem name the extension _hola followed by two numbers: the layer and homogenization indicators. Thus, for the metal matrix composite example, HOMAT generates the following two files with microscopic displacements: fas30_fine_hola11.vtk and fas30_fine_hola12.vtk. Note that these vtk files have a common part composed by a header, the definition of the undeformed finite element mesh and the local nodal values of the material_ID parameter. This parameter allows the drawing of the different domains or grains of the RVE. In Fig. 18 (left) the implicit flux in X direction is drawn for the UD reinforced metal matrix composite and in Fig. 18 (right) the corresponding microscopic displacement. Due to the higher thermal conductivity of the FG75 matrix than the sapphire fiber, axial implicit fluxes are generated in X direction, which tends the fiber to be compressed along this axis. As the implicit fluxes are usually small, they are multiplied by a constant factor of 10³ in the flux result files.

  

  Figure 18: Implicit fluxes in X direction generated at fibre/matrix interface due to jump of thermal conductivities (left) and the corresponding microscopic displacement component X (right).

Mechanical field results¶

If the parameter U or ALL is specified in the /OUTPUT keyword, each thermo-elastic mono- or multilayer homogenization analysis produces the output of an implicit forces file and of an microscopic displacement field file on the undeformed RVE.

• The 1^st file contains all the generated implicit forces at the domain interfaces, characterized by a jump in the material properties. If the temperature varies inside each domain, it contains also a volumetric contribution (see Appendix Homogenization of Materials with Periodic Microstructure). Its nomenclature is similar to the implicit flux file, except, of course, its extension differs: _hfor instead of _flux. In the present example, as two temperatures are specified ($T=800$°C and $T=1100$°C) on one layer, HOMAT writes following two implicit force files: fas30_fine_hfor11.vtk and fas30_fine_hfor12.vtk. The 1^st corresponds to the thermoelastic homogenization at 800°C and the 2^d one to the homogenization at 1100°C.
• The second file contains again the periodic microscopic displacement fields $\zeta^{rs}$ with $\text{rs}=1, 6$ in the 3-D case, on the undeformed RVE. This file is characterized by the extension _hoel. Otherwise, its nomenclature is identical to the implicit forces file. Therefore, following two periodic displacement files are written: fas30_fine_hoel11.vtk and fas30_fine_hoel12.vtk.

The number and nomenclature of different vector solutions in each of these files depends not only of the nature of the homogenization analysis (pure mechanical or thermoelastic) but also on the homogenization model type: 3-D or 2-D (plane strain or plane stress model):

• 3-D: in the case of a pure mechanical homogenization, we have six forces and microscopic displacements, corresponding each to a unitary macro-strain: either $E_{XX}$, $E_{YY}$, $E_{ZZ}$, $E_{XY}$, $E_{XZ}$ or $E_{YZ}$. Each force/displacement vector solution has 3 components in the VTK file. The implicit forces have the identifier IMPL_FORC_ followed by the extension of the nominal strain (e.g. XX). For example, HOMAT defines in the force file fas30_fine_hfor11.vtk, the solutions IMPL_FORC_XY and IMPL_FORC_YY, which means that an implicit force vector of 3 components is generated by the unitary strains $E_{XY}$ and $E_{YY}$ respectively. For the microscopic displacements on the undeformed RVE, the nomenclature is similar except, of course, the identifier name becomes DISPL_. Therefore, in the file fas30_fine_hoel11.vtk, following solution identifiers of microscopic displacements are then specified: DISPL_XY, DISPL_YY. In the case of thermoelastic homogenization, an additional implicit force solution due to the thermal expansion coefficient and the corresponding microscopic displacement solution are saved in the respective VTK files. They are characterized by specific identifier names: TH_EXP_FORCE for the implicit force and DISP_TH for the microscopic displacement field. Moreover, if volumetric eigenstrains are active, a second additional implicit force solution, with identifier TH_EIG_FORCE is written in the VTK file PROBLEM_NAME_hfor.vtk and the corresponding microscopic displacement field, with identifier DISP_EIG is written in the VTK file PROBLEM_NAME_hoel.vtk. In this special thermoelastic case, eight vector solutions are output.
• 2-D plane strain model: in the case of a purely mechanical homogenization, only four implicit forces, namely $F_{XX}$, $F_{YY}$, $F_{ZZ}$ and $F_{XY}$ are evaluated and the two planar components X and Y are saved in the implicit force file with the same identifier IMPL_FORC as in 3-D. In presence of a thermoelastic homogenization, one additional implicit force vector with two components, named TH_EXP_FORC and corresponding to the thermal expansion effect, is added. Moreover, if volumetric eigenstrains are active, a second additional implicit force vector solution, named TH_EIG_FORCE is introduced. The corresponding microscopic displacement file contents then either four (pure mechanical case), five (thermoelastic case, only thermal expansion coefficient) or six solutions (thermoelastic case with eigenstrain and th. expan. coeff.). The identifiers of these solutions are the same as for the 3-D case, namely: DISPL_XX, DISP_YY, DISP_ZZ, DISP_XY, DISPL_TH and DISPL_EIG.

• 2-D plane stress model: three implicit forces, namely $F_{XX}$, $F_{YY}$ and $F_{XY}$ are expressed in a pure mechanical homogenization analysis. Again only the planar components are saved per force or displacement, leading to three stored solutions with 2 components per solution in the VTK file. The nomenclature and the identifier names are the same as for the 3-D case. In a thermoelastic homogenization analysis, again additional in-plane solutions are expressed for the thermal expansion coefficient and, if specified, for the volumetric eigenstrain in each VTK file, leading to a total number of either four or five vector solutions. In Fig. 19 (left) the microscopic displacement in X direction, induced by the unitary strain $E_{XX}$ is drawn on the deformed RVE of the UD reinforced metal matrix composite and in Fig. 19 (right) the microscopic displacement in Y direction, induced by the unitary strain $E_{YY}$ is outlined. As the elasticity of the sapphire fibre is stronger than the elasticity of the FG75 NiAl matrix, the implicit forces in X direction, induced by the macro-strain $E_{XX}$, tends to extend the fibre in X direction (formation of an ellipse). Otherwise, the implicit force in Y direction under macro-strain $E_{YY}$ tends also to distort the fibre elliptically with the larger semi-axis in Y direction.

  

  Figure 19: Microscopic displacement in X direction, induced by unitary macro-strain $E_{XX}$ (left) and microscopic displacement in Y direction under macro-strain $E_{YY}$ (right).

Remarks:

• Due to the smallness of the thermal expansion coefficients, the corresponding microscopic displacements are also small and in order to better visualize these displacements they are scaled by a factor 10³. Note that the volumetric eigenstrains are usually larger, so that their corresponding displacements are not scaled by default.

• In order to represent a 2-D vector solution, the user must create an additional solution by using the calculator option in Paraview. For example, for the microscopic displacement solution DISPL_XX, the user introduces in the calculator following expression:

  iHat*DISPL_XX_X + jHat*DISPL_XX_Y + kHat*0
  

  This additional solution permits to use the “wrap by vector” option in order to visualize the desired planar vector solution, for example, DISP_XX.

• Since version 6.0, the user can also specify the parameters P and E in the /OUTPUT keyword. The parameter P, corresponding to the local pressure field, is set, by default, if stabilized mixed finite elements are defined in the thermo-elastic homogenization 3-D or 2-D analysis to model quasi-incompressible material behaviour. The additional output parameter E allows to display the periodic micro-strains $\vec{\varepsilon}(\vec{\varsigma}^{\text{rs}})$ induced by the variation of the microscopic displacements $\vec{\varsigma}^{\text{rs}}$ over the RVE. This way, the user can more easily analyse the evaluation of the effective Hooke matrix via expression \eqref{eq:eff_hoke_mat} and the effective eigenstress via expression \eqref{eq:eff_eigenstress}.

• In the presence of mixed finite elements, a scalar pressure field is saved for each resolved linear problem: The number of saved pressure fields in quasi-incompressible materials varies from 6 to 8 in 3-D homogenization, from 4 to 6 in 2-D plane strain or from 3 to 5 in 2-D plane stress homogenization. The file of periodic pressure fields is characterized by the extension: _hdru. Otherwise, its definition is identical to the nomenclature of the microscopic displacement file. For the example of the crystalline/amorphous lamella design at the nanoscale (see Fig. 9 in ²), the following pressure file is written: bubN_draw80_hdru11.vtk.

• As a purely mechanical homogenization is performed in 3-D in this example, this file contains the following scalar pressure solutions, named: PRES_XX, PRES_YY, PRES_ZZ, PRES_XY, PRES_XZ, and PRES_YZ. Moreover, as the parameter ALL was specified in the command file also the microscopic strain tensors are written for visualization purposes. The file of microscopic strains is characterized by the extension _hdefo and in the present example it is bubN_draw80_hdefo11.vtk. This file contains in 3-D the six micro-strain tensor compnents STRAIN_XX, STRAIN_YY, STRAIN_ZZ, STRAIN_XY, STRAIN_XZ and STRAIN_YZ. In the case of thermo-elastic homogenization, one or two additional strain tensors are written in the PROBLEM_NAME_HDEFO.vtk file. They are characterized by the identifiers STRAIN_TH and STRAIN_EIG respectively.

Note: By analogy to ABAQUS, the microscopic strain tensors contain engineering shear strains with $\gamma_{xy} =2 \varepsilon_{xy}$. Paraview easily interprets the six components of each microscopic strain tensor as tensorial components. However, for plane strain homogenization analysis, Paraview denotes the four strain components as 0, 1, 2, 3 which leads to following correspondence: $0 = XX$, $1 = YY$, $2 = ZZ$, $3 = XY$. If 2-D plane stress analysis is adopted, Paraview denotes the three stress/strain components X, Y and Z. In this case, we have simply following correspondence: $X = XX$, $Y = YY$ and $Z = XY$.

Stokes flow field results¶

By default, for a fluid flow homogenization analysis the parameters V and P are specified on the /OUTPUT keyword and HOMAT v6.0 prints for each mono- or multilayer homogenization three different VTK files:

• The first file contains all the generated implicit volumetric boost forces, corresponding to an unitary pressure difference applied in the structural direction X, Y, and Z respectively. The number of boost terms depends on the homogenization type: either 3 in the 3-D case or 2 in the 2-D case. The corresponding keyword in the VTK file is BOOST_B_*. The star * stands either for X, Y or Z in 3-D, specifying the structural direction of the boost forces. These boost forces are save as vector solutions in the corresponding VTK file. Its nomenclature is the following: to the problem name, FibBL_r2_bubLi, is added extension _boost11. To preserve the same definition as for the VTK definition of thermal and mechanical homogenizations, two numbers are also added: the first one specifies the considered layer and the second one, the considered homogenization. Note that they are always 1 in the case of Stokes flow homogenization.
• The second file contains the periodic microscopic velocities $\omega^{r}$, with $r=1, 3$ in the 3-D case or $r=1,2$ in the 2-D in fluid domains of the RVE. Of course, these velocities are zero in the solid domains of the RVE. Each velocity vector has three (or two) components. They are the solution of the special Stokes problems. HOMAT adds to the problem name the extension _hvel11 to specify the corresponding VTK file. The three or two velocity vector solutions in this VTK file have as identifier: VELO_B_*, where again * stands for X, Y or Z, respectively. Note that B in the identifier recalls that boost forces in the respective direction X, Y or Z are responsible for these velocities.

• The third VTK file contains the periodic pressure fields in the fluid domains of the heterogeneous material. A scalar pressure field is saved for each resolved Stokes flow problem. The number of saved pressure fields is either 3 in 3-D or 2 in 2-D fluid flow homogenization. The file of periodic pressure fields is characterized by the extension: _hdru11. In this VTK file, the different scalar pressure fields are denoted via identifier PRES_XX, PRESS_YY and PRESS_ZZ in 3-D. For the example of the flow transverse to a cylindrical fiber of radius 2 mm in a cubic box (10*10*10 mm³) the pressure fields PRES_XX and PRES_YY are drawn in Figs. 20 and 21, respectively.

  

  Figure 20: Pressure fields in a section of the cylindrical fiber due to boost forces in X direction.

  

  Figure 21: Pressure fields in a section of the cylindrical fiber due to boost forces in Y direction.

Localization Results¶

Micro-strains and stresses are evaluated on the RVE in selected, critical points of the thermomechanical structural analysis of the sample at the macro-scale. For each of these locations, total micro-strains are evaluated via expression \eqref{eq:loc_str} in the Gauss points of the elements and elastic micro-stresses are calculated via expression \eqref{eq:loc_str_1} in the same Gauss points. Independent of the finite element discretization, a linear extrapolation scheme and a dedicated averaging technique are used to express their corner node values. As strains are continuous fields through material interfaces, all the nodal contributions on both interface sides are taken into account in the weight function; whereas for discontinuous fields, like stresses, only contributions from the same material domain are considered.

For each localization run, the total strain tensor, an equivalent strain as well as the elastic stress tensor and the equivalent von Mises stress are saved in the element nodes of the RVE or of the selected layer. By specifying the monolayer option using the keyword /DEF LAYER, it is possible to select the layers for which micro-stresses and -strains are evaluated. This feature increases the efficiency of the program, as less CPU time is required to gain the desired microstructure results.

Note that, by analogy to ABAQUS, the total strain tensor contains engineering shear strains with $\gamma_{xy} =2 \epsilon_{xy}$.

The number of stress/strain components depends on the localization analysis type: six in 3-D plus the equivalent values. In 2-D plane strain analysis they are four plus the equivalent term and in 2-D plane stress analysis, three plus the equivalent stress/strain. In the VTK files appear either the identifier STRAIN or STRESS for both tensors and STRAIN_EQ or STRESS_EQ for both scalar equivalent measures. Note that in 3-D, Paraview interprets easily the six components of both tensors as tensorial components. But, for a 2-D plane strain homogenization analysis, Paraview denotes the four stress/strain components as 0, 1, 2, 3, which leads to following correspondence:

The micro-strain and -stress files have the following nomenclature: after the problem name, the extension _loeps specifies a micro-strain VTK file, which contains the components of the strain tensor as well as the equivalent strain; the extension _lostr specifies a micro-stress VTK file, which contains the component of the stress tensor, the equivalent stress and, for mixed elements, the pressure. To this extension is then added the number of the localization analysis. In the header of each VTK file, HOMAT writes the selected critical point of the macro-simulation (Gauss point of the considered macro-element, part of specified element set).

As an example, consider the localization analysis performed after the U-O forming simulation of a ferrite-pearlite pipeline tube (see Fig. 1). In the critical area ²² two Gauss points nearest the upper and lower surface of the X65 steel sheet are selected, where total micro-strains and -stresses are expressed in the ferrite-pearlite microstructure via the localization procedure. The VTK file names are respectively: rve2_loeps_1.vtk, rve2_loeps_2.vtk for the micro-strains and rve2_lostr_1.vtk, rve2_lostr_2.vtk for the micro-stresses.

In Fig. 22 (left) the most significant shear strain component is drawn on the RVE and in Fig. 22 (right) the equivalent stress distribution in the pearlite inclusions is illustrated. This strain distribution shows clearly that some elongated pearlite inclusions, oriented in rolling direction, are strongly sheared. As pure elastic localization analysis is performed and no plasticity is taken into account, the equivalent stresses in the pearlite inclusions are overestimated.

Log File Output¶

In this paragraph a brief overview of the most relevant output on the screen or - if redirected - in the log file, named either PROBLEM_NAME.hout, for a homogenization analysis, or PROBLEM_NAME.lout for a localization run, will be presented briefly. This output is so variable and depends on all possible combinations of input data, that a rigorous, detailed description of the log file output is impossible. Here some characteristic messages and information are highlighted. These messages will help the user to verify the correctness of its analysis. First, some common messages are discussed, then specific messages of a thermal homogenization run are presented, followed by messages of a thermo-elastic homogenization analysis and of the Stokes flow homogenization and, finally, by the description of the specific output of a localization run.

Common output¶

Each HOMAT run starts by the interpretation of the command file. The corresponding log file messages are self-explanatory. The first important check is performed in HLIGEO to verify if the bounding box around the RVE is correctly and that its dimensions in LENCHAR are correct expressed in the adopted unit length. The determination of nodes belonging to a periodic surface pair is a delicate operation, mainly in presence of different finite element (FE) meshes on both linked periodic surfaces. The user can check at first the correctness of the mean and minimum element wedge length used in the projection algorithm of the slave nodes on the master surface. Then, the output of routine PROPER for each periodic surface is important, as it gives the number of directly corresponding nodes and the number of nodes, which require, after projection on the master surface, an interpolation within the element containing the projection point. Globally, we can say that the more nodes need to be interpolated, the slower is the convergence of the iterative solver. Therefore, in order to reduce the number of nodes to be interpolated for incompatible surface meshes, the user can specify a larger value of the EPSCOORD parameter (see keyword /NUMERICS), for example 5.E-02 instead of 1.E-02. Sometimes, it can be interesting to check the number of detected double nodes per wedge of the RVE, provided by routine HLIECK. For the dimension of the homogenization or localization problem, the total number of Gauss points provided by routine FLAEPER is an important quantity. If only one layer is specified in the model (either multi- or monolayer homogenization), the volume and its fraction are expressed for each phase or material. These results provide important information as they measure, for example, the discrepancy between the discretized phases in the RVE/layer and the experimental or theoretical volume fraction of each phase, mainly if the MICRESS file .korn file is used to separate the individual grains of the microstructure. Independent of the considered homogenization analysis, the total surface area of each periodic surface is printed and help the user to verify the model's correctness. For each non-zero right-hand side vector, the reduced global linear system is solved iteratively. The output frequency of intermediate solver results, like the value of the convergence criterion and the actual norm of residuum, is governed by the OUTFREQ parameter of keyword /NUMERICS. The user can check via these outputs the convergence progress. The value of residuum and relative convergence criterion at the best iteration has to be analysed carefully. A too large relative convergence criterion ($\epsilon > 0.01$) is an indication of a badly configured global system. Bad systems occurs mainly in presence of large implicit forces or fluxes due to strong material discontinuity (e.g. crystalline/amorphous lamella) and sometimes in presence of 2-D periodicity inducing large microscopic fluctuation displacements on the RVE: Please, check if a 3-D periodicity B.C. cannot help to circumvent these too large fluctuation displacements on the RVE.

Thermal homogenization output¶

The norm of the implicit source terms is expressed at the end of the element generation in the routine GENVOTH. If a constant temperature pro domain is specified, no source terms are generated and their values have to be zero. Then, only implicit forces at the boundary between two different materials or phase are generated. These fluxes have to be in equilibrium in each structural direction. Please control their equilibrium. For each domain, the volume and the mean temperature are provided. The size of the reduced system, defined by the number of Degrees Of Freedom (DOF) is given by the routine RANDIR, please check this relevant information for the memory allocation. For each structural direction the maximum and minimum value of the implicit fluxes and their nodal location are outlined. If the normal at the material interfaces has no component in a specific direction, for example, in the third one (case of UD reinforced fiber composite), then the corresponding min. and max. values are zero and no contribution to the total implicit forces is expressed. Similarly, after the resolutions of the global reduced system, the min. and max. values of the predicted microscopic displacement field are written for each component. Eventually, the effective thermal conductivity tensor and the other effective thermophysical properties (thermal capacity, specific heat, density) are provided. The corresponding volume averaging values are also printed.

Mechanical homogenization output¶

The mechanical finite element generation starts with the output of the Hooke matrix in the structural axis system of the 1^st element of each domain. Please carefully verify these Hooke matrices. Again implicit volume forces are only evaluated in presence of temperature variations within a RVE domain or when the spherulite model is used. In this case, different lamella orientations in each Gauss point induce implicit volume forces. Then, for each identified interface with a jump of mechanical properties, the evaluated implicit mechanical forces are given for each component. In the case of thermo-elastic homogenization, also the total implicit forces due to unit thermal stresses, are printed. For each structural direction the maximum and minimum value of the implicit forces and their nodal location are given for the mechanical homogenization. Please check if the number of solved microscopic problems corresponds to the desired one for the defined homogenization analysis. Please check the number of DOF of the reduced system, printed by routine RANFIX; it directly impacts the memory allocation for the considered homogenization run. After the resolution of the different global systems, the min. and max. values of each predicted microscopic displacement component are given. Subsequently, the effective Hooke matrix before symmetrisation, the symmetric Hooke matrix and its inverse, the symmetric compliance matrix, from which the effective engineering moduli are derived, are provided. Thereafter, the Hooke matrix in the specified local axis system and/or the volume averaging engineering moduli are printed.

Stokes flow homogenization output¶

For each fluid domain, the total volumetric boost forces for each structural direction are provided by the GENFLOW routine. In addition, for each domain (solid or fluid) its volume is written. Further, for each material or phase, the volume fraction is printed, specifying the porosity of the porous medium. Please check this porosity for its correctness. After redistribution of the master nodal values, via parameter VISLAV=YES, the max. and min. boost forces are provided. Due to non-slip boundary conditions at fluid-solid interfaces, the reduced system formed of velocity and pressure DOF in the fluid domains is further reduced and the corresponding total number of DOF of the Stokes flow problem is printed. Then, after resolution of the three or two special Stokes problems, the min. and max. values of each predicted microscopic velocity solution and of the microscopic pressure field are written. Finally, the predicted Darcy permeability tensor in the structural axis system is written in the log file.

Localization output¶

The begin of the localization output is identical with the log output of the mechanical homogenization run, as the microscopic displacement fields have to be determined for the calculation of the microscopic strains and stresses respectively. After the resolution of the different microscopic problems, starts the real localization analysis for the different specified macro-strain states. For each localization analysis, the mean internal work per layer or for the whole multilayer is evaluated and printed. In the case, if the internal elastic work is expressed for the complete RVE, its value can be compared with the deformation work at the macro-scale and can thus be used to validate the Hill-Mandel principle ²³. Moreover, the mean internal work for each phase or material is given, showing so the contribution of each phase/material to the internal work of the layer or the RVE. Finally, the routine LOCWRITE performs a similar output of min. and max. nodal values, for each component, of the total micro-strain field and microstress field on the RVE or on the considered layer.

Transfer of macro-scale results from ABAQUS to HOMAT¶

In order to perform a linear elastic localization analysis with HOMAT, the user has to provide via the keyword /INPUT FILES, the name of the ASCII file of the linked thermo-mechanical macro-simulation of the part with ABAQUS. This ASCII file has always the suffix .txt. A dedicated python script, named Aba2Hom.py, has been written in order to extract from the data base of macro-simulation (.odb) field results (strains, stresses, temperature, von Mises stresses,...) in critical regions of the sample at selected load steps and iterations. Please look at the corresponding user guide of this python script to perform the desired extraction in critical regions of the specimen.

Here, the focus is given on the description of the transfer file NAME.txt, which is read by HOMAT. Example 36 shows such a file extracted from a plane strain thermo-elastic analysis of a 10° degree sector of a cooled cylindrical combustion chamber wall subjected to a stationary thermal gradient (see Fig. 23, left). The final deformation and the radial strain are plotted in Fig. 23 (right). The region of interest is magnified in Fig. 24. Strain, stress and temperature results are extracted by applying Aba2Hom.py with the analysis type -CritSelect. for the elements, labeled 10, 30, 50, 79 and 80. The generated transfer file is shown in Example 36.

Figure 24: Zoom of the central sector part. The 2-D finite element mesh of TBC and bondcoat layers is labelled there.

/HEADER
macro-strains extracted in the center of combustion chamber wall
Step-1
ODB Name:orth_lin.odb
ABAQUS Version 2019
HOMAT Version 6.0
Thu April 2 10:23:30 2020
/PARA, ABA_MOD=1, NEL=5, NPG=4
/ELSET=CRIT_TBC+BOND
/STRAIN
10,1,-0.00033,0.01487,0.00000,-0.00013
10,2,-0.00033,0.01487,0.00000,-0.00003
10,3,-0.00047,0.01486,0.00000,-0.00014
10,4,-0.00047,0.01486,0.00000,-0.00004
......................................
79,1,0.01610,0.01481,0.00000,0.00001
79,2,0.01617,0.01481,0.00000,0.00001
79,3,0.01610,0.01481,0.00000,0.00000
79,4,0.01617,0.01481,0.00000,0.00000
.....................................
/TEMP
10,1,931.26550293
10,2,931.26550293
10,3,907.40246582
10,4,907.40246582
..................

Aba2Hom.py six header lines. In the first two lines the user can specify a title describing the performed multiscale investigation. On the third line, the selected ABAQUS analysis step is given, whereas, in 4^th line the name of the extracted ABAQUS data base is written, here orth_lin.odb. Eventually, in line 5^th and 6^th the used program version of ABAQUS and HOMAT are provided.

The keyword /PARA specifies the used ABAQUS discretization model, via the mandatory model parameter ABA_MOD:

• ABA_MOD=0: 3-D macro-scale analysis with six strain components
• ABA_MOD=1: 2-D plane strain, axisymmetric or plane stress analysis with four strain components
• ABA_MOD=2: 3-D membrane or shell analysis with three strain components
• ABA_MOD=3: axisymmetric shell elements having two strain components
• ABA_MOD=4: truss elements having only one axial strain component

This keyword also provides the number of selected elements via the parameter NEL and the number of Gauß points per element via the parameter NPG. Then, the field results in the Gauß points of the selected elements are written in the transfer file. In Ex. 36, only the total strain tensor and the temperature field are transferred to HOMAT. But, note that in this .txt file other field results could have been written by the Python script aba2hom.py, such as the stress tensor, equivalent von Mises stress, etc. HOMAT neglects their presence and reads only the mandatory macro-strain results and in presence of thermal localization (parameter LOCATH on keyword /PARA of HOMAT command file) the temperature state. For each Gauß point, the script writes, the element label, then the Gauß point number followed by the field values in this point. As a 2-D plane strain analysis is defined in Example 36, four strain components are written: $E_{XX}$, $E_{YY}$, $E_{XY}$ and $E_{ZZ}$ in the structural axis system of ABAQUS, leading to following specification:

10, 3, -0.00047, 0.01486, 0.00000, -0.00014

Remark: If the structural axis systems differ at macro- and micro-scale, the user must specify the rotation matrix of the macro-strains in the micro-scale axis system. This mandatory rotation will be defined by a dedicated orientation definition via the optional parameter ORIENTATION on the /INPUT FILES keyword, placed there behind the MACRO parameter, specifying the macro to micro transfer file.

Homogenization of Materials with Periodic Microstructure¶

The asymptotic homogenization method will be briefly outlined here for a thermo-elastic material with a periodic microstructure subjected to volumetric eigenstrains. A more detailed and theoretical description of the method can be found in the literature ²⁴.

Static equilibrium of a heterogeneous material¶

The asymptotic homogenization method (AHM) assumes that the heterogeneous material is built at the micro-scale by the periodic repetition of a Representative Volume Element (RVE), At the macro-scale this heterogeneous material with volume $\Omega$ is in a static equilibrium with the body forces $b_{i}$, the prescribed surface tractions $\overline{t}_{i}$ on $\Gamma_{t}$ and the prescribed displacements $\overline{u}_{i}$ on $\Gamma_{u}$. The sample is further subjected to a temperature variation $\Delta T = T^{h} - T_{r}$, where $T_{r}$ is a reference temperature being set to 0°C without loss of generality.

Figure 25: Definition of a periodic unit cell for the example of an effusion cooled multilayer plate built of a CMSX-4 substrate, a MCrAlY bondcoat and Y₂O₃ stabilized ZrO₂ thermal barrier coating as used in gas turbines [^8].

Figure 26: Static equilibrium of a heterogeneous material.

In the macroscopic reference coordinate system $\vec{z} = (z_{1}, z_{2}, z_{3})$, see Fig. 25, the static equilibrium of the heterogeneous material is given by

As a pure elastic material behaviour of isotropic, orthotropic or anisotropic type is assumed for the heterogeneous material, the stress tensor $\sigma^{h}_{ij}$ is related to the total strain $\varepsilon^{h}_{kl}$ by following generalized Hooke law

The Hooke matrix $H^{h}_{ijkl}$ and the thermal expansion coefficient $\alpha^{h}_{kl}$ are both temperature dependent quantities. The volumetric eigenstrains $\kappa^{h}_{kl}$ depend only on the position Z and are generally present in polycrystalline materials. To complete the boundary value problem, the applied surface tractions $\overline{t}_{i}$ and displacements $\overline{u}_{i}$ are subject to the following conditions on $\Gamma_{t}$ and $\Gamma_{u}$, respectively

Periodicity and two-scale description¶

A heterogeneous material has a regular periodicity if its field variables $\varphi$ satisfy:

with

where $n_{1}$, $n_{2}$, $n_{3}$ are integers and $\vec{Y} = (Y_{1}, Y_{2}, Y_{3})$ is a constant vector, whose components describe the material periodicity (see Fig. 28).

In the homogenization theory, the periodicity $\vec{Y}$ is small compared to the specimen's dimensions. As the function $\varphi$, like e.g. the Hooke matrix, varies strongly in the neighbourhood of an arbitrary point $Z_{0}$, (see Fig. 27), two scales are introduced in order to describe the dependency of each physical quantity:

• a macroscopic, global variable $\vec{x}$, which measures the variation of the mean value of function $\varphi$ from unit cell to unit cell,
• a microscopic, periodic variable $\vec{y}$, which represents the strong variation of the function within each unit cell.

Figure 27: A strong oscillating function on the macro-scale (left) and a zoom of one period of this function on the micro-scale (right).

Figure 28: Periodic building of the heterogeneous material (a) and the corresponding RVE (b).

The scale ratio $\epsilon$ between the micro and macro-scale, defined by $y=x/\epsilon$, is usually a small factor, which specifies also the size ratio between the extensions of RVE and specimen. The material at the macro-scale can thus be described as the combination of identical unit cells of extension ($\epsilon Y_{1}$, $\epsilon Y_{2}$, $\epsilon Y_{3}$) (see Fig. 28).

The asymptotic homogenization method¶

Starting point of this method is a formal asymptotic expansion of the periodic displacement and temperature fields, $\vec{u}^{h}$ and $T^{h}$ on the heterogeneous material with the aspect ratio $\epsilon$:

As the displacement and the temperature field expansions are both limited to the 1^st order here, the corresponding homogenization scheme is named 1^st order homogenization method. Expansions up to the 2^nd order have been introduced by Koutznetsova et al. ²³ in order to capture microstructural size effects. This formalism is complex and leads to the introduction of an enriched macro-continuum of Cosserat type not being discussed further here.

The Hooke matrix $H^{h}_{ijkl}$, the thermal expansion coefficients $\alpha^{h}_{kl}$ and volumetric eigenstrains $\kappa^{h}_{kl}$ are assumed to be periodic function on the unit cell and their 1^st order asymptotic expansions are given by:

The series expansions Eqs. \eqref{eq:asy_expansion} and (\ref{eq:asy_expansion_props}) and the differential operator

imply the following stress definitions:

where $\vec{u}^{r}$ specifies either $\vec{u}^{0}$ or $\vec{u}^{1}$ and $E_{kl}$ and $e_{kl}$ are the macroscopic and microscopic total strains, expressed either in terms of $\vec{x}$ or $\vec{y}$ respectively:

These strains are introduced in the static equilibrium condition Eq. \eqref{eq:equil_condition}, leading thus to a system of differential equations, which can be sorted with respect to the scale ratio $\epsilon$. The differential equations related to the different powers of $\epsilon$ are discussed in the following:

• $\epsilon^{-2}$ term:

  $$\begin{equation} \frac{\partial}{\partial y_{j}} s^{0}_{ij}\left(\vec{x},\vec{y},T^{0}\right) = 0\,. \label{eq:eps-2} \end{equation}$$

  These differential equations are homogeneous. As the Hooke matrix is positive definite, $e_{ij}(\vec{u}^{0}) = 0$ and $u^{0}(\vec{x},\vec{y}) = u(\vec{x})$. Thus, the first term of the serial expansion, $\vec{u}^{0}$, depends only on the macroscopic variable $\vec{x}$ and not on the microscopic one. - $\epsilon^{-1}$ term:

  $$\begin{equation} \frac{\partial}{\partial y_{j}} S^{0}_{ij}\left(\vec{x},\vec{y},T^{0}\right) + \frac{\partial}{\partial y_{j}} s^{1}_{ij}\left(\vec{x},\vec{y},T^{0}\right) - \frac{\partial}{\partial y_{j}} p^{0}_{ij}\left(\vec{x},\vec{y},T^{0}\right) = 0\,. \label{eq:eps-1} \end{equation}$$

  This differential system is linear in $\vec{y}$. Its solution is composed of the global solution $\vec{u}^{1}(\vec{x},\vec{y})=\bar{\vec{u}}^{1}(\vec{x})$ of the homogeneous system plus a special solution of type

  $$\begin{equation} \label{eq:spec_solution} \vec{u}^{1}(\vec{x},\vec{y}) = \zeta^{rs}_{i} E^{0}_{rs}(\vec{x}) + \psi_{i}(\vec{y})\,, \end{equation}$$

  where $\vec{\zeta}^{rs} = ( \zeta^{rs}_{x}, \zeta^{rs}_{y}, \zeta^{rs}_{z})$, with $rs = xx, yy, zz, xy, xz, yz$ are periodic microscopic displacement fields induced by the initial strains $E^{0}_{rs}$. Furthermore, $\psi(\vec{y})$ are microscopic displacements generated by the eigenstress tensor $p^{0}_{ij}$ due to thermal expansion and/or volume changes of the material phases. The special solution \eqref{eq:spec_solution} is introduced into the system \eqref{eq:eps-1}. Then, the initial uniform strains $E^{0}_{rs}(\vec{x})$ are extracted from this system, leading to:

  $$\begin{multline} \label{eq:solu_1} E^{0}_{rs} \frac{\partial}{\partial y_{j}} \left[ H_{ijkl}\left(\vec{x},\vec{y},T^{0}\right) - H_{ijkl}\left(\vec{x},\vec{y},T^{0}\right)\left[ e_{kl}\left( \vec{\zeta}^{rs} \right) + e_{kl}\left( \vec{\psi} \right) \right] \right] + \\ + \frac{\partial}{\partial y_{j}} p^{0}_{ij}\left(\vec{x},\vec{y},T^{0}\right) = 0 \end{multline}$$

  As the initial strains and generalized eigenstrains can be specified on the unit cell separately, the modified system \eqref{eq:solu_1} can be split into:

  $$\begin{equation} \begin{aligned} & \frac{\partial}{\partial y_{j}} H_{ijkl}\left(\vec{x},\vec{y},T^{0}\right) = H_{ijkl}\left(\vec{x},\vec{y},T^{0}\right) e_{kl}\left( \vec{\zeta}^{rs} \right)\\ & \frac{\partial}{\partial y_{j}} p^{0}_{ij}\left(\vec{x},\vec{y},T^{0}\right) = H_{ijkl}\left(\vec{x},\vec{y},T^{0}\right) e_{kl}\left( \vec{\psi} \right) \end{aligned} \label{eq:solu_2} \end{equation}$$

  A variational formulation is applied to both systems \eqref{eq:solu_2} in order to determine the unknown displacement fields $\vec{\zeta}^{rs}$ and $\vec{\psi}$ ⁸. Then, the 3-D RVE of the heterogeneous material is discretized by isoparametric volume elements leading to following linearized systems in the nodal values $\vec{Z}^{rs}$ and $\Psi$ of the unknown displacement fields $\vec{\zeta}^{rs}$ and $\vec{\psi}$:

  $$\begin{equation} \begin{aligned} & \vec{K}\left( T^{0} \right) \vec{Z}^{rs} = \vec{g}^{rs}_{\text{im,H}}\left( T^{0} \right) + \vec{g}^{rs}_{\Delta T, \text{H}}\left( T^{0} \right)\\ & \vec{K}\left( T^{0} \right) \vec{\Psi} = \vec{g}^{rs}_{\text{im,p}}\left( T^{0} \right) + \vec{g}^{rs}_{\Delta T, \text{p}}\left( T^{0} \right) \end{aligned} \label{eq:solu_3} \end{equation}$$

  where:

  • $\vec{K}\left( T^{0} \right)$ is the classical elastic stiffness matrix of the RVE;

  • $\vec{g}^{rs}_{\text{im,H}}$ and $\vec{g}^{rs}_{\text{im,p}}$ are implicit surface tensions, which are induced either by the jump of the Hooke matrix or of the eigenstress at the phase/material interfaces;

  • $\vec{g}^{rs}_{\Delta T, \text{H}}$ and $\vec{g}^{rs}_{\Delta T, \text{p}}$ are implicit body forces, which are induced either by the temperature dependence of the Hooke matrix or of the eigenstress in each phase/material of the RVE.

• $\epsilon^{0}$ term: The solution of the $\epsilon_{0}$ term leads to the definition of following effective thermo-elastic properties of the heterogeneous material:

  1. effective Hooke matrix:

     $$\begin{equation} \label{eq:eff_hoke_mat} H^{\text{hom}}_{ijkl}\left( T^{0} \right) = \frac{1}{|Y|}\int_{Y}\left[ H_{ijkl}\left( y, T^{0} \right) - H_{ijrs}\left( y, T^{0} \right) e_{rs}\left( \vec{\zeta}^{kl} \right)\right]\,\text{d}y \end{equation}$$

  2. effective eigenstress:

     $$\begin{equation} \label{eq:eff_eigenstress} p^{\text{hom}}_{ij}\left( T^{0} \right) = \frac{1}{|Y|}\int_{Y}\left[ p_{ij}\left( y, T^{0} \right) - H_{ijrs}\left( y, T^{0} \right) e_{rs}\left( \vec{\psi}\,, \right)\right]\,\text{d}y \end{equation}$$

  where $|Y|$ is the volume of the RVE.

The generalized effective eigenstrains $\bar{\kappa}^{\text{hom}}_{ij}$ are obtained easily by inverting the relation:

Expressed for $\Delta T^{0} = 1$ Kelvin, relation \eqref{eq:gen_eff_eigenstrain} cannot simultaneously provide the effective thermal expansion coefficients $\alpha^{\text{hom}}_{ij}$ and the effective transformation strains $\kappa^{\text{hom}}_{ij}$. Therefore, two effective eigenstresses are evaluated successively: one including the volumetric strain $\kappa_{ij}$ and one considering only the thermal expansion coefficient. In addition to the purely mechanical problem (determination of the 6 periodic displacement fields $\vec{\chi}^{\text{rs}}$), the new 3-D thermo-elastic procedure will consider either 2 additional effective eigenstresses (in the presence of eigenstrain and thermal expansion) or one additional eigenstress field (only thermal expansion). This procedure is implemented in HOMAT since version 5.2 and for 3-D thermo-elastic homogenization the solution of either 8 or 7 linear systems.

Finally, in order to get the effective engineering properties, the predicted Hooke matrix is made symmetric and then inverted. From the effective compliance matrix the orthotropic Young ($E_{x}$, $E_{y}$, $E_{z}$) and shear moduli ($G_{xy}$, $G_{xz}$, $G_{yz}$) and the Poisson coefficients ($\nu_{xy}$, $\nu_{xz}$, $\nu_{yz}$) can be derived. In the special case of a quasi cubic global behavior, not only the effective Young and shear moduli and the Poisson coefficient $\nu$ can be extracted, but also an effective Zener anisotropy coefficient, defined by

Post-processing of macro-scale results: the localization step¶

The coupling to the macro-scale simulation is realized by a localization step comprising the evaluation of the micro-strains and -stresses on the RVE. The thermo-mechanical analysis using effective properties on the macro-scale allows the evaluation of macro-strains and -stresses in the material. In regions with high stress concentrations it is also important to accurately know the micro-stress and -strain levels on the RVE. They are needed to predict the damage and/or failure initiation criterion on the RVE and to localize critical zones of the microstructure (see Figs. 29 and 30).

Figure 29: Composite plate under uniaxial tensile loading

Therefore, at selected load increments the macro-strain state is extracted and transferred to the RVE. It specifies there the boundary conditions. In the framework of the asymptotic homogenization method (see Appendix Homogenization of Materials with Periodic Microstructure), the asymptotic development of the displacement field \eqref{eq:asy_expansion} is introduced in the total strain expression

The special solution \eqref{eq:spec_solution} for the displacement term $\vec{u}^{1}(\vec{x},\vec{y})$ allows to transform expression \eqref{eq:loc_str}, as follows

This expression specifies the total micro-strains in the RVE for the given macro-strains $E_{ij}^{0}(\vec{x})$. Similarly, by taking into account the definitions [Eq. \eqref{eq:stresses}] of the macro- and micro-stresses, $S_{ij}$ and $s_{ij}$ and the eigenstress $p_{ij}^{0}\left( \vec{x},\vec{y},T \right)$, the microscopic stresses are given by:

Based on these local stress/strain results on the RVE, a material dependent failure criterion can be applied in order to localize failure initiation in the heterogeneous material's microstructure at a specific load level.

Periodic boundary conditions via a master-slave projection algorithm¶

The application of the periodic boundary conditions on the RVE via the penalty method works but leads to a poorly conditioned system matrix with high differences between the stiffness terms and induces locally additional artificial stiffness. In order to more accurately predict effective mechanical properties, a master-slave projection algorithm, detailed in Zhang and Oskay ²⁵, has been implemented in HOMAT version 6.0. Complex RVE's, like the cross-hatched PP bi-lamella mesh are often at best close to but not strictly periodic. In strictly periodic microstructures, the periodicity is simply enforced by assigning master-slave pairs on opposite RVE surfaces and eliminate the slave nodes from the global systems \eqref{eq:solu_3}, as they have the same displacements at the corresponding master nodes. For nonmatching meshes on opposite periodic surfaces, we specify the surface with the finer mesh density to be the slave surface and the other one to be the master one. Then, each node, belonging to the slave surface element is projected on the master surface (see Fig. 31). Its specific projection point in a master-element is identified, allowing its interpolation in this specific element. Performing this nodal projection for all nodes on a slave surface, yields the following relationship between slave and master displacement vectors

where $\overline{Z}_{s}^{rs,e}$ and $\overline{Z}_{m}^{rs,e}$ denote the vectors of slave and master nodes respectively and $\overline{C}$ is the rectangular projection matrix build by the local interpolation coefficients. Expression \eqref{eq:mast_slave} is then used to eliminate the slave nodes per surface element of the slave surface. The global matrix systems \eqref{eq:solu_3} are thus reduced and contain only bulk and surface master nodes.