Stress Simulation

The Hooke is a PoligonSoft solver designed for calculating the stress-strain state (SSS) of a cast and a mold. The basis for calculating the SSS is the files of the results of the previously performed solidification calculation. The solidification and porosity solver Fourier creates .cst and .mld files that contain the temperature fields of the casting and mold, respectively. These files are part of the input data for the Hooke stress solver.

As a result of calculating SSS, the Hooke solver generates files *_stress.nds and *_stress.u3d. The main *_stress.nds results file contains several fields:

• Effective Strain;
• Effective Stress;
• Displacements;
• Cracks.

Files with *_stress.u3d extension are optional, there can be several, each file contains only one scalar field. Files of both types can be visualized and viewed in the Mirage postprocessor, which allows you to view the results not only after the end of the calculation, but also directly during the work of the Hooke solver, automatically determining the current last calculated step.

Stress Simulation Parameters

The parameters required to run the solver are located on the Stress tab (see the figure below).

Stress Simulation Parameters

• Enable Stress Model. This parameter enables Stress simulation. If the parameter is disabled, the calculation of stresses is not performed and the parameters of the Stress tab are ignored. The parameter is disabled by default.



• Calculation Timestep, s. This parameter indicates how often the stresses simulations should be performed. It is not related to the temperature simulation step. This means that if the CST-file records the results every 10 seconds, the stress simulation can be carried out with an independent step of 1, 7, 10, 15, etc. seconds. If the step with the thermal simulation does not coincide, the temperatures at the nodes are the linear interpolation between two steps closets by the time. If automatic step selection is enabled parameter is used as the base value (see the description of the UseVariableTimeStep parameter in the Additional Parameters of the Hooke Solver - Automatic Timestep Tab section).
  Enter a float value in seconds.
  Default value: 50.0.



• Contact Interface with Mold. The parameter specifies the mode of the contact interaction of the casting with the mold:


• Not simulated. In this case, the form does not participate in the calculation and does not in any way affect the formation of the SSS of the casting. With this value of the parameter, it is required to set additional boundary conditions for fixing the casting in space. This is done using the Boundary Conditions File or Automatic Model Fixation parameters (see below).


• Rigid Mold. In this case, the mold affects the formation of the SSS of the casting, because resists changing its (casting) dimensions. In this case, the mold itself is not deformed and its SSS is not calculated. This is the most commonly used type of contact between casting and mold, because it does not require setting additional boundary conditions for fixing the casting in space and significantly increases the reliability of calculations. With this type of contact, the Boundary Conditions File and Automatic Model Fixation parameters are not required.


• Linear Elastic Mold. More complex case of interaction of the casting with the mold than the previous one. The mold not only influences the stress-strain state of the casting, but itself deforms and changes dimensions under the influence of temperature. Since the calculation of SSS is carried out not only in the casting, but also in the form, its duration can significantly increase. In addition, it is required to set additional boundary conditions for fixing the shape in space. This is done using the Boundary Conditions File or Auto Fix Model parameters (see below).


• Elasto-plastic Mold. The most complex and accurate model of form behavior. An Elasto-Plastic Mold model requires elastic properties of the material, as well as the Yield Strength and Hardening. In addition, boundary conditions are required to fix the mold in space. This is done using the Boundary Conditions File or Auto Fix Model parameters (see below).

If contact algorythm is enabled, the input cst file must contain a finite element mold mesh. And in the case of a contact with a linear elastic or elasto-plastic mold, the mld file with temperature fields is also required to run the solver. In the absence of this data, the calculation will be completed with an error.
Default value: Not simulated.


• Mold Stress Properties File. The parameter specifies the full stress properties filename of the mold material .bdd.The content of the file can be changed in the editor, which is started from the parameter line (for more details see the Introduction section). The user can change the mechanical properties of the mold and save them (including with a new name). After closing the editor, the file information in the parameter line will be updated.



• Boundary Conditions File. The field specifies the path to the mcst file containing the information about casting fixed nodes. This file is created in the postprocessor Mirage (for more details see the Nodes Fixation section). The parameter is ignored when the simulation of the contact interaction of the casting with the rigid mold is included as well as when using the automatic fixation of the model in space.



• Mold Boundary Conditions File. The field specifies the path to the mcst file containing the information about mold fixed nodes. This file is created in the postprocessor Mirage (for more details see the Nodes Fixation section). The parameter is ignored if: there is no Mold in model, the simulation of the contact interaction of the casting with the rigid mold is included as well as when using the automatic fixation of the model in space.



• Automatic Model Fixation. This parameter controls the automatic assignment of boundary conditions for fixing the nodes of the model, excluding the possibility of moving it as an absolutely rigid body (for more details see the Nodes Fixation section). This is an alternative to the Boundary Condition File parameter that requires a decision from the user. The parameter is disabled by default.



• Output Timestep, s. This parameter specifies how often the stress simulation results should be saved to the nds file. The stress simulation files contain much more information than the other results files of the PoligonSoft. Therefore, their size can be quite large. To be able to manage the size of the stress simulation file, you can specify a saving step different from the simulation step (but not less than itself). If the saving step is set lower than the simulation step, the simulation will be performed at the frequency set by the saving step.
  Enter an integer value in seconds from 1 to 10000.
  Default value: 100.



• Threshhold Tension, MPa. The parameter controlling the start of saving the results to disk. If the maximum stress in the casting is less than the value specified by the parameter, saving to the disc is not performed. Using this parameter allows you to "postpone" the start of saving the results to disk.
  Enter a foat value in MPa.
  Default value: 0.



• Simulate Stress after Mold and Gating System Removal. The parameter includes the simulation of the stress on the additional technological stages: after demolding (knockout, disassembly of the chill mold, etc.) and after cutting the elements of the gating and feeding system (see Stress Simulation at subsequent technological stages section for more details). The parameter is disabled by default.



• Mold Removal Temperature, ˚C. The parameter specifies the temperature value at which the casting is taken from the mold (see Stress Simulation at subsequent technological stages section for more details).
  Enter an integer value in Celsius degrees. A negative value of the parameter means that the temperature field of the casting is taken from the last step of the cst temperature file.
  Default value: 300 ˚C.



• Gating System Removal Temperature, ˚C. The parameter specifies the temperature value at which the gating and feeding system is cut (see Stress Simulation at subsequent technological stages section for more details).
  Enter an integer value in Celsius degrees. A negative value of the parameter means that the temperature field of the casting is taken from the last step of the cst temperature file.
  Default value: 20 ˚C.



• Gating System Removal Order. The parameter that determines the removal order of the elements of the gating and feeding system. A text string consisting of numbers (indices) of the deleted volumes separated by the spaces (for example: 2 6 4) is entered in the field. The order of the specified numbers specifies the order of removal of the corresponding volumes from the stress simulation. Volume index 1 is ignored because this is the casting index.Volume index 1 is ignored because this is the casting index. An empty string mean the refusal to simulate the stress when removing gating system (see Stress Simulation at subsequent technological stages section for more details).
  Enter:
  • text string, for example: 2 6 4 - for the removal order of the gating system volumes;
  • an empty string – refusal to simulate the stage.
  Default value: empty string.



• Write x-conponent of Displacement (in mm). The parameter allows saving the x-component of the nodes displacement vector to .U3D universal file in millimeters. The file gets the name [NDS_NAME]_Ux_mm.u3d, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The parameter is disabled by default.



• Write y-conponent of Displacement (in mm). The parameter allows saving the x-component of the nodes displacement vector to .U3D universal file in millimeters. The file gets the name [NDS_NAME]_Uy_mm.u3d, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The parameter is disabled by default.



• Write z-conponent of Displacement (in mm). The parameter allows saving the x-component of the nodes displacement vector to .U3D universal file in millimeters. The file gets the name [NDS_NAME]_Uz_mm.u3d, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The parameter is disabled by default.



• Write Effective Plastic Strain. The parameter allows saving the Effective Plastic Strain to .U3D universal file. The file gets the name [NDS_NAME]_Eip.u3d, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The parameter is enabled by default.



• Write Principal Stress 1. The parameter allows saving the Principal Stress 1 to .U3D universal file. The file gets the name [NDS_NAME]_S1.u3d, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The parameter is disabled by default.



• Write Principal Stress 2. The parameter allows saving the Principal Stress 2 to .U3D universal file. The file gets the name [NDS_NAME]_S2.u3d, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The parameter is disabled by default.



• Write Principal Stress 3. The parameter allows saving the Principal Stress 3 to .U3D universal file. The file gets the name [NDS_NAME]_S3.u3d, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The parameter is disabled by default.



• Write Gap beetwen Casting and Mold. The parameter allows saving the Gap value beetwen casting and mold to .U3D universal file.in millimeters. The file gets the name [NDS_NAME]_gap_mm.u3d, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The parameter is disabled by default.



• Write Hot Tearing Indicator (for Mold Casting). The parameter allows saving the Hot Tearing Indicator to .U3D universal file. The file gets the name [NDS_NAME]_HTI.u3d, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The parameter is disabled by default.

Additional Parameters of the Hooke Solver

The Hooke solver has many parameters that allow the user to control flexibly the work of the solver. These parameters are not simulation parameters, so they are not on the Stress Tab. To access the solver settings, click the Advanced button on the Stress Tab (see the figure below).

Advanced Hooke Solver Settings button

The button opens the preprocessor settings window in the Stress Model section (see the figure below). There are many solver parameters that are grouped in a different tabs.

Hooke Solver Settings

In general, you may not change these parameters and use the default values, which in most cases will be optimal. However, getting the experience with the Hooke module, the user can look at the processor settings to have an access to the additional (advanced) capabilities of the solver.

General Tab

• Required Accuracy of Stress-strain State Calculatiob. Parameter sets the required accuracy of the simulation of the stress (the maximum modulus of the components of the generalized forces brought to the mesh nodes, the constituents corresponding to the fixed components of the corresponding nodal displacements are not taken into account in the boundary nodes).
  Default Value: 10-7 (recommended).



• Limit of Newton's Method Iterations. Parameter sets the maximum allowed number of the iterations of the Newton method. Consideration of the contact interaction in the stress simulations usually requires a significantly larger number of the iterations for the convergence of the Newton's method as compared with the simulation without taking into account the contact interaction. With the simulation of the contact interaction enabled, the maximum number of the iterations should be at least 1000.
  Default Value: 1000 (recommended).



• Criterion Used for Determination of "Damaged" Nodes. The parameter sets the conditions under which the destruction is recorded at the mesh node of the casting. Valid values:


• Turn off Damage Modelling.


• Effective Stress Criterion. A criterion σ_u ≥ σ_B is used, that means that the crack is fixed at the mesh node if the effective stress at the node exceeds the tensil strength at the current temperature at the node. Tensile strength as a function of temperature should be set in the alloy stress properties file.


• Maximum Principal Stress Criterion. A criterion σ_I ≥σ_B is used, that means that the crack is fixed at the mesh node if the maximum principal stress at the node exceeds the tensil strength at the current temperature at the node. Tensile strength as a function of temperature should be set in the alloy stress properties file.


Default Value: Turn off Damage Modelling.


• Upper Bound of Poisson's Ratio. Parameter determines the maximum permissible value of the Poisson's ratio. The Poisson's ratio cannot take the values of 0.5 and above in the file of mechanical properties of the alloy * .bdd. If this condition is violated during the simulation, the value of the Poisson's ratio automatically becomes equal to the PoissonMax parameter value.
  Default Value: 0.49.

Calculation stoppage Tab

• Termination of Computation when the Specified Percentage of "Damaged" Nodes is Reached. Parameter sets the fraction of "damaged" mesh nodes of the casting, when the simulation stops abnormally.
  Enter an integer value in % from 0 to 100.
  Default Value: 200 (i.e. the parameter is ignored)



• Termination of Computation when Specified Timeis Reached. Parameter sets the time point (specified in seconds) upon which the simulation is interrupted (use it if necessary to limit the time interval in which the simulation is performed). If the value is negative, it is ignored.
  Default Value: -1.



• Maximum Number of Computation Steps. Parameter sets the maximum allowed number of simulation steps (if set to negative, the parameter is ignored).
  Default Value: -1.

Diagnostic file creation Tab

• Write *_StepStat.txt. The parameter enables writing file [NDS_NAME]_StepStat.txt , where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The option is disabled by default.



• Write *_NewtonStat.txt. The parameter enables writing file [NDS_NAME]_NewtonStat.txt, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file). The option is disabled by default.



• Write *_IDataStat.txt. The parameter enables writing file [NDS_NAME]_IDataStat.txt, where NDS_NAME is the name of the stress simulation file (the file is created in the directory containing the stress simulation file).The file IDataStat.txt contains he diagnostic information about step selection by the time based on the analysis of absolute and relative changes of the input task data with the transition from the current (already “simulated”) time moment to the next one. The option is disabled by default.

Automatic Timestep Tab

• Use Variable Time Step. The parameter enables simulation with a variable time step, selected automatically based on the analysis of the rate of change in time of the input data of the task. The option is enabled by default.



• Minimum Value of Time Step, s. The parameter sets the smallest permissible value of the step in time, in seconds.
  Default value: 1.



• Contraction Multiplier for Time Step. The parameter sets the value, which (in case of dissatisfaction with the conditions of limitation of changes in the input data) is multiplied by the trial value of the time step.
  Enter a float value from 0 to 1.
  Default value: 0.5 (is recommended).

The following there are the parameters that determine the maximum allowable values of the absolute and relative changes in the input data of the task. The parameters, which names end with "atol" and "rtol", correspond to the absolute and the relative changes, respectively. In the event that the value of any of these parameters is set to negative, then the corresponding check will not be performed when selecting the time step value.

• Constraint for Absolute Time Variation of Temperature. Limitation for the absolute change of the temperature value.
  Default value: -1



• Constraint for Relative Time Variation of Temperature. Limitation for the relative change of the temperature value.
  Default value: 0.1 (10%).



• Constraint for Absolute Time Variation of Young’s Modulus. Limitation for the absolute change of Young’s modulus value.
  Default value: -1



• Constraint for Relative Time Variation of Young’s Modulus. Limitation for the relative change of Young’s modulus value.
  Default value: 0.1 (10%).



• Constraint for Absolute Time Variation of Poisson’s Ratio. Limitation for the absolute change of Poisson’s ratio value.
  Default value: -1



• Constraint for Relative Time Variation of Poisson’s Ratio. Limitation for the relative change of Poisson’s ratio value.
  Default value: 0.1 (10%).



• Constraint for Absolute Time Variation of Yield Stress. Limitation for the absolute change of the Yield Stress value.
  Default value: -1



• Constraint for Relative Time Variation of Yield Stress. Limitation for the relative change of the Yield Stress value.
  Default value: 0.1 (10%).



• Constraint for Absolute Time Variation of Hardening Coefficiernt. Limitation for the absolute change of the Hardening value.
  Default value: -1



• Constraint for Relative Time Variation of Hardening Coefficiernt. Limitation for the relative change of the Hardening value.
  Default value: 0.1 (10%).



• Constraint for Absolute Time Variation of Linear Expansion. Limitation for the absolute change of the thermal linear expansion value.
  Default value: -1



• Constraint for Relative Time Variation of Linear Expansion. Limitation for the relative change of the thermal linear expansion value.
  Default value: 0.1 (10%).

Contact Modeling Tab

• Penalty Magnitude. Parameter sets the value of the penalty multiplier in the method of penalty functions. Increasing the value of the multiplier in a number of cases makes it possible to increase the accuracy of the modelling of the contact interaction, but with increasing this parameter, the simulation time increases (the number of iterations necessary to achieve Newton's convergence increases).
  Default value: 5.0 (is recommended).

Stress Simulation at subsequent technological stages

The stress-strain state of the casting at the end of solidification is due, in general, to the history of the temperature field change of the casting and its contact with the mold. This condition is not final, and can significantly differ from the stress simulation of the casting demolded and cleaned from the elements of the gating-feed system (GFS). Thus, in order to understand what the final geometry of the casting will be and how the residual stresses will be distributed in it, it is necessary to simulate the stress after removing the casting block from the mold and then after removing the GFS elements. The Hooke solver provides necessary functionality for this.

The parameters responsible for the stress simulation at the technological stages, after the casting has solidified, are set on the Stress tab. The Simulate Stress after Mold and Gating System Removal parameter enables the additional Stress simulation at the stages of demolding and removal of the gating system.

The extraction of the casting block from the mold is carried out at a certain point in time or at a certain temperature. For the correct modeling of the stress during the demolding, the dynamics of the change in the temperature field of the casting is important first of all. Therefore, the simulation closest to the actual situation can be performed if the *.cst temperature field file simulated in the Fourier solver contains a complete history of temperature changes in the casting nodes from the beginning of solidification to the point of extraction from the mold. However, in the most cases, the demolding of the casting block from the mold occurs at the temperatures that are well below the solidus temperature, and sometimes at ambient temperature. Then, the simulation of the cooling down to the room temperature will take a long time, the casting temperature files (*.cst), the mold (*.mld) and others will be quite large. To avoid such a situation, the Mold Removal Temperature parameter sets the average temperature of the casting when it is removed from the mold. Using this parameter allows you to do the solidification simulation in the Fourier solver as usual, to the solidus temperature. When simulating in the Hooke solver, first the stress of the casting (with or without the mold) is simulated by the temperature fields from the *.cst file (until the last moment of time), and then the stress of casting without a mold is simulated at the temperature specified by the indicated parameter. The task of searching the temperature fields and the time required to achieve this average temperature is not considered here. If to set the parameter Mold Removal Temperature to a negative value (for example, -1), the simulation of the stress of the casting block upon demolding will be made on the basis of the temperature field saved in the last simulation step in the *.cst file.

At this stage the results of the stress simulation are saved in a separate file *_MLDREM.nds. If the creation of the additional files is enabled on the "Show the additional fields" tab in the solver settings, they will also be created when simulating this stage. The names of the additional files will have a "_MLDREM_" prefix.

Note: Stress simulation after the demolding of the block from the mold is carried out only if the contact interaction of the casting with the mold is switched on. With the contact interaction turned off, the stress simulation after the block is removed from the mold is not performed.

Cutting off elements of the gating-feed system usually occurs at ambient temperature. However, the user can set a different average casting temperature for this operation using the Gating System Removal Temperature parameter. The temperature of the casting varies step-like from the temperature of demolding to the temperature of the gating-feed system cutting of. The case of searching the temperature fields and the time required to reach this average temperature is not considered here. If to set the Gating System Removal Temperature parameter to a negative value (for example, -1), the simulation of the stress will be made on the basis of the temperature field saved in the last simulation step in the *.cst file.

The separation order of the gating-feed elements from the casting is specified by the Gating System Elements Removal Order parameter. The indices of the separated volumes are specified via a space in the required order. It means that the gating-feed elements, the removal of which is planned to be simulated, should be assigned the separate volume indices in advance (except 1). The casting itself must have a volume index of 1. This is done in the Master module before simulating the temperature field in the Fourier solver. If the parameter is specified by an empty line, the stress simulation after gating-feed system elements cutting off will not be performed.

At this stage, the results of the NDS simulation are saved in a separate file *_REMx.NDS, where x is the indices of the deleted volume. If writing additional fields to *.u3d files is enabled, they will also be created when simulating this stage for each deleted volume. The names of the additional files will have a prefix "_ REMx _".

Note: The files containing the gap field between the casting and the mold (*_gap.u3d and *_gap_mm.u3d) are not recorded when simulating the stress at the stage of mold and gating-feed system elements removal.

About the method of analysis of stress simulation

This section contains the brief data of the simulation methods of stresses and deformations used in the Hooke solver. More detailed can be found in the followings publications:

• Alexandrovich A.I., Monastirsky A.V., Solovjev M.B. Simulation of the thermo-elastic-plastic strength condition of the 3D parts // Report for applied mathematics of computation center (CC) named after Dorodnitsin A.A., the RAS, М., 2006, 31 p.
• Monastirsky V.P., Alexandrovich A.I., Monastirsky A.V., Solovjev M.B., Tihomirov M.D. Simulation of stress-strain state of the casting during the crystallization. Foundry.- 2007, № 8, p.45-47.
• Monastirsky A.V., Alexandrovich A.I., Solovjev M.B. Simulation of the stresses and deformations in the cast sections during crystallization // The proceedings of the XV International conference on computational mechanics and modern applied software systems (CMMASS’2007) May 25-31, 2007, Alushta. – M.: Vuzkniga, 2007. – 544 p.
• Monastirsky A.V., Smikov A.F., Pankratov V.A., Alexandrovich A.I., Solovjev M.B. Estimation of creation of hot cracks and simulation of the casting twisting in the CSS «PoligonSoft». Foundry.- 2009, № 5, p. 46-49.

and other.

The simulation of stress, appearing at the moment when the casting is cooled, is made on the basis of the theory of small elastoplastic deformations created by AA. Ilyushin (see, for example, [1]) using the linear finite element (FE) approximation on the unstructured tetrahedral meshes. The FE approximation of the elastoplastic task leads, at each step of the simulation, to a system of nonlinear equations with respect to the components of the nodal movements. To solve this system, the Hooke solver uses the Newton's method (for example, see N.S. Bahvalov, N.P. Zhidkov, G.M. Kobelkov. Numerical methods – third edition. – M.: BINOM. The laboratory of the knowledge, 2004).

The Newton’s method has a high rate of iterations convergence, but it requires a good initial approximation to the solution. In case if the initial approximation is far from the solution, the convergence may be violated (moreover, the iterative process may begin to diverge, which will be accompanied by a catastrophic increase in the error of the solution).

At each step of the simulation, the movement field, simulated in the previous step, is taken as the initial approximation of Newton's method. Thus, the quality of the initial approximation will be determined by the "proximity" of the movement fields corresponding to the current and previous steps of the simulation. This is crucially influenced by the rate of change in time of the input data of the task (the temperature fields and physical-mechanical properties of the environment) and the magnitude of the time step. When you manually select the value (Use Variable Time Step parameter is disabled) of this step, you should guide how intensively these data change over time. Using the variable time step (Use Variable Time Step parameter is enabled) allows to automate this process: in this case, the step is selected automatically according to the set restrictions on the values of the time variation of the input data of the task (parameters * _atol and * _rtol). In the case if the magnitude of the time step is set too large, or the limits of the rate of change in the time of the input data of the task are set unsuccessfully, the convergence of the Newton’s method may be violated. The Hooke solver indicates that the application of the Newton’s method is failed, if either the magnitude of the discrepancy value of the system increased monotonically on three consecutive iterations, or if the required accuracy of the solution (set by the Required Accuracy of Stress-strain State Calculation parameter) was not reached after the maximum allowed number of iterations of the Newton’s method (set by the Limit of the Newton's Method Iterations parameter). In the case of such a situation, the initial approximation for the Newton’s method at the next step of the simulation (which is taken as the movement field simulated at the current step) can be too rough, and the application of the Newton’s method may fail again. Thus, the "unsuccessful" application of the Newton's method at any step of the simulation, as a rule, leads to "spoilage" of the rest of the simulation. Therefore, after the first "failure" of the Newton's method, an emergency completion of the simulation is performed. After the emergency termination, you can, for example, specify the other simulation parameters (reduce the time step size, or set more stringent limits on the rate of change of input data) and try to restore it based on the data of the state file.

When simulating the stress for the initial time, a zero movement field is taken as the initial approximation of the Newton’s method, which corresponds to a homogeneous temperature field in the absence of the external loading and plastic deformations. For this reason, to ensure the convergence of the Newton’s method, it is desirable to start the simulation from a homogeneous temperature field (or close to that).

Note: Starting with the version 13.4, the Hooke solver uses the modified Newton’s method, which is characterized by higher stability and much less demanding for the quality of the initial approximations.

Nodes Fixation

When simulation the stress of the casting without taking into account the contact interaction of the casting and the mold, it is necessary to set the special boundary conditions for some mesh nodes - the fixing conditions (FC). The purpose of the setting the FC is to exclude the possibility of shifting the entire cast as an absolutely rigid body. Otherwise, the task turns out to be undetermined, and the movement field acquires a certain degree of arbitrariness (the corresponding elastic-plastic task has an infinite set of solutions). In particular, this can lead to the fact that at each step of the simulation the part will be moved or rotated in an unpredictable manner, or even to the impossibility of the stress simulation on individual simulation steps.

The FC setting can be performed either automatically or manually; the corresponding mode is determined by the value of the parameter Automatic model fixing. When automatic setting of FC (Automatic model fixing parameter is enabled), the number and spatial configuration of the symmetry planes are taken into account of the task under consideration. The constraining of the additional restrictions on any nodes of the casting is not allowed.

When manually setting the FC nodes (Automatic model fixing parameter is disabled), the correct setting of these conditions is entirely the user's responsibility, since an automated check of the correctness of the setting of these conditions is not possible. When specifying these conditions, it is necessary to take into account the number and spatial arrangement of the symmetry planes of the task under consideration (see below). However, it should be noted that in the presence of the symmetry planes that are not parallel to the coordinate planes, in general, it is impossible to correctly determine the FC task, since in this case it is necessary to limit the movements of the individual model nodes along the directions other than the coordinate axes.

Setting of the FC including with the use of the automatic fixation of the model in space, is not allowed in the stress simulation considering the contact interaction of the casting and the mold. This is due to the fact that the FC implied on the casting can conflict with the geometric constraints on the contact boundaries, in the result of which the simulated fields can acquire the nonphysical distortions.

Important! The nodes fixing that are not on the surface of the body is not allowed.

The conditions for the nodes fixation are specified in the Mirage module and saved to the mcst file (see the section Creating the boundary conditions files for the stress simulation). The number of the fixed nodes cannot exceed five thousand.

In case if the considered thermoelastoplastic task does not have any plane of symmetry, then the movement field is determined to within an arbitrary movement of the entire model as an absolutely rigid body, that is, the model has six degrees of freedom. For this case, there is a traditional scheme of imposing the restrictions on 3 casting nodes (see the figure below), the application of which completely excludes the arbitrariness in determining the movement field. Its meaning is that, despite the fact that all six degrees of freedom of casting are limited, there is no obstacle to the movement of its individual nodes (with the exception of one node). As a result, a reliable stress field without peak values at the fixing nodes is obtained.

Scheme of fixing the casting nodes

Consider the proposed scheme. The first node is prohibited from moving across all three axes. This restriction fixes three degrees of freedom of the casting (moving along the axes x, y, z). The second node is moved relative to the first along the x-axis (the orientation of the fixed nodes relative to the coordinate axes is important), it is prohibited to move along the vertical axis (fixing the casting rotation around the z-axis) and along the axis perpendicular to the segment node1-node2 (prohibition of the casting rotation around the vertical axis). At that, the second node can move freely along the x-axis. And, finally, the third node fixes the last degree of freedom - prohibits the rotation of the casting around the x axis.

Let us consider the possible schemes of a "manual" setting of FC in the case of the presence of symmetry planes. It is necessary to distinguish the following cases:

• One symmetry plane,
• Two intersecting symmetry planes,
• Three intersecting symmetry planes.

These cases are characterized by a different number of the degrees of freedom, and, as a consequence, the various schemes for fixing the nodes.

The case of one symmetry plane. For the clearness let the plane of the symmetry be the plane y = 0 (see the figure below). Then, due to the symmetry of the movement field with respect to this plane, the component uy of the movement field becomes zero on this plane. This circumstance excludes the rotations around any axes parallel to the coordinate axes Oz and Ox, as well as moving the entire body as the absolutely rigid along the direction of the Oy axis. Thus, the number of the degrees of freedom is reduced to three: the movements of the body are possible as the absolutely rigid in the directions of the axes Oz and Ox (i.e. parallel to the plane of the symmetry) and the rotations of the body around a certain (arbitrarily chosen) axis parallel to the Oy axis. It is possible to fix these degrees of freedom by imposing the following additional conditions on the movement field. Let’s take an arbitrary segment AB parallel to the Ox axis (see the figure below). At the A node, we set uz = ux = 0 (which fixes the translational degrees of freedom), and at the B node - uz = 0 (which fixes the remaining rotational degree of freedom). We can achieve the same result in the following way: taking an arbitrary segment CD parallel to the Oz axis, at the C node we require uz = ux = 0, and at the D node - ux = 0.

Scheme of fixing with one symmetry plane

The case with two intersecting symmetry planes.For clearness, let the planes z = 0 and x = 0 be the planes of symmetry. In this case, under the conditions: uz = 0 for z = 0 and ux = 0 for x = 0, only one degree of freedom remains: the movement of the body as the absolutely rigid in the direction of the line of intersection of the planes of symmetry (that is, the Oy axis). In order to fix this degree of freedom, it is enough to put uy = 0 in one arbitrarily chosen node. In the general case, that is, when there are two intersecting symmetry planes oriented in a space in an arbitrary manner, it is necessary to require that the projection of the movement vector to the direction of the intersection line of these planes vanishes at some nodes.

The case with three intersecting symmetry planes. There are possible two situations:

1. The lines of the intersection of these planes are parallel to each other;
2. Among the intersection lines of these planes there are the intersecting.

In fact the case a) does not differ from the case of two intersecting planes of symmetry already considered. In the case b) all six degrees of freedom of the body are fixed, and no additional conditions imposing on the movement field is not required.

Stabilized Newton’s method

The possibility of the casting movement as an absolutely rigid body, as a rule, does not lead to a significant loss of the quality of the solution in simulation taking into account the contact interaction due to the presence of the geometric constraints on the contact boundaries. Nevertheless, the arising nonphysical movement because of the incomplete certainty of the problem may turn out to be significant in some cases. The procedure for fixing the casting described in the previous section cannot be used to eliminate such movements in the simulations taking into account the contact interaction, since the conditions of fixing imposed on the casting can conflict with the geometric constraints on the contact boundaries, and as a result, the stress simulated on the casting can acquire an unphysical character.

Completely eliminate the occurrence of such movements in the simulations, taking into account the contact interaction, allows the use of a special method for solving equilibrium equations - the stabilized Newton’s method. This method requires about twice as much memory and on average 2-3 times as much time for simulations as the Newton's standard method, and therefore it is recommended to use it only when significant non-physical movements of the casting occurs as an absolutely rigid body. In other cases it is recommended to use the standard Newton’s method.

Macrostructure Simulation

Computational Domain

The casting macrostructure is modeled on a 2D mesh of a cellular automaton (CA) with a square cell.

The simulation uses temperature fields calculated in the Fourier module on a finite element mesh (FE). To represent the temperature fields, the CA cells are grouped into blocks of the same size. The choice of block sizes depends on the temperature gradient in the selected section of the casting. Usually, the block size is at least 10x10 cells.

The nucleation and growth of grains occurs in the casting mushy zone. The dimensions and configuration of the mushy zone change over time. Therefore, in the modeling of the macrostructure at each moment of time, only a part of the CA blocks participate, which in the course of the calculation can become active for a while and then be excluded from the calculation.

To organize this process, CA blocks are united into local computational domains (LCD), which are connected to the calculation and removed from the calculation, depending on whether they have active blocks.

Initial Data of the Thermal Problem

The casting macrostructure is modeled using temperature fields calculated in the Fourier thermal solver on a FE mesh. When the structure solver is launched, the temperature fields are processed and transferred from the finite element mesh to the cellular automaton mesh. As a result of this processing, the temperature distributions for the selected section of the casting are extracted from the *.cst file, and a temperature history file (*.tbl) is formed, recorded for the CA blocks (see a, b in the figure below).

At the stage of creating a tbl-file, a global computational domain (GCD) is formed, which includes the selected area of the casting.

Simultaneously with the temperature history file, a mask file (*.flags) is created containing information about the boundaries of the area where the casting alloy is located. The file marks local computational domains (LCD) containing metal and LCD, located outside the casting and not participating in the calculation, as well as boundary blocks located at the of casting-environment and casting-mold boundaries . Based on these data, a file (*_temper_loc.fld) is created, which contains fragments of the temperature distribution for LCD that are active at a given time (see in the figure с below).

a	b

c

Thermal problem preparation : a) temperature field calculation in the Fourier solver; b) selection of the casting section and determination of the boundaries of the region (GCD) for the macrostructure modeling ; c) splitting the GCD into local computational domains.
1 – casting; 2 – melt (liquid phase); 3 – casting mushy zone; 4 – mold; 5 – global computational domain (GCD); 6 – boundaries of one of the local computational domains (LCD); 7 – active LCD, where there is metal in a solid-liquid state; 8 – inactive LCD, the metal temperature is higher than the liquidus temperature, or the metal has already solidified

Dividing the selected section of the casting into local computational domains allows optimizing computational resources and reducing the computation time. LCD become active when the temperature in at least one of the blocks drops below the liquidus temperature. LCD become inactive when the temperature in all blocks becomes less than the solidus temperature.

The number and dimensions of LCA and the sequence of their activation and deactivation are automatically determined by the solver based on the geometry and the selected casting section.

Basic Model Provisions

As a rule, metal alloys in the solid state have a body-centered or face-centered cubic crystal cell.

Initially, the solid phase nucleus has a shape close to spherical, however, the anisotropy of the growth of the crystalline material very quickly affects and the nucleus acquires the shape of a cube. In the coordinate system associated with the nucleus and oriented along the edges of the cubic crystal cell, these growth directions (i.e., the normals to the faces) are denoted as [100], [010], [001] (see Fig. a below). In case of the cubic crystal cell, the faces of the cube have the maximum growth rate. All these directions are equal, so we say that a nucleus with a cubic crystal cell grows along directions of the [001] type.

