Shape Optimization of the Precision Die


The purpose of this simulation was to optimize the shape of a precision die. These kind of dies are used in extrusion process to create objects with fixed cross-sectional profiles such as window frames or bars, for example. The unmodified shape of the die is shown in the Figure 1a. Note that the outlet section has a varying thickness. Hence the velocity in the thinner section will be significantly lower than in the thicker one, as shown in Figure 1b. However, only the uniformed velocity at the die’s outlet warranties that the cast is straight. Hence one needs to modify the shape of the die so that more flow is directed to its thinner part.

Figure 1a. Undeformed shape of the die 

 Figure 1b. Velocity profile at outlet.

In order to perform the shape optimization process we used the QuickerSim CFD Toolbox for MATLAB in conjunction with the MATLAB Optimization Toolbox.


In the current case the flow of polystyrene (Styron 663) was simulated. Its properties are given in Table 1.

Property Value
Density [\frac{kg}{m^3}] 1040
Specific heat [\frac{J}{kgK}] 1200
Thermal conductivity  [\frac{W}{mK}] 0.1231
Thermal volumetric coefficient [\frac{m}{mK}] 6.60e-5



Table 1. Styron 663 properties

The polymer is a non-Newtonian fluid – we employed the Carreau-Yasuda rheological model,  where apparent fluid viscosity is defined by the following equation:

    \[\mu(\gamma)= \mu_\infty +( \mu_0 - \mu_\infty )[ 1+ (\lambda\gamma)^a ]^\frac{n-1}{a} ,\]

where \gamma denotes local shear rate and \mu_0 , \mu_\infty, \lambda , n, a denote model parameters describing fluid rheology. All coefficient values are given in Table 2. Note that the polymer exhibits a shear-thinning behaviour – viscosity drops with increasing shear rate. Also note that the value of asymptotic viscosity \mu_\infty is zero. Clearly, viscosity never drops to zero – that is just a mathematical model and the deformation rate in the extrusion process is limited.

Coefficient Value
\mu_0 [Pa \cdot s] 13 400
\mu_\infty [Pa \cdot s] 0
\lambda 0.527
a 0.845
n 0.351

Table 2. Carreau-Yasuda model  coefficients.

The variation of viscosity \mu with shear rate \gamma is presented in Figure 2.

Figure 2. Styron viscosity data.


Grids of different refinement levels were used to study covergence of the method. The number of points and elements in each mesh is given in Table 3.

Grid symbol Number of points Number of elements
1k 972 3440
2k 2152 7985
4k 4116 17722
7k 7024 32017
10k 10286 48912
20k 19831 99163

Table 3.

Impact of mesh parameters on results was investigated by means of comparison outlet velocity profiles along the middle line of the thicker cross section (the line that crosses the area of the highest speed), which is presented in Figure 3.

velocity profile

Figure 3. Velocity magnitude for different mesh sizes.

Note that the coarsest grid is the only one, for which the discrepancy is relatively high. For other cases the relative discrepancy doesn’t exceed 3 %. Since the time required to perform simulation significantly increases with the mesh size, the chosen number of nodes and elements were 2152 and 7985, respectively. The mesh was unstructured and the code was based on the Finite Element Method.

The shape of the undeformed geometry is shown in Figure 4.

Figure 4.  Undeformed geometry



To optimize the shape of the extruder the grid of control points was created. It contained three layers, all of which included the same number of points. The first and the last layer were contained in inlet and outlet plane, respectively and the middle one was at a height of 0.1 m. Each plane consisted of nine control points (see Figure 5).

Figure 5. Control points distribution.

They were translated only in horizontal planes so their z coordinates remained constant during optimization. Control points of the first and the last layer were not modified. Lagrange polynomials were used as interpolation functions in each plane. In order to interpolate points between layers linear polynomial were used due to technological aspects of manufacturing process (surface of the geometry should be ruled with straight lines). There are also other mesh deformation techniques like Radial Basis Function [see tutorial 6]

The objective function was defined as follows:


    \[ J(\textbf{x} ) = \int_S (v_z - v_{AVG})^2 ds,\]

where S denotes outlet surface v_z, v_{AVG denote z-component of the velocity vector and it’s average value, respectively and x denotes vector of decision variables. Vector \textbf{x} consists of translations of coordinates of the control points. Two types of constraints were applied. The first one is inequality with constant  values and the second one prevents from creating inconsistent mesh. Clearly, we wanted to minimize the discrepancy between the outflow-average velocity and the local one. Variation of the objective function value with gradient descent steps is shown in Figure 6.

Figure 6. Objective function convergence.


The outcome of the optimization procedure is the reduction of the cost function by 39%.  Simulation parameters are summarized in Table 4.  However, not the exact number is important – note that it is impossible to to reduce this function to zero – velocity will always be non-uniform because of the no-slip condition on the channel walls. The crucial fact is that the maximal velocity in the thinner section of the outlet is very similar to the one in the thicker case, what is clearly visible by comparing the outflow velocity profiles at Figures 7b and 8b.

Number of decision variables 18
Reduction of the cost function 39.6 %
Number of iterations 266

Table 4. Results

Figure 7a. and 7b. Original shape wiht velocity profile

Figure 8a. and 8b. Deformed shape with uniformed velocity profile.

This case also demonstrates that the Quickersim CFD Toolbox for MATLAB is a very flexible application for some CFD purposes. The strong point is a seamless integration with the MATLAB Optimization Toolbox.  Also due to simple data format for storing mesh, grid deformation is straightforward.

For other examples how  to combine CFD and optimization algorithms, see [Tutorial 22] of the QuickerSim CFD Toolbox for MATLAB.