logo simtec
company software case studies literature events main
Numerical methods of the SIMTEC simulation software

1. Non stationary temperature field calculation (TFB)
    1.1 Basic proceeding
    1.2 Internal Heat Transfer
        1.2.1 Validation
    1.3 Latent Heat
    1.4 Solidification Morphology - Undercooling - Microstructure

2. The Calculations of Thermal Stress, Outer Load and Distortion (SPA3)
    2.1 Basic Proceeding
    2.2 The Plastic Approach
        2.2.1 Validation

3. The Calculation of the Form Filling (FILL)
    3.1. Fundamental Proceeding

1. Non stationary temperature field calculation (TFB)

1.1 Basic proceeding

The basis for the calculation of the heat conduction problem is the heat conduction equation of Fourier:

In this partial differential equation, the temperature T is described as a function of the space coordinates X,Y,Z and the time t. The temporal change of the temperature at a given location results of the flow of heat quantities, plus a term for the heat origin and depression. For the unambiguous description of a temperature field, there are required additional boundary conditions and starting conditions. The starting condition

yields a description of the temperature field to the time t=0. The boundary conditions supply information on the temperature on the surface of the object at a time t>0. Usually there are two different functions used for that. The Dirichlet - condition

contains the statement of an actual temperature as a function of location and time, whereas the Cauchy - condition

includes the spatial derivation of the temperature orthogonal to the surface into the description.

The solution of the differential equation is achieved by solving an equivalent variation problem

In order to solve this functional, the geometry is subdivided into small (finite) elements, in our case into prisms. All prisms are mapped by a transformation onto an unity prisms. Instead of integrating over the whole geometry, we only integrate over the standardized prism. Here we choose a linear approach for the calculation of the temperature distribution in the interior of the prism. This means we assume that the temperature changes in the interior of the element along a straight line. The more the actual temperature field differs to this linear approach, the greater the inaccuracy in the numerical calculation will be. But the integration for those prisms can only be done with such a particular assumption. The minimization will be solved, by demanding that the first partial differentiation vanishes. Thus we yield a system of ordinary differential equations of the for

The solution to this system of equations is achieved by applying backwards the Euler - procedure [2],

which provides reliable results also on bigger time steps. In this way we are able to reduce the partial differential equation (1.1) to a system of linear equations (1.7), which is solved by the Gauŝ - Seidel - procedure [3]. When stating and solving this system of equations, we also have regard for minimizing the used disk space as well as the number of numerical operations as possible. Here we can make use of the fact that the matrix

from equation (1.7) is symmetric, positive definite and contains much zeros. There are possibly 1,000,000 x 1,000,000 elements in the matrix (according to the number of nodal points), but in one line there are only as many elements that are non-vanishing, as the corresponding nodal point has direct neighbors, that means adjacent nodes. By a skillful organization of the matrix we are able to avoid the storing of many of the zeros and the corresponding operations.

1.2 Internal heat transfer

Whereas Stafford et al.
[9] 1987 suggested the usage of contact elements of the thickness 0 to calculate the heat transfer between different materials, we consider another method. In general to consider the heat transfer it is assumed that the heat quantity that is transferred at a phase boundary is proportional to the difference of the temperatures

The factor of proportionality a (alpha) considers all the material- and technique depending parameters, that influence the heat flow. Comparing equation (1.8) with (1.4), the according to g (gamma) is the heat flow Q. The term

that is the dependency of the heat flow to the temperature gradient orthogonal to the surface, is neglected in (1.8). The introduction of an internal boundary at the place of contact of casting to mold enables us to treat equation (1.8) as a Cauchy - condition [10] with

1.2.1 Validation


1.3 Latent heat