The growth anisotropy leads to the fact that the cubic nucleus very soon takes the form of a dendrite. A freely growing dendrite in a supercooled melt with a uniform temperature field has the shape of a three-dimensional cross with six primary arms (trunks) of growth. As the primary arms grow, their surface loses stability and secondary branches form on it. They also grow in directions like [001], i.e. perpendicular to the primary arms (see fig. b below).

a	b

Schematic representation of a solid phase nucleus (a) and a globular grain growing in a uniform temperature field (b). The coordinate system (X, Y, Z) is linked to the crystal cell

In the 2D approximation, a grain growing in a uniform temperature field has four first order arms of the [01] type of equal length. On fig. below, the grain is shown in the (X, Y) coordinate system associated with the casting. The grain has a random orientation in space, which is characterized by the angle a between the nearest vector [01] and the Y axis. The lines connecting the vertices of the primary arms form an encircling figure that outlines the boundaries of the space occupied by the grain. In a uniform temperature field, all primary axes grow at the same rate and, therefore, the grain enclosure is a square. The correct shape of the enclosure is realized only under conditions of a uniform temperature field. The presence of a temperature gradient leads to different growth rates of primary dendrite arms, which distorts grain boundaries. With the growth of many grains, their growth stops when they collide with neighbors, which also violates the correct shape of the grain.

Schematic representation of a solid phase grain in 2D approximation. The coordinate system (X, Y) is linked to the casting. Angle a characterizes the crystallographic orientation relative to the casting.
1 – liquid hase; 2 – solid phase; 3 – enclosure

Grain Growth Model

To simulate the growth of a solid-phase nucleus with a cubic crystal cell, the cellular automaton method is used. It is well known that CA have anisotropy of properties. The anisotropy of CA leads to the fact that, if special measures are not taken, cubic nuclei with an arbitrary orientation lose their original orientation after three generations.

The model implements a growth algorithm that does not require grain shape adjustment. The cellular automaton rules are based on the following principle. Having built a dendrite environment of scale t₁, on the field of the cellular automaton, at the next step we can obtain a similar figure of scale t₂ = t₁ + Dt by constructing a set of figures of scale Dt at each point of the current figure. n this way, consistent grain growth from smaller to larger scales can be ensured. The current time is used as the scale, and Dt is the time step for solving the growth problem.

The figure below shows the algorithm for changing grain size A of a square shape rotated at some angle relative to the CA grid in a uniform temperature field. The gray color shows the CA cells located within the boundaries of grain A, i.e. inside the line marked as t₁ = const. Since t₁ is the current process time, the grain boundary is a constant time line. Let us assume that at the next time step, under conditions of a uniform temperature field, the grain boundaries should represent a similar figure - a square with a half-diagonal, the length of which increased by segment L. The length of this segment depends on the growth rate of the primary dendrite axes and the time step Dt. The new boundary t₂ = const is obtained by constructing a set of new type B grains on the boundary gray cells of the CA (only some grains В₁ – В₅ are shown in the figure below).

For the correct display on the CA mesh of figures of scale Dt (type B), it is necessary that the number of cells representing these figures on the CA mesh be sufficiently large.

The optimal time step for solving the growth problem is defined as:

Algorithm for growth of grain of arbitrary orientation in a uniform temperature field

Empirically, it was found that the shape and orientation of the square when expanding its boundaries is not lost if the length of the segment L is 4-5 cells. On fig. below the result of simulation of the growth of two grains of arbitrary orientation in a uniform temperature field is shown, demonstrating the absence of dependence of the simulation results on the CA mesh.

Growth of two grains in a uniform temperature field at different times τ₁, τ₂, τ₃

The growth rate of the primary dendrite arms depends on the current undercooling of the melt:

Grain Nucleation Model

This model assumes heterogeneous nucleation of grains, both on the surface of the mold and in the volume of the melt.

The solid phase nuclei formation occurs almost instantaneously when critical undercooling is reached. Random deviations in the nucleation process are described by the Gaussian distribution:

In the interval the nucleus formation probability is p > 0.68.

Distribution of the solid phase nuclei during nucleation in the volume of the melt and at the “casting-mold” boundary

The maximum density of nuclei is limited by the size of the CA mesh. It is assumed that only one solid-phase nucleus can form in a CA cell. Therefore, the density of nuclei that can potentially arise in the volume of the melt must satisfy the condition , where h – is a CA cell size. Accordingly, at the “casting-mold” boundary .

Gaussian distribution characteristics (3) - average undercooling DT_max, standard  deviation DT_s and nucleus density n_max an only be determined experimentally. It should be taken into account that the experimental data on n_max can refer to a three-dimensional casting. The values for the 2D section of the casting are obtained from the following stereological relations:

Macrostructure Simulation Parameters

The parameters required to run the solver are on the Structure tab (see figure below). The parameters for the grain growth model are also required (see Equation 2 in the Grain Growth Model section). They are set in the material editor (see the Grain Growth Tab section for more details).

• Enable Macrostructure Model. The parameter enables the calculation of the macrostructure. If the parameter is disabled, the macrostructure calculation is not performed and the parameters of the Structure tab are ignored. Option disabled by default.

Below are the parameters required for the solver, they are divided into the following groups:

• Computational Domain and Time Range
• Model Parameters
• Local Computational Domains (LCD)
• Results Storage
• Solver Settings

Parameters for macrostructure calculation

Computational Domain and Time Range Group

• Computational Domain. Macrostructure calculation is performed on the casting temperature fields. Some casting parts may take a long time to solidify and, at the same time, they may not be part of the macrostructure calculation area. This parameter is used to limit the calculation time to the solidification time of the desired area only. To set the boundaries of the simulation area, double-click the left mouse button in the parameter line. Then you should select the desired casting area. To do this, turn on the Elements by Frame mode and select the desired area (for more details, see the Editing Element Indices in the Arbitrary Areas section). The casting mesh elements of the selected area will be highlighted (see the fig. below). The current selection can be copied to the parameter line with the button located at the end of the line. After that, the parameter line will display the coordinate ranges of the selected area, and the area will be selected by the described box (see the figure below). If the parameter is not set, the temperature field of the entire casting will be considered.

a)	b)

Selection of casting area for macrostructure calculation:
a) selection of casting elements using a frame; b) macrostructure calculation area

• Section Plane. The calculation of the macrostructure is performed in the selected section of the casting. To set it, double-click the left mouse button in the parameter line. In this case, the model will switch to the section plane editing mode, and the parameter line will display the coordinates of the normal to the plane (in brackets) and the shift value (d). The position of the section plane is adjusted as usual (for more details, see the Section Plane section), while the direction of the normal, which usually determines the direction for the section of the model, does not matter in this procedure (see the figure below). The current position of the plane can be copied to the parameter line with the button located at the end of the line. The final position of the plane as a parameter value is performed by the Apply button at the end of the parameter line.

Selection of the model section for macrostructure calculation

• Block Size. The temperature fields in the calculation are set on the blocks. At the beginning of the calculation, the temperature field of the casting is transferred from the FE mesh to the cellular automaton cells. CA cells are grouped into blocks of n × n cells (see next parameter). Within the block, the temperature in all cells of the cellular automaton will be the same.
  Specified in micrometers.
  Default value: 150 µm.



• Number of Cells per Block side. CA cells are grouped into blocks of n × n яof a given size (see parameter above), so the number of cells in a block determines the size of the CA cell. The following considerations should be taken into account when choosing the size of the CA cell. In order to correctly display the shape and crystallographic orientation of a grain, it must be represented on the automaton CA grid by a sufficiently large number of cells. The macrostructure of the casting is largely the result of the competitive growth of columnar grains, during which grains with a predominant crystallographic orientation displace other grains with a less favorable orientation. For adequate simulation of competitive growth, the size of the CA cell should be commensurate with the distance between the secondary axes of the solid phase dendrites l₂. There are diagrams in the literature that relate l₂ to the melt cooling rate in the casting mashy zone. It must be taken into account that l₂ value at the moment of melt crystallization is significantly less than such diagrams show, since they are built for a completely solidified casting. For castings weighing 10-20 kg, crystallizing within minutes (not hours), the recommended CA cell size is 15-20 microns, therefore, for a block with a size of 150 microns, 10 cells per block side should be specified. Further, all recommended parameters of the problem are given for such castings.
  Default value: 10.

Usually, the Fourier thermal solver calculates the temperature field starting from the moment the mold is filled until the casting is completely solidified. To model the structure in a selected region (section), it is necessary to have a temperature field in this region from the moment the liquidus temperature is reached to the moment the solidus temperature is reached. This time range can be substantially narrower than the total solidification time of the entire casting. The correct choice of the time range significantly reduces the size of the saved files and the calculation time. If the macrostructure calculation is performed separately from the thermal calculation, the necessary values can be determined in the Mirage postprocessor by analyzing the temperature fields in the desired area of the casting.

• Modelling Start Time. The point in time from which the solver will start processing temperature fields and calculating the macrostructure.
  Specified in seconds.
  Default value: 0 sec.



• Modelling End Time. The point in time at which the solver has finished processing the temperature fields and calculating the macrostructure. If the calculation of the macrostructure is performed simultaneously with the solution of the thermal problem and there is no way to narrow the time range (it is not yet known), the value of the parameter “-1” should be set, which means the processing of the temperature field of the entire casting.
  Specified in seconds.
  Default value: -1.

Model Parameters Group

• Minimal Growth Front Advancement. In order not to lose the shape and crystallographic orientation of the grain, the grain growth at each time step should occur by capturing a sufficiently large number of cells. The parameter specifies the length of the half-diagonal L, (see the section above).
  Specified by the number of CA cells.
  Default value: 5 cells (recommended).

The appearance of solid phase nuclei in a supercooled melt is a random process. The number of nuclei that appear in the melt when the supercooling changes by a small d(ΔT) value is described by the Gaussian distribution (3). The parameters of solid phase nucleation are set separately for the melt volume and for the casting-mold interface.

• Maximum Density of Nuclei. The parameter (corresponds to n_max in equation (3)) sets the density of solid phase nuclei formed in the melt upon reaching undercooling many times greater than DT_max (see equation (3)).
  Specified in m^-2 for volume and in m^-1 for surface, respectively.
  Default value: 6∙10⁶ for volume and 3∙10³ for surface, respectively.



• Average Undercooling. The parameter determines the position of the maximum of the Gaussian curve (corresponds to DT_max in equation (3)).
  Specified in degrees.
  Default value: 6 for volume and 0.1 for surface, respectively.



• Standard Deviation. The parameter defines the peak width of the Gaussian curve (corresponding to DT_s in equation (3)).
  Specified in degrees..
  Default value: 0.5 for volume and 0.1 for surface, respectively.



• Integration Range Parameter. The parameter sets the melt undercooling range, at which the nucleation of the solid phase can occur. Outside this range, the appearance of solid phase nuclei is impossible. Thus, the subcooling range is a multiple of the standard deviation DT_s.
  Default value: 5.

Local Computational Domains (LCD) Group

In the presence of large temperature gradients, a high rate of melt cooling, and/or an unreasonably large time step for calculating the thermal problem, errors in the automatic processing of the LCD are possible. To prevent this situation, the user is given the opportunity to independently determine (extend) the temperature interval for the existence of LCD using the following parameters.

• LCD Lower Temperature Limit. The parameter defines the minimum temperature at which LCD is included in the calculation. If the temperature of all LCD blocks is less than the specified one, LCD ceases to participate in the calculation.
  Specified in Celsius Degrees.
  Recommended value: T_sol - 50°C.



• LCD Upper Temperature Limit. The parameter defines the maximum temperature at which LCD is involved in the calculation. If the temperature of at least one LCD block is less than the specified one, the LCD begins to participate in the calculation.
  Specified in Celsius Degrees.
  Recommended value: T_liq + 50°C.

Results Storage Group

• Store last step result only. During the calculation process, a file with the only last calculation step is created and periodically updated. The data in this file provides a quick visualization of the final casting structure. It is many times smaller than the file with the full calculation history that is recorded to view the dynamics of the formation of the macrostructure. If the parameter is set, then the history of macrostructure formation is not recorded. The option is enabled by default.



• Resulta storage frequency. Simulation results can be saved to a file with a user-defined interval. This file is used to visualize the dynamics of grain nucleation and growth in the selected casting section. Storage frequency is specified in the number of time steps. If the result recording step is equal to 1, all calculation steps are recorded. With a recording step of 2, every 2nd step is recorded, and so on.
  Specified in steps.
  Default value: 1.

Solver Settings Group

• Random Number Generation Method. The moment of grain nucleation and their crystallographic orientation are determined randomly. There are two ways to generate random numbers that determine this process.
  • Use System Timer. When using the system timer, the results of solving the same problem will differ somewhat in the number of grains, their shape, and crystallographic orientation. In this case, the type of structure, for example, the presence of columnar grains on the periphery and globular grains in the center of the casting, will be preserved.
  • Use Fixed Value. If a fixed value is used to run the random number generator, the results of solving the problem will be absolutely the same.
  Default value: Use System Timer.

Heat Treatment Simulation

The structural and mechanical properties models in the Heat Treatment module are suitable both for cast and deformable steels and can be used to simulate the following types of heat treatment: annealing, normalizing, quenching and tempering. The basic model is developed in the Creusot-Loire laboratory [1, 2], including:

• austenite decomposition model, which is used to calculate the distribution of structure and hardness after the part cooling;
• tempering model, which is used to calculate the hardness field as a result of tempering steel (calculation results based on the austenite decomposition model are used as input data).

1. Maynier P., Jungmann B., Dollet J. Creusot-Loire System for the Prediction of the Mechanical Properties of Low Alloy Steel Products. Hardenability Concepts with Applications to Steel. //The Metallurgical Society of AIME, 1978, pp. 518-545.
2. Maynier P., Dollet J., Bastein P. Influence des elements d'ailliage sur Ia trempabilite des aciers faiblement allies. //Revue de Metallrgie, Avril, 1970.

The basic model is supplemented with its own statistical formulas for tensile strength, yield stress and elongation (percentage elongation) prediction.

To calculate the structure and properties obtained after cooling from the austenitization temperature, the following input data are required:

• chemical composition of the alloy;
• time of isothermal exposure at the temperature of austenitization;
• calculation results of part cooling from austenitization temperature (i.e. heating temperature for normalizing, annealing or quenching).

Heat treatment modeling is performed in conjunction with the calculation of workpiece cooling according to a given mode. In this case, it is necessary to set the austenitization temperature as the initial one (the workpiece is considered to be heated evenly). The cooling rates required to predict heat treatment results will be calculated based on the data obtained from solving the thermal problem.

The result of the calculation using the austenite decomposition model will be the following set of data:

• cooling rate in the temperature range 800-500°C, °C/s;
• volume fraction of ferrite-pearlite mixture, %;
• volume fraction of bainite, %;
• volume fraction of martensite, %;
• Vickers hardness;
• yield stress, MPa;
• tensile strength, MPa;
• elongation (for five-fold samples), %.

To calculate properties after tempering the following data are required:

• chemical composition of the alloy;
• tempering temperature;
• holding time at tempering temperature.

The result of the calculation will be the following data set:

• Vickers hardness;
• yield stress, MPa;
• tensile strength, MPa;
• elongation (for five-fold samples), %.

The Heat Treatment model is integrated into the Fourier solver and can be used by the user when setting the appropriate calculation problems. The figure below shows a diagram of the heat treatment model operation coupled with thermal analysis performed by the Fourier solver.

Scheme of the Heat Treatment model functioning

Upon completion of the calculation, in addition to the temperature field file with the cst extension, universal format files with the u3d extension will be written, containing heat treatment simulation results. When simulating heat treatment with austenite transformation:

• *_ht_crs.u3d – cooling rate;
• *_ht_fp.u3d – ferrite-pearlite structure;
• *_ht_b.u3d – banite structure;
• *_ht_m.u3d – martensite structure;
• *_ht_hv.u3d – Vickers hardness;
• *_ht_ys.u3d – yield stress;
• *_ht_uts.u3d – tensile strength;
• *_ht_elong.u3d – elongation.

After further tempering (if provided):

• *_ht_t_hv.u3d – Vickers hardness after tempering;
• *_ht_t_ys.u3d – yield stress after tempering;
• *_ht_t_uts.u3d – tensile strength after tempering;
• *_ht_t_elong.u3d – elongation after tempering.

Limitations and Assumptions

Chemical Composition Limitations

Calculations can be made for low carbon, medium carbon and low alloy steels with alloy content according to the table below.

Heat Treatment Technologies Limitations

The models underlying the module are designed to analyze heat treatment processes, including heating to the austenitization temperature and subsequent cooling down to ambient temperature, as well as heating and holding for high tempering. Using the model, the structure and properties after annealing, normalization, quenching, quenching followed by tempering, some cases of surface quenching, etc. can be assessed.