At the solidification of metals there is a internal energy that is freed in correlation with the assembling of the crystal grid, which affects the behavior of the temperature of the casting in a decisive way
[4]. The temperature correction method regards this energy. To this purpose the quantity of energy that is set free has to be known as a function of temperature. You can calculate the heat quantity that is added or subtracted at an element from its temperature fields to the time t and t+(t. This heat quantity is then divided into one part for the change of temperature and another part for the change of latent heat. Afterwards the temperatures are corrected according to this division. This method guarantees the preservation of the total of energy and can possibly be improved even more for bigger quantities of heat by an a-priori estimation.

1.4 Solidification morphology: under-cooling - microstructure

Depending on the alloy composition the part of the freed quantity of latent heat in the liquidus - solidus - interval may differ to a degree from the assumed linear model. Such a distribution function can be implemented as a material dependent parameter and be regarded. If the distribution function is not known, a standardized one will be chosen. Another factor that influences the solidification is the Keim - state of an alloy. Depending on the gradient with which the material enters the liquidus - solidus - interval, there may occur a remarkable under-cooling during the pre-solidification, as well as during the eutectic solidification. This also can be calculated.

For the calculation of the microstructure it is not enough to regard only the temperature behavior up to reaching solidus, but there perhaps may be a change of the crystal grid even at considerably lower temperatures. For the calculation of microstructure we analogize the cooling curves with the
time-temperature-change diagrams.

The intersections of the microstructure lines and the cooling curves yield a specific microstructure for each nodal point. The energy freed during the change of the crystal grid can be calculated using the temperature correction method, as explained above for the latent heat.

2. The calculation of thermal stress, outer load and distortion (SPA3)

2.1 Basic proceeding

The method of Finite Element is based on the fact that in each point of an object that is in no motion, the sum of all forces has to be zero. This concludes a system of equations that provides for each point three equations, that are the forces in the three coordinate directions, the parameters are nodal translation. First the distortion of the elements coming from the nodal translations are calculated. This is done under the assumption, that the tension in one element is constant. Then one is able to calculate the tension in one element with the translations of the nodal points. (For simplicity we developed the following equations for two dimensions.) So we yield the following equations

These tensions are the summation of the elastic and thermal tensions. We calculate the stress with the elastic tensions following Hookís theorem [6]

The forces at the nodal points are determined following the principle of virtual labor out of the stress. The labor, done by all the exterior forces on the object, has to be at an arbitrary (virtual) translation the equal to the labor, done by the internal forces:

As all values in one element are considered constant, the integral can be calculated easily

In this way it is possible to calculate the forces of one element in its nodal points out of the nodal translations.

2.2 The plastic approach

Hook's Theorem indicates the linear correlation between the stresses and the tensions of one element. This correlation however may not necessarily be valid always, foremost at the cooling down of a casting, since the interval described by Hookís Theorem is becoming very small at higher temperatures. Usually one tends to apply the fluid theory to describe the elastic-plastic material properties. We characterize the stress/tension correlation by three further material properties
  1. a flow condition, determining the multi-axial stress condition, at which starts the plastically flow
  2. a flow rule, linking the plastically distortion increments with the actual stresses and the stress increments due to the flow and
  3. a solidification rule, that states in which way the flow condition is modified during the plastically flow.
As the plastically material behavior is introduced into a system of finite element equations, we yield a non linear system of equations instead of equation (2.4). During the constitution of the system of equations it is not distinct whether an element shows elastic or plastic material behavior. This decision has to be made afterwards by examination of the flow condition. Thus it is not yet known, when building the stiffness matrix (elasticity matrix), which material rule to follow. If it turns out for one element to show plastically material behavior, there wonít be additional linearities, as the stress is included as a parameter in the stress/distortion matrix. To solve such a system of equations one usually employs the incremental method. This method assumes, that the solution is known at a time t. Now you achieve the solution at a time t+(t with the superposition

The tangential stiffness matrix is calculated out of the geometrical and material dependent premises at a time t. Dujis corresponding to the incremental increase of the nodal translations. We yield a system of equations with which we can calculate those incremental nodal translations. Thus we are able to calculate the translation to the time t = tD.

This on the other hand results that we can calculate the stresses at the time t = tD and we can proceed with the next time increment. For the calculation of the tangential stiffness matrix we make use of the stress/tension relations, that differ to a high degree, dependent on the exterior load. In the elastic case there are only as material parameters the E-module (module of elasticity) and the density reduction factor. If however there is a plastic flow, we have to use the stress/tension matrix. In praxis the research for the necessary data is very difficult, because the stress diviator as well as the distortion/solidification module are not known beforehand, especially in the non-isothermal case. For example Bathe[7] reduces the applicability of the stress/tension matrix expressis verbis to isothermal conditions. For the calculation of the thermal stress, the SIMTEC software package makes use of a modified model of the tangential stiffness method, bypassing this difficulty. If the load respectively time increments are sufficiently small, the stress / distortion relation can be linearised gradually by employing the tangential module rather than the E - module in the increment. For every time increment t = tD we can calculate the tangential module using the stress comparison of van Mises and the temperature T from the data coming from a strain experiment. But also on this method it is important to make an inaccuracy estimation. Especially the reduction of load on an element evolves parallel (as is generally accepted) to the original module of the stress - tension - diagram. This means for this case, that instead of using the tangential module we can take the original module for calculating the stress and there may be more iterations necessary. In fact the practical calculation of different cooling downs could show that usually this inversion of stress is done once or twice per element, and usually almost simultaneous in each element of an object.

2.2.1 Validation

To proof this plastic approach we made use of the semianalytical solution of H.G. Landau
[11] for tempering an infinite plate. According to this, one is able to calculate the process of the residual stress of a plate with ideal elastic-plastic material behavior, which has a known temperature distribution. For the present uni-dimensional heat conduction problem, the temperature distribution was determined numerically with a differentiation method and from this the process of residual stresses was calculated analytically. See also Honsel - Graduation [6] p. 74

3. The Calculation of the Form Filling (FILL)

3.1. Fundamental proceeding

In accordance with the &#quote;classical&#quote; balance consideration for the nodal point of an finite element, it is

In this case the sum of interior forces considers the following physical single forces

The individual forces can be computed as follows:

The calculation of the forces of friction is based on the stresses that effect the fluid - elements.

This equation corresponds formally to the tri-dimensional Hook's Theorem, from which we determine the forces at the nodal points in the linear - elastic stress calculation. For the calculation of the form filling, the forces caused by friction can be calculated following the principle of virtual labor for linear changing velocities within an element. Leaving aside the forces of inertia at the first step, this approach leads to a system of linear equations, that for each nodal point calculates the three unknown components of velocity from the three balance conditions

Starting point of a form filling calculation is for example a completely filled sprue. In order to keep the expenditure for the calculation low, the balance consideration is only carried out at the fluid front. This leads to a considerable decrease in the expenditure, whereas a detailed observation of the filling behavior however, is further given. At the filling front, we receive from the FEM calculation a vector for the velocity at each nodal point at the filling front, with which we can calculate the new front for the next time increment.


Richter, W. :
Numerische Lösung partieller Differentialgleichungen mit der Finite-Elemente-Methode

Vieweg Braunschweig (1986)

Grigorieff, R. D.:
Numerik gewöhnlicher Differentialgleichungen Band 1

Teubner (1972)

Schwartz, H. R.:
Numerische Mathematik

Stuttgart (1986)

Weiŝ, K.:
Temperaturfeldberechnung bei Erstarrungsvorgängen unter Ber¸cksichtigung des Einf¸llvorganges
Dissertation Aachen (1986)

Hahn, H. G.:
Einf¸hrung in die Methode der finiten Elemente in der Festigkeitslehre
Frankfurt (1975)

Honsel, Ch.:
Die Berechnung von W”rme- und Eigenspannungen infolge von Abk¸hlprozessen mit der Methode der tangentialen Steifigkeiten
Dissertation Aachen (1992)

Bathe, K. J.:
Berlin (1986)

Truckenbrodt E.:
Fluidmechanik Band 1
Springer-Verlag Berlin 1980

Stafford, R. O.; Rice, A. B.; Pinella D. F.:
Investment Casting Process Design - Part 2: Solidification Simulation

Proc. ASME Intl. Computers in Engng. Conf. (1987), 435-445

Prinz B.; Neite G.; Schulze T.; Klodt, J.; Honsel C.
Messung und Modellierung des Einflusses von Schlichten auf den W”rm¸ebergang zwischen Schmelze/Guŝst¸ck und Kokille beim Dauerformguŝ von Aluminiumlegierungen
Abschluŝbericht zum BMFT Vorhaben 03 K 4301 (COST504) 1995

Landau H.G.; Weiner J.H.
Transient and Residual Stresses in Heat-Treated Plates

J. of Appl. Mech. 25 (1958), 459 - 465

Landau H.G.; Weiner J.H.; Zwicky E.E.
Thermal Stress in a Viscoelastic - Plastic Plate with Temperature Dependent Yield Stress

J. of Appl. Mech. 27 (1958), 459 - 465