The models is not intended to predict properties based on the results of heat treatment processes that include isothermal holding during the cooling stage or interrupted cooling. The module is not designed to simulate the properties generated by heat treatment processes such as step quenching, interrupted quenching, self-tempering quenching, and isothermal quenching.

The allowable range of tempering parameters is shown in the table below.

Limitations in the Differentiation of Structural Components

The structural components "ferrite" and "pearlite" in the results of the model calculation are not separated in any way and their separate ratios are not calculated. The model also does not provide for the separation of various morphologies of bainite.

Model Assumptions

At the beginning of the calculation, the part is considered to be heated uniformly: all mesh nodes have the same temperature and the same thermal history. The isothermal holding time at the austenitizing temperature is set by the user.

The calculation of the tensile strength is carried out according to known dependencies based on the obtained Vickers hardness values. Yield stress and elongation values are calculated using statistical models built on the b30asis of data from mechanical tests in production laboratories. The input parameters of these models are the tensile strength values and the chemical composition of the steel.

Model Accuracy

Below in fig. a comparison of the calculated and actual values of Vickers hardness according to the authors of the method [1] is presented: the measured value of hardness is plotted along the X-axis, the calculated value is plotted along the Y-axis. As can be seen, the method provides an acceptable accuracy of calculations.

Comparison of calculated and actual Vickers hardness values

Simulation Parameters

Since the heat treatment analysis uses the temperature history of the cooling of the workpiece, it is necessary to specify all the necessary data for the thermal calculation on the General and Solidification tabs in the Start Simulation and Parameters window (in this case, the name of the Solidification tab should be understood as “Thermal Calculation”). On the fig. below is the Shaft model prepared for heat treatment simulation. Cooling conditions are set as usual at the outer boundary. Depending on the mode, either the default boundary with index 17 (casting-environment) or any other combination of boundaries can be used, for example, if it is necessary to set a variable temperature of the quenching medium, designate air-cooled surface areas, etc.

Boundary conditions for Shaft quenching into oil simulation

• On the General tab, the following data is set (see the figure below):
  
  
  • Alloy Temperature. Sets the initial temperature value in Celsius degrees, the part cooling will be started.
  
  
  
  • Ambient Temperature. The temperature of the air or quench medium is set (when using boundaries 17 or/and 18).
  
  
  
  • Alloy properties. The part material is specified (imported from the database or its thermophysical properties are entered independently).
  
  
  
  • Heat Transfer Parameters. Heat transfer coefficients should be specified for all external boundaries. If necessary - the emissivity and the temperature of the medium (for more details, see the Heat Transfer Parameters section).
  
  
  
  • Temperature Modes (optional).

General Tab

• On the Solidification tab, the following data is set (see the figure below):
  
  
  • Calculation Step. The greater the expected cooling rate of the part, the smaller the calculation step should be, for example, 0.5 seconds. It is recommended to choose it in such a way that up to several hundred steps are performed during the cooling calculation.
  
  
  
  • Output Step. Should not be less than the calculation time step.
  
  
  
  • Simulation Stop Condition. Set the criterion for stopping the calculation Stop Temperature.
  
  
  
  • Stop Temperature. The parameter sets the temperature when the calculation will be stopped. Set the value to 500 °C or lower.

Solidification Tab

The parameters required to run the Heat Treatment solver are located on the Heat Treatment tab (see figure below).

Heat Treatment Tab

• Enable Austenite Decomposition Model. The parameter enables of the structure and properties model resulting from the decomposition of austenite during cooling. The parameter is disabled by default.



• Austenitization Time, h. The parameter sets the duration of holding the part before cooling.
  Specified in hours by decimal number.



• Enable Tempering Model. The parameter enables the property calculation model after tempering.



• Tempering Temperature, °C. The parameter sets the tempering temperature within the allowable range of 500-700°C.
  Specified in Celsius, as a decimal number.



• Tempering Time, h. The parameter sets the duration of isothermal holding at the tempering temperature in hours. The duration of exposure should be in the range of 1-200 hours.
  Specified in hours by decimal number.



• Composition, %. The table specifies the chemical composition of the alloy. In the corresponding lines, you need to enter the content of chemical elements in mass percent. Exceeding the specified range for each element is not recommended (in this case the model will work in extrapolation mode). If the alloy chemistry is set in the Composition tab of the Material Properties editor, it can be imported using the Extract Composition From Assigned Alloy button.

When click on the Advanced button, which is located under the data table, you go to the solver settings window (see the figure below). It contains the Use Adaptive Timestep algorithm parameter, which enables automatic control of the calculation step (enabled by default). In this case, the value of the Calculation Timestep parameter on the Solidification tab determines the maximum calculation timestep, which, if necessary, will automatically decreased as the calculation progresses to ensure that the result is obtained with sufficient accuracy. It is recommended to use the adaptive timestep in all cases, keeping the default setting.

Heat Treatment solver settings

Simulation and Results

To calculate the structure and properties, you should click on the Simulation button. The cooling of the workpiece during the heat treatment process and the resulting structure and properties will be calculated.

Note. The calculation of mechanical properties after tempering (if required) is always carried out in conjunction with the simulation of austenite decomposition and cannot be run separately.

As a result of calculations, u3d files with names consisting of the project name and abbreviated names of structural components and properties will appear in the project folder. These files (see figure below) can be opened in the Mirage postprocessor when using the universal format file load dialog.

Files with the results of the calculation of the structure and properties obtained during the analysis of the Shaft quenching

On fig. below some results of the calculation of oil quenching of a 600 mm shaft made of 40Cr steel and subsequent tempering at 550°C are shown.

The martensite structure distribution over the volume of the shaft - the result of the calculation according to the model of austenite decomposition (shaft_ht_m.u3d file)

Yield Stress distribution over the shaft volume - the result of the calculation according to the tempering model (shaft_ht_t_ys.u3d file)

Run of the Simulation

Running the Simulation

After the parameters of all solvers are specified, you can start the simulation. To do this, click the Simulation button in the lower right corner of the Start Simulation and Parameters window (see the figure below).

Running the simulation

Before the calculation is started, the preprocessor Master will require saving all the changes made in the model during the preparation of the calculation. The standard procedure for saving the project will be launched, where the user can specify the folder in which the following will be written:

• geometry file «.g3d»;
• DATA folder with project data and files required to run the solvers;
• DUMP folder for creating files needed to restore an interrupted calculation;
• solver source data configuration files Fourier (.cfo), Euler, (.cnh) and Hooke (.cnf);
• a file .bat with the same name as the geometry file and the calculation status file .status.

If a project that has already been run for calculation is overwritten on save, a warning message appears indicating that calculation results exist. (see figure below). If you answer Yes, the previous results will be overwritten. If you answer No, the save project will be canceled.

Warning about overwriting simulation results

After saving the model, the .bat file will be launched from the calculation folder (the folder where the .g3d file is saved), which is a script for launching the solvers. If desired, the user can edit it independently and start the calculation without the participation of the Master preprocessor.

After starting the calculation, the user sees the DOS console on the screen, which displays the progress of the calculation. To interrupt the calculation, just close the console.

During the calculation, the results are written to the project folder and can be viewed in the Mirage postprocessor (see the chapter Viewing Results).

Continuing an Interrupted Simulation

If the calculation was interrupted for any reason, it can be continued. When the solvers are running, their current state is constantly written to the DUMP folder, which is created in the project folder when the calculation is first started. When restoring, files from the DUMP folder are read and the calculation continues from the moment it was stopped. In this case, the results are appended to the previously created files.

To continue the interrupted calculation, in the Start Simulation and Parameters window, next to the Simulation button, there is the Continue Simulation option (see the figure below). If the calculation is started for the first time, the option is disabled and unavailable for the user. If the calculation was started and then interrupted, the option is enabled and can be used. When the calculation is started with the checked Continue Simulation, it will continue from the time point that was saved in temporary files in the DUMP folder.

Continue Simulation option

It should be note that when restoring an interrupted calculation, changes made to the project will be taken into account. The calculation will continue with new settings and parameters. It is not recommended to change the geometry of the model and the properties of materials, this can lead to malfunctioning solvers and / or unwanted results.

At startup, a message will be displayed (see the figure below) warning about possible changes to a previously started calculation. If you answer Yes, the calculation will continue. If the answer is No, the start of the calculation will be canceled.

Warning of changes in the project

If the Save Calculation State option is turned off in the Master module settings, continuation of an interrupted calculation is possible only from the start of the next solver (see details in Module Settings / General section). For example, if mold filling and then solidification and porosity are simulated, and the calculation is interrupted during the solidification stage, it can only be continued from the beginning of the solidification calculation. This limitation of functionality is indicated by the yellow icon next to the Continue Simulation option (see figure below).

Limited calculation continuation mode

If you unckeck the Continue Simulation option, then when you start the calculation, as usual, you will be prompted to save the project data and the results of the interrupted calculation will be overwritten.

New Pipe Shrinkage and Macroporosity Model

In the standard model available in the PoligonSoft, the creation of the pipe shrinkage and macroporosity is associated with the movement of the melt mirror. Moving the mirror is a consequence of feeding the mushy zone. In the standard model of macroporosity, it is assumed that the flow of liquid metal through the solid phase skeleton is unimpeded right up to the overlapping of the flow channels. At the moment when the melt mirror passes through this point of the casting, the formation of the porosity necessarily takes place, regardless of the magnitude of the external pressure.

In the new model, the formation of the porosity occurs at a time when the creation of a new section interface becomes energetically favorable, i. e. when the resultant pressure in the melt, formed from the external pressure, the metallostatic pressure, the capillary pressure and the pressure drop due to the shrinkage, is able to create a new pore. Thus, in the new model, the macroporosity prediction depends on the external pressure and the rate of metal cooling during crystallization. The cooling rate determines the structure of the casting (the grain size and parameters of the dendritic structure) and, thus, the capillary forces that play an important role in feeding the two-phase zone, the formation and growth of the pores.

It should be emphasized that the standard model of the porosity is an extreme case of a new model corresponding to low melt strength at break, high modulus of compressibility of liquid metal and zero surface tension. Under these conditions, the models demonstrate very close results, the differences in which are associated with the modernization of computational algorithms that ensure more accurate observance of the mass balance.

General Provisions

In general, the solidification process of the casting goes through the following stages: the formation of a pipe during the solidification of an open thermal node, the formation of a closed thermal node, the formation of an internal pipe shrinkage or dispersed porosity during the solidification of a closed thermal node (see the figure below).

A thermal node is considered open if the melt is in contact with the surrounding environment (Figure b). In the event that there is no direct contact of the melt with the surrounding environment, the thermal node is considered to be closed (Fig. c).

Schematic image of the casting solidification stages. М – melt mirror; L - melt; S – solid phase; F - mold. a) – initial conditions; b) – pipe shrinkage formation; c) – solidification of the closed volume of the melt in the thermal node with the creation of the internal pipe shrinkage or dispersed porosity.

Solidification of the Open Hot Spot

The solidification of the melt is accompanied by the shrinkage of the metal. In an open thermal node, the crystallization does not lead to a drop in pressure, since the shrinkage is compensated by lowering the free surface (or mirror) of the melt. The surface of the melt is able to move, i.e. is free if it does not have a fixed framework of the solid phase formed by the closed dendrites of the growing grains. As in the standard macroporosity model adopted in PoligonSoft, a fixed solid-phase framework appears when the fraction of the liquid phase decreases to a critical value f_L^**, i.e. at f_L ≤ f_L^**.

In the standard porosity model, the movement of the melt mirror in the stationary solid-phase framework leads to dehydration of the interdendritic spaces and the appearance of a porosity numerically equal to the fraction of the liquid phase f_L.

In a new model of the porosity, the melt mirror exists only where there is no solid-phase framework, and the interdendritic spaces above the melt mirror are filled with a melt that is retained there due to the capillary effect.

Pipe Shrinkage Formation In the Closed Hot Spote

At some point in time, because of the decrease in the fraction of the liquid phase, the boundaries of the thermal node become impermeable, and the thermal node becomes closed (above in the Fig. в). At this point, the free surface of the melt no longer exists in the thermal node and, consequently, the shrinkage of the metal during crystallization is no longer compensated by a change in the position of the mirror. This leads to a decrease in pressure in the thermal node, the intensity of which depends on the modulus of compressibility of the melt E, the shrinkage volume of the metal at a given time step, and the volume of the melt in the thermal node.

It should be noted that, the compressibility modulus of liquids is extremely large (the water is about 2000 MPa). Therefore, in the absence of communication with the environment, the crystallization of the metal leads to a rapid pressure drop in the thermal node to a critical value of P_crit.

If, at this point, there is a region of melt in the thermal node that is not bound by the fixed dendritic framework (i.e., the region where f_L > f_L^**), a new free surface of the melt appears, which subsequently leads to the formation of a closed pipe shrinkage. If there is no such region, the scattered macroporosity is formed.

Macroporosity Formation In the Closed Hot Spote

In the event when a fixed framework of the solid phase is present everywhere in a closed thermal node, i.e. everywhere f_L ≤ f_L^**, the formation of a flat free surface of the melt is impossible.

Breaking the melt integrity due to the pressure drop occurs by the formation of the macroporosity. The first pore occurs when the pressure drops below the critical value P_crit. As a rule, the point of origin of the pore is the point at which the pressure is minimal.

After the formation of pores began in the zone, the shrinkage that occurs at each time step is most likely to be compensated by the growth of already existing pores. The growth of the pore depends on the ratio of the capillary forces P_σ = 2σ/r, which tend to close the pore, and the pressure in the melt (see the figure below).

Formation of the macroporosity in the closed thermal node. а) – scheme of the pore creation in the fixed dendrite framwork; b) – the forces implied on the pore – melt section surface; L - melt; S – solid phase; – pore; – meniscus radius of curvature; - the distance between the first and the second dendrite axes .

The existence and growth of the pores depend on the equilibrium of forces on its surface, determined by the Laplace equation P_σ = -2σ/r (see the figure below). In this equation, P is the pressure in the melt; σ is the surface tension coefficient; r is the minimum radius of curvature of the pore, determined by the distance between the secondary branches of the dendrites λ_II, i.e. one of the parameters characterizing the framework of the solid phase.

Selecting Model Parameters

The formation of the shrink pore occurs when the pressure in the melt drops to a critical value P_crit, which is determined by the melt strength. In this model only shrink porosity is considered, for the formation of which a tensile stress must appeae in the liquid, i.e. the pressure must be negative.

Based on the surface tension for the nickel and iron melt σ ≈ 1.6 N/m, and the pore diameter observed in the castings 3.5÷ 60 µm, it is evident that the pressure of the melt at the moment of the pore formation can be from -0.1 to -1 MPa.

It is known that the theoretical strength of the melts is very high and can exceed 100 MPa. The actual strength of the melt is much lower and strongly depends on the purity of the metal and the way new interfaces are created, i.e. the pores. It can be assumed that the critical pressure at which the melt integrity is disturbed is within -0.01 - 1 MPa.

The meniscus radius r depends on the size of interdendritic spaces, which are determined by the rate of cooling of the melt in the two-phase zone. When solidification of the shaped castings weighing up to 100 kg, depending on the cooling rate of the melt the possible distance between the secondary dendrite axis constitute 20-100 μm. The distances between the primary axes, or between the grains, are 200 ÷ 600 μm.

The melt compression modulus can be estimated from the formula E = a²ρ_i, where a is the sound velocity in the melt; ρ_i is the melt density. Based on the known data on the speed of sound for pure metal melts, we can assume that Е = 3.104 ÷ 105 MPa. It should be taken into account that the sound velocity is much lower in the alloys than in pure metals.

In the model, the compression modulus characterizes the process of pressure drop in the thermal node. Under the ideal conditions, the rate of pressure drop in a closed thermal node is proportional to Е. During the crystallization of a real casting, the solid skin of metal surrounding the thermal node may be permeable, which reduces the rate of pressure drop in the node. There is also the possibility of deformation of the skin under the influence of the difference in the pressures of the environment and inside the thermal node, which also reduces the rate of pressure drop, since a part of the crystallization shrinkage is compensated by the deformation. This indicates that the effective compression modulus of the melt should be much less than the theoretical estimation possibly is equal to 200 ÷ 2000 MPa.

Within the frames of this model, the listed phenomena are not considered, and therefore the compression modulus, like the critical pressure, are the adjustable parameters that must be determined based on the experimental data.

To adjust the model, it is necessary to perform a series of simulations for the existing castings at different values of the parameters Е, P_crit, σ and λ_II and to compare the results with the porosity data. In this case, it makes sense to verify the adequacy of the specified thermophysical properties of the alloy and the mold materials, heat transfer conditions, etc., comparing the indications of the thermocouples installed on the casting block with the results of the thermal calculation.

Below there are some reference data on the surface tension for pure metals.