Inverse Design Based on Nonlinear Thermoelastic Material Models Applied to Injection Molding

TL;DR

This paper introduces an inverse shape design method for thermoelastic bodies that determines initial shapes and temperature distributions from a known equilibrium shape, using nonlinear PDEs solved via finite element or isogeometric analysis.

Contribution

It presents a novel inverse design approach capable of estimating initial shapes and temperature fields for thermoelastic bodies, applicable to injection molding and other fields.

Findings

Successfully applied to inverse cavity design in injection molding

Capable of prescribing initial temperature and stress fields

Uses iterative fine-tuning for improved accuracy

Abstract

This paper describes an inverse shape design method for thermoelastic bodies. With a known equilibrium shape as input, the focus of this paper is the determination of the corresponding initial shape of a body undergoing thermal expansion or contraction, as well as nonlinear elastic deformations. A distinguishing feature of the described method lies in its capability to approximately prescribe an initial heterogeneous temperature distribution as well as an initial stress field even though the initial shape is unknown. At the core of the method, there is a system of nonlinear partial differential equations. They are discretized and solved with the finite element method or isogeometric analysis. In order to better integrate the method with application-oriented simulations, an iterative procedure is described that allows fine-tuning of the results. The method was motivated by an inverse cavity design problem in injection molding applications. Its use in this field is specifically highlighted, but the general description is kept independent of the application to simplify its adaptation to a wider range of use cases.

Figures

Tables

Table 1. Table 1: Forward simulation

	initial state	equilibrium state
prescribed quantities	1. boundary shape 2. internal stresses 3. temperature	1. temperature 2. external and internal forces
unknown quantities		1. boundary shape 2. internal stresses

Table 2. Table 2: Inverse simulation

	initial state	equilibrium state
prescribed quantities	1. internal stresses 2. temperature	1. boundary shape 2. temperature 3. external and internal forces
unknown quantities	1. boundary shape	1. internal stresses

Table 3. Table 3: Overview of stress measures

Symbol	relates forces in	to surface elements in
$𝝈$	equilibrium configuration $\overline{𝕩}$	equilibrium configuration $\overline{𝕩}$
$\widehat{𝐏}$	equilibrium configuration $\overline{𝕩}$	stress-free configuration $\widehat{𝕩}$
$\widehat{𝐒}$	stress-free configuration $\widehat{𝕩}$	stress-free configuration $\widehat{𝕩}$
$𝐏$	equilibrium configuration $\overline{𝕩}$	reference configuration $𝕩$
$𝐒$	reference configuration $𝕩$	reference configuration $𝕩$

Table 4. Table 4: Simulation parameters

Elasticity model	St. Venant-Kirchhoff
$\lambda$	$2\times {10}^6$
$\mu$	$0.5\times {10}^6$
$\overline{𝕗}=\rho 𝕘$	${(0,-2000)}^T$

Table 5. Table 5: Heat equation parameters

$\kappa$	$0.41$
Initial temperature	$0$
Boundary temperature	$-50$
$\Delta t$	$0.01$
Number of time steps	$5$

Table 6. Table 6: Elasticity parameters

Elasticity model	St. Venant-Kirchhoff
$\lambda$	$0.01$
$\mu$	$100.0$
$\alpha$	$1.0$
$\kappa$	$0.2$

Equations

$${\mathbf{\sigma }}^{\phi }(r)$$

$$-\nabla \cdot \mathbf{P}(r,u,\theta )-\bar{f}\bar{J}(u)$$

$$u$$

$$\mathbf{P}(r,u,\theta )n-\bar{h}\bar{J}(u)\left|{\bar{\mathbf{F}}}^{-T}(u)n\right|$$

$$-\nabla \cdot {\mathbf{\sigma }}^{\phi }(r)$$

$$r$$

$$-\nabla \cdot q(r,u,\theta )-\bar{g}\bar{J}(u)$$

$$\theta$$

$$q(r,u,\theta )\cdot n-\bar{k}\bar{J}(u)\left|{\bar{\mathbf{F}}}^{-T}(u)n\right|$$

$$-\nabla \cdot \mathbf{P}(r,u,\theta )-\bar{f}\bar{J}(u)$$

$$-\nabla \cdot {\mathbf{\sigma }}^{\phi }(r)$$

$$-\nabla \cdot q(r,u,\theta )-\bar{g}\bar{J}(u)$$

$$u$$

$$\mathbf{P}(r,u,\theta )n-\bar{h}\bar{J}(u)\left|{\bar{\mathbf{F}}}^{-T}(u)n\right|$$

$$r$$

$$\theta$$

$$q(r,u,\theta )\cdot n-\bar{k}\bar{J}(u)\left|{\bar{\mathbf{F}}}^{-T}(u)n\right|$$

$$u$$

$$\mathbf{P}(r,u,\theta )n-\bar{h}\bar{J}(u)\left|{\bar{\mathbf{F}}}^{-T}(u)n\right|$$

$$-\nabla \cdot \mathbf{P}(r,u,\theta )-\bar{f}\bar{J}(u)$$

$$-\nabla \cdot {\mathbf{\sigma }}^{\phi }(r)$$

$$-\nabla \cdot q(r,u,\theta )-\bar{g}\bar{J}(u)$$

$$u$$

$$\mathbf{P}(r,u,\theta )n-\bar{h}\bar{J}(u)\left|{\bar{\mathbf{F}}}^{-T}(u)n\right|$$

$$\theta$$

$$q(r,u,\theta )\cdot n-\bar{k}\bar{J}(u)\left|{\bar{\mathbf{F}}}^{-T}(u)n\right|$$

$${\mathcal{P}}_e:\hat{\Omega }\to {\Omega }_e,$$

$${\mathcal{P}}_e\in {\left[{\mathbb{P}}_1(\hat{\Omega })\right]}^d,$$

$$\lambda \left({\Omega }_e\cap {\Omega }_f\right)$$

$$\bigcup _{e\in E}{\bar{\Omega }}_e$$

$${\mathcal{I}}_{\text{lin}}=\left\{f\in C(\Omega )|f{|}_{{\Omega }_e}\circ {\mathcal{P}}_e\in {\mathbb{P}}_1(\hat{\Omega }),\forall e\in E\right\},$$

$${\mathcal{I}}_{\text{spline}}$$

$${\mathcal{S}}_u$$

$${\mathcal{S}}_r$$

$${\mathcal{S}}_{\theta }$$

$${\mathcal{T}}_{\mathbb{U}}$$

$${\mathcal{T}}_{\mathbb{R}}$$

$${\mathcal{T}}_{\Theta }$$

$$0=$$

$$+{\int }_{\Omega }\nabla {\mathbb{R}}^h:{\mathbf{\sigma }}^{\phi }(r^h)\text{d}\Omega$$

$$+{\int }_{\Omega }\left[\nabla {\Theta }^h\cdot q(r^h,u^h,{\theta }^h)-{\Theta }^h\bar{g}\bar{J}(u^h)\right]\text{d}\Omega -{\int }_{{\Gamma }_N^{\theta }}{\Theta }^h\bar{k}\bar{J}(u^h)\left|{\bar{\mathbf{F}}}^{-T}(u^h)n\right|\text{d}\Gamma ,$$

$${\bar{\mathcal{S}}}_u$$

Full text

††journal: Finite Elements in Analysis and Design

1 Introduction

One of the principal strengths of numerical simulations is that they can give insights into the behavior of products before they are actually built. These insights can even be useful in improving product design before the first prototype is produced. In a best-case scenario, numerical simulations could actually be used to find the perfect product design on the computer, without having to produce any prototypes.

In this paper, we deal with inverse design problems for thermoelastic materials with nonlinear elastic behavior. Our goal is to determine an initial shape of a body, when all other relevant parameters, as well as the body’s shape under equilibrium conditions, are known. These types of problems are also called shape optimization problems.

In order to solve these types of problems, we have developed a method that is general enough to be applicable to a broad range of applications, mostly in the field of production engineering. We will specifically highlight the method’s usefulness in the field of injection molding, for which it has been developed. However, in the way the method is presented, or with just minor adjustments, it can also be applied to processes such as high pressure die casting or additive manufacturing (e.g., selective laser melting).

Injection molding is a process where a liquid polymer melt is injected into a cavity where it is cooled down such that it solidifies [1]. Both the cooling and solidification cause a decrease in the specific volume of the material. Since this volume change happens inhomogeneously, the shape of the body changes such that the end result is not just a scaled version of the cavity shape. Figure 1 is a simple sketch of how a molding’s shape could deviate from the cavity shape.

Of course, we would like to be able to prescribe the molding’s shape exactly. Thus, what we are looking for is an inverse method that can yield a cavity shape that produces the correct molding shape. This is illustrated in Figure 2.

The injection molding process can be subdivided into several stages. In this work we will focus on the part of the process that happens after the material has solidified enough to be considered fully elastic. This means that our simulation starts from a point in time when there is a solid body that has a shape that is still identical to the cavity shape. This body has inhomogeneous temperature and stress distributions – these can for example be determined from a filling and solidification simulation or also from experiments. The body eventually cools down to the environment temperature, and the stresses relax to a point where the body is in equilibrium, i.e., the normal stresses on the boundary become zero. During this process, the body undergoes significant shape changes. It is the simulation of this process that we have to invert to determine a useful cavity shape from a prescribed molding shape. For a detailed description of how this method can be applied to injection molding and what other approaches are available, we refer to [2]. To treat the problem in a more general fashion, independent of the specific application of injection molding, we will now refer to the cavity shape as the body’s initial shape, and the molding shape as the equilibrium shape.

In general, when a forward simulation is available, i.e., a simulation that can predict the body’s equilibrium shape given an initial shape, this problem can be dealt with in the context of mathematical optimization. One further requirement for this is the availability of an objective function, i.e., a scalar measure that can describe the geometric differences between some output shape of the simulation and the desired equilibrium shape. This would make it possible to use, e.g., methods of quasi-Newton type for the optimization, such as BFGS [3], or derivative-free optimization methods. Optimization methods have been used in the past to solve problems very similar to ours [4, 5, 6, 7].

Optimization methods can be used as soon as a working forward simulation for the problem is available. This makes them very versatile, but they also present challenges. Most importantly, quasi-Newton methods can only be applied if the derivative of the objective function with respect to the design parameters is available. In shape optimization problems, these parameters describe the unknown shape. Therefore, their number can be quite large for complex shapes. In cases, such as ours, where the objective function depends on the result of a rather complex simulation, the calculation of the derivative is difficult. This means that it either has to be approximated, or derivative-free optimization methods have to be used. Both options are time-consuming for large numbers of design parameters.

The optimization problem we are dealing with is a very special type of shape optimization problem, where the objective function is also a shape. Additionally, in contrast to other optimization problems, where the goal is to minimize some objective function without prior knowledge of the simulation result or minimum value, we actually know the result that we wish to obtain. This means that what we are trying to achieve is an actual inversion of the simulation. Such an inverse simulation can be much more efficient than the aforementioned shape optimization methods. Furthermore, this can remove the need to parameterize the unknown shape — a step that limits the design space to less-than-optimal solutions.

For isothermal elasticity, inverse formulations have been researched already for some time. In 1967, Schield [8] showed that an inverse nonlinear elasticity problem, where the equations are solved for the reference rather than the deformed state, has the same general form as a regular (i.e., forward) elasticity problem. This was later extended by Carlson and Schield [9]. Govindjee et al. [10, 11] showed how the standard equilibrium equations can be re-parameterized to formulate the inverse problem. In this case, the Lagrangian viewpoint typical of elastostatics is exchanged by a Eulerian viewpoint, where the known equilibrium configuration is used as the reference and the unknown initial configuration is solved for. [12] is another example of a similar formulation.

As we will point out in Section 2, we require the reference configuration to be different from both the equilibrium and initial configurations. Such an ALE-type (Arbitrary Lagrangian-Eulerian) approach was used by Yamada in 1998, where this was applied to incompressible hyperelasticity [13].

Such inverse formulations have been used, among other things, for the optimization of turbine blade shapes [14, 15, 16], large beams [17], compliant mechanisms [18, 19, 20], and airplane wings [21].

In contrast to many inverse elasticity problems, the main source of deformations are, in our case, not external forces but temperature changes, i.e., thermal expansion or contraction. The notion of inverse formulations in thermoelasticity has appeared in literature before, although sometimes in a slightly different context. Dennis et al. have applied the term to problems where the initial state is known, but parts of the boundary conditions in the equilibrium state are unknown, i.e., neither displacements nor forces are known on parts of the boundary. This is made up for by over-specifying boundary conditions on other parts of the boundary, which requires both displacements and forces to be prescribed [22, 23]. This means that both Dirichlet and Neumann boundary conditions are applied simultaneously in the same place in some cases. Such an over-specification of Dirichlet and Neumann boundary conditions also happens in our case, but on the entire boundary, to allow the prescription of an equilibrium shape at the same time as zero normal stresses on the boundary. The major difference to the work of Dennis et al. lies in the fact that we use this over-specification to solve for displacements in the initial state, which creates the need for multiple displacement fields. This is explained in more detail in Section 3.3.

2 Solution Strategy

2.1 Description of the Inverse Problem

Our problem has a slightly different setup from many other inverse design problems in elastostatics. In many cases, two states of a body are considered. There is one state where the body is subjected to external loads and body forces such as gravity. In a second state, the body is considered without any forces acting on it. We have a distinct problem where a body transitions from a constrained to an unconstrained state. We will call the chronologically first state the initial state, and the latter state the equilibrium state, since the body will be in thermodynamic equilibrium.

Before we detail the setup of the inverse problem, we will first fully describe the corresponding forward simulation. In the initial state, we are dealing with an inhomogeneously heated body that is physically constrained. The state of the body is fully described by its shape, the field of internal stresses and the temperature field. In the equilibrium state, we have knowledge of the environment temperature, gravitational forces and external forces, which may be zero in many cases. The principal quantity of interest is the shape of the body in the equilibrium state. However, the field of internal stresses in this state is also needed, since the stresses need to stay below a certain threshold for the elastic material models to remain valid. Table 1 shows an overview of the prescribed and unknown quantities in both states.

The inverse problem results from the wish to prescribe the equilibrium shape of the body. One particularity of this inverse problem is the necessity to prescribe a temperature and stress distribution on an initial shape that is still to be computed. Table 2 shows the changed distribution of prescribed and unknown quantities for the inverse simulation.

2.2 Handling of Auxiliary Fields

In Table 2, the internal stresses and temperature are italicized since they require special treatment. It is difficult to prescribe, for instance, an inhomogeneous temperature field, if the shape of the body is unknown. A mechanism needs to be devised that can produce the correct temperature field for a certain shape. An important requirement of this mechanism is that it has to be integrated seamlessly with the system of equations for the inverse problem to ensure that it can still be solved in a monolithic, i.e., non-iterative, fashion.

To initialize the inverse simulation, the temperature and stress fields are determined for a reference configuration chosen by the user (cf. also Section 2.3): In injection molding, one might perform a filling simulation for a cavity resembling a scaled-up version of the product shape. These fields, from now on referred to as reference temperature and reference stress, will then serve as a basis for estimating the corresponding fields in arbitrary shapes; a step which is of course only necessary if the fields cannot be computed directly. For this estimation, we use a mechanism that is identical for both the temperature and stress field, and could be applied to other auxiliary fields. Therefore, we will only describe it for the temperature field.

The temperature field in the initial state entirely depends on the application that is considered. Even small changes in the boundary shape could mean an entirely different temperature field. For our application of injection molding, we make the assumption that small changes in the boundary shape will also cause small changes in the temperature field. Following this logic, we look at the difference between a certain shape and the reference shape, and try to translate this shape difference to a difference in the temperature distribution. The only requirement for this translation is that for identical shapes, the temperature fields should also be identical.

As long as we do not wish to incorporate any knowledge of how the temperature field is produced in the first place into the system of equations, we should always aim to keep the difference between the estimated and reference temperature fields as small as possible, as long as we fulfill the constraint of the adjusted boundary shape.

In order to achieve this, we borrow an idea from interface tracking methods in free-surface flow simulations. In such contexts, the simulation mesh needs to be updated to fit a changed boundary in a way that avoids damage to the mesh in terms of deteriorating simulation properties [24]. A popular method for this purpose is called EMUM, the Elastic Mesh Update Method [25]. In this method, an elasticity problem is solved with Dirichlet boundary conditions to achieve, in a sense, minimal movement of interior mesh nodes. We apply the same method to move temperature nodes such that they fit a new shape. This means that we assume no actual changes in temperature values, but only in their distribution. One should note a small imperfection in this method, which stems from the fact that tangential sliding of boundary nodes, even if this does not contribute to an actual shape change, will induce unwanted changes in the temperature distribution. We currently neglect this issue in favor of efficiency. It is possible to use a linear elasticity model for this purpose in cases where the deformation is small enough. In other cases, more sophisticated equations, such as nonlinear elasticity, should be considered. In our description, we will refer to the equation for the node movements as smoothing equation.

2.3 Reference Configuration

We have so far introduced three different states of a body: the initial state, the equilibrium state, and a reference state. We will now explain in more detail the connection between the configurations of the body in these different states. For a structural finite element simulation, we need to define one reference configuration. This defines the simulation mesh as well as the coordinate system used to formulate the differential equations. We will define the set of points that belong to the body in the reference configuration as $\Omega\subset\mathbb{R}^{3}$. This allows us to describe the other configurations using displacement fields ${\mathbb{d}}\in\Omega$ that store the displacement vectors of all material points.

In a regular, i.e., forward-in-time, structural simulation, the initial configuration may be chosen as the reference. This is convenient, since the initial configuration is known and the displacement field for the equilibrium configuration immediately translates to the unknown quantity that is solved for. This approach is referred to as the Lagrangian formulation.

If we consider an inverse structural simulation, where the equilibrium configuration is fully prescribed, the opposite approach can be applied. This means that the equilibrium configuration serves as the reference, and a displacement field for the initial configuration is calculated. This is called the Eulerian formulation. As we already mentioned, this approach has been used in many cases to formulate inverse design problems [8, 10].

For our simulation, we have already mentioned a reference configuration that is distinct from both the initial and the equilibrium configurations. This is the configuration for which the auxiliary fields, such as temperature and internal stresses, that are used to estimate these fields in the initial configuration, are provided. We introduce the symbols

${\mathbb{r}}$ for the displacement field describing the initial configuration, and

• 2.

${\mathbb{u}}$ for the displacement field describing the equilibrium configuration,

both defined on the reference configuration $\Omega$.

The aforementioned smoothing equation for the transformation of the auxiliary fields determines the interior of the displacement field ${\mathbb{r}}\in\Omega\setminus\partial\Omega$ from the boundary displacement ${\mathbb{r}}\in\partial\Omega$. For this reason, the interior values of the displacement field ${\mathbb{u}}\in\Omega\setminus\partial\Omega$ need to be left unspecified. This means that neither ${\mathbb{r}}\in\Omega$ nor ${\mathbb{u}}\in\Omega$ are fully prescribed, so a different reference must be used. This third approach is called ALE, the Arbitrary-Lagrangian-Eulerian formulation. In the context of inverse design problems, this formulation can be found in, e.g., [13].

Figure 3 shows a sketch of how the different configurations relate to each other and which symbols are used to describe these relations. We have introduced the additional symbol $\mathring{\mathbf{B}}$ to describe the deformation from a theoretical, stress-free configuration to the initial configuration. This field indirectly stores the initial stresses. It makes sense to keep track of the initial deformation rather than the stress resulting from it, since this can more easily be integrated with existing stress-strain relationships. Furthermore, we only require the symmetric tensor $\mathring{\mathbf{B}}$ instead of $\mathring{\mathbf{F}}$, since this suffices in connection with the material laws that we are using and since it allows us to use results from simulations that do not explicitly keep track of $\mathring{\mathbf{F}}$.

2.4 Iteration Scheme

The method for handling auxiliary fields that we presented in Section 2.2 only approximates these fields based on an assumption of minimal change. However, the true effects of the geometry changes that are proposed by the inverse method on these fields are unknown. If we intend to incorporate these effects, we have to rerun the simulations that yield these auxiliary fields with the changed geometry. Such a rerun will offer new initial data for the inverse method, so the initial shape can in turn be improved. Repetition of these steps leads to an iteration scheme, such as the following:

Select best guess for initial shape (e.g. desired shape). 2. 2.

Determine auxiliary fields for current shape. 3. 3.

Run forward simulation to determine the equilibrium shape. 4. 4.

If the quality of the equilibrium shape is good enough (or another stopping criterion, such as maximum number of iterations, is reached), stop the iteration. 5. 5.

Run inverse simulation to determine improved guess for initial shape. 6. 6.

Continue at Step 2.

This iteration scheme is illustrated in Figure 4. The correction step that is mentioned in the illustration could be, as previously mentioned, another simulation, or some other method to determine the auxiliary fields.

There is no guarantee that such an iteration will converge. This completely depends on how well the true changes in the auxiliary fields match with the assumptions made in the smoothing equation.

In the specific case of cavity shape design in injection molding, the correction step consists of a viscoelastic fluid simulation, including polymer solidification models, that can generate a temperature and stress field for a provided cavity shape.

3 Conservation Laws in Regular and Inverse Formulations

After having described the inverse problem, we will now describe the equations that are needed to solve it. These include the equilibrium equation for the elastic body as well as the smoothing equation for the initial displacement field. We will also show how these equations can be discretized to be solved using the finite element method or isogeometric analysis.

3.1 Smoothing Equation

The equation that is needed to find an approximation for the changes in the auxiliary fields should lead to minimal changes in the interior of the field to match the changes on the boundary. Since we only have information about the geometric changes of the boundary, we can also only make assumptions about geometric changes or displacements in the interior. As we pointed out in Section 2.2, we borrow the idea from the Elastic Mesh Update Method (EMUM) to use an elastic material law to determine the interior displacements.

To keep this description simple, we chose linear elasticity for the smoothing equation, which was sufficient for all of our test cases. However, more elaborate equations could be used in its place. We use the symbol $\varphi$ for quantities relating to the smoothing equation. The stress tensor $\boldsymbol{\sigma}^{\varphi}$ that is used in this context is defined as follows:

[TABLE]

with the virtual Lamé parameters $\lambda^{\varphi}$ and $\mu^{\varphi}$.

3.2 Forward simulation

We will now formulate the differential equations that are required to run a forward simulation of thermoelasticity with initial displacements. We use the solution space $\mathcal{S}=C^{2}(\Omega)$, with our reference domain $\Omega$. The spatial dimensionality will be referred to as $d$, such that $\mathcal{S}^{d}$ is a multi-dimensional solution space. The momentum conservation is ensured by equation

[TABLE]

with the first Piola-Kirchhoff stress tensor in the reference configuration $\mathbf{P}:\mathcal{S}^{d}\times\mathcal{S}^{d}\times\mathcal{S}\to\left[C^{1}(\Omega)\right]^{d\times d}$ as a function of the initial displacements ${\mathbb{r}}\in\mathcal{S}^{d}$, the equilibrium displacements ${\mathbb{u}}\in\mathcal{S}^{d}$, and the temperature $\theta\in\mathcal{S}$.

$\bar{{\mathbb{f}}}$ is a field of forces acting in the equilibrium configuration, which is defined in the reference configuration.

The Dirichlet and Neumann boundary conditions for this equation are given by

[TABLE]

where $\Gamma_{D}\cup\Gamma_{N}=\partial\Omega$, which is the full boundary of $\Omega$. ${\mathbb{u}}^{D}$ is a prescribed boundary displacement and $\bar{{\mathbb{h}}}$ is a force acting on the boundary in the equilibrium configuration.

The smoothing equation for ${\mathbb{r}}$ is in this case also a type of momentum conservation equation, but for the virtual smoothing stress $\boldsymbol{\sigma}^{\varphi}$:

[TABLE]

There are no Neumann boundary conditions for the smoothing equation, as the boundary shape is completely prescribed. Therefore, with the boundary displacements ${\mathbb{r}}^{D}$ given in the initial configuration, we have

[TABLE]

Stationary energy conservation is fulfilled by the heat conduction equation

[TABLE]

with heat flux vector in the reference configuration ${\mathbb{q}}:\mathcal{S}^{d}\times\mathcal{S}^{d}\times\mathcal{S}\to\left[C^{1}(\Omega)\right]^{d}$. $\bar{g}$ is a heat source field in the equilibrium configuration. The Dirichlet and Neumann boundary conditions,

[TABLE]

can be used to set either fixed temperature values in $\theta^{D}$ or a heat flux normal to the surface in the equilibrium configuration ($\bar{k}$). Analogous to the boundary conditions of the momentum equation, we have $\Gamma^{\theta}_{D}\cup\Gamma^{\theta}_{N}=\partial\Omega$.

In summary, our system of equations becomes

[TABLE]

3.3 Inverse simulation

In order to run an inverse simulation instead of the forward simulation, very few things need to be changed. Notably, the differential equations for the conservation laws, i.e., Equations (2), (5), and (7), stay in place. The boundary conditions for the heat equation, Equations (8) and (9), also remain unchanged. In order to switch from prescribing the boundary in the initial configuration to prescribing it in the equilibrium configuration, Equations (3), (4), and (6) are replaced by the new boundary conditions

[TABLE]

This specifies both Dirichlet and Neumann boundary conditions for the elasticity equations on the whole boundary. This over-specification makes it possible to under-specify boundary conditions for the smoothing equation. Notably, both ${\mathbb{u}}^{D}$ and $\bar{{\mathbb{h}}}$ now have to be available on the full boundary $\partial\Omega$. If the shape ${\mathbb{u}}^{D}$ is known for an unconstrained state, $\bar{{\mathbb{h}}}$ is zero on the full boundary and Equation (11) simplifies to $\mathbf{P}({\mathbb{r}},{\mathbb{u}},\theta){\mathbb{n}}=0$.

The full system of equations for the inverse formulation reads

[TABLE]

3.4 Discretized Form of Equations

We solve the equations using the finite element method. As a preparation, we need to find suitable interpolation spaces and formulate the discretized weak form of the equations. We will sketch the interpolation spaces for both the standard finite element method with linear basis functions as well as isogeometric analysis with B-spline or NURBS basis functions. Other spaces can also be used but have not been tested.

We will first construct the interpolation space for the finite element method with simplex P1 (linear) basis functions. We use a reference element defined by a compact set $\hat{\Omega}\subset\mathbb{R}^{d}$. For the individual elements $e\in E$, where $E$ is the discrete set of all element indices, we define a linear projector into physical space

[TABLE]

with the space of linear polynomials $\mathbb{P}_{1}$, where $\Omega_{e}\subseteq\Omega$, and especially

[TABLE]

where $\lambda(\square)$ is the Lebesgue measure. This means that the physical elements $\Omega_{e}$ do not overlap except on their boundaries, and when they are combined, they make up the full domain $\Omega$. The finite element interpolation space can be constructed using the projectors as

[TABLE]

which makes it a space of piecewise linear continuous functions on $\Omega$.

As an alternative, we define the spline interpolation space for isogeometric analysis as

[TABLE]

with a compact reference space $\tilde{\Omega}\subset\mathbb{R}^{d}$ and the spline basis functions $B_{i}:\tilde{\Omega}\to\mathbb{R}$, which could be, e.g., B-Spline or NURBS basis functions [26].

In the following, we will refer to the interpolation space simply as $\mathcal{I}$, assuming that one of the presented options is used. Based on this we define the solution spaces

[TABLE]

and the test spaces

[TABLE]

For simulations in forward mode, we obtain the following discretized weak formulation: Find ${\mathbb{u}}^{h}\in\mathcal{S}_{{\mathbb{u}}}$, ${\mathbb{r}}^{h}\in\mathcal{S}_{{\mathbb{r}}}$, and $\theta^{h}\in\mathcal{S}_{\theta}$, such that

[TABLE]

for all ${\mathbb{U}}^{h}\in\mathcal{T}_{{\mathbb{U}}}$, ${\mathbb{R}}^{h}\in\mathcal{T}_{{\mathbb{R}}}$, and $\Theta^{h}\in\mathcal{T}_{\Theta}$.

For the inverse simulation, we have to make changes to the solution spaces for the displacement fields:

[TABLE]

Additionally, the test space for the momentum conservation equation has to be adjusted to

[TABLE]

to incorporate the full domain including its boundary.

The discretized weak formulation of the inverse problem is defined as: Find ${\mathbb{u}}^{h}\in\bar{\mathcal{S}}_{{\mathbb{u}}}$, ${\mathbb{r}}^{h}\in\bar{\mathcal{S}}_{{\mathbb{r}}}$, and $\theta^{h}\in\mathcal{S}_{\theta}$, such that

[TABLE]

for all ${\mathbb{U}}^{h}\in\bar{\mathcal{T}}_{{\mathbb{U}}}$, ${\mathbb{R}}^{h}\in\mathcal{T}_{{\mathbb{R}}}$, and $\Theta^{h}\in\mathcal{T}_{\Theta}$.

When compared to the equation for the forward simulation in (24), Equation (28) seems almost identical. The only change is in the integral limits of the boundary term for the momentum conservation equation, which now contain the full boundary $\partial\Omega$. In cases where $\bar{h}\equiv 0$, this term is removed and the equations become identical.

4 Derivation of Constitutive Equations

So far, we have left the first Piola-Kirchhoff stress tensor $\mathbf{P}$ and the heat flux ${\mathbb{q}}$ in the reference configuration undefined. We will provide definitions for both of these quantities in the following sections.

We have already provided some definitions of deformation measures in Figure 3. In addition, we will now introduce the deformation gradient tensor $\mathbf{F}$, which describes the full deformation in the equilibrium configuration with respect to the stress-free configuration:

[TABLE]

Analogously, we define $J:=\text{det}\,\mathbf{F}$.

Due to the multiple coordinate systems that we use, we also have to define multiple stress measures. Table 3 gives an overview of the stress tensors that we use. The tensors $\hat{\mathbf{P}}$ and $\hat{\mathbf{S}}$ correspond to the first and second Piola-Kirchhoff stress tensors, respectively, as they are commonly used.

4.1 Constitutive Laws for the First Piola-Kirchhoff Stress

In many cases where nonlinear elasticity laws are used, the reference configuration coincides with the stress-free configuration, i.e., $\hat{\mathbf{P}}=\mathbf{P}$, as well as, $\hat{\mathbf{S}}=\mathbf{S}$. For us, this is not the case. Therefore, when we use constitutive laws that are defined with a stress-free configuration used as reference, we need to carry out a coordinate transformation to the reference configuration.

Based on St. Venant-Kirchhoff

The following is a simple constitutive law for thermoelasticity that is based on the St. Venant-Kirchhoff material:

[TABLE]

where the Lamé parameters $\lambda$ and $\mu$ control the elastic behavior and $\alpha$ controls the thermal expansion (cf. [27] for a derivation of this law). It uses a simple linear thermal expansion law. Such a law should only be used with small temperature variations or for demonstration purposes. For this constitutive law, we can derive the following expression in $\mathbf{P}$, which can be inserted into the presented differential equations111For a full derivation of this equation, see A.1.:

[TABLE]

with $\mathbf{B}=\bar{\mathbf{F}}\tilde{\mathbf{F}}^{-1}\mathring{\mathbf{B}}\tilde{\mathbf{F}}^{-T}\bar{\mathbf{F}}^{T}$.

Based on Neo-Hooke

We also use a constitutive law that is based on a Neo-Hooke material (cf., e.g., [28]), with additional thermal expansion terms. This is given as an expression that relates the Cauchy stress and the left Cauchy-Green deformation tensor $\mathbf{B}$:

[TABLE]

The elastic behavior is controlled through the material parameters $D_{1}$ and $C_{1}$, and the thermal expansion behavior through the parameter $\alpha$. The thermal expansion term has been chosen this way for consistency with the law based on St. Venant-Kirchhoff that was presented in the previous section. For the stress tensor $\mathbf{P}$ that we require for our differential equations, we obtain the following expression222For a full derivation of this equation, see A.2.:

[TABLE]

with $J=\bar{J}\tilde{J}^{-1}\mathring{J}$ and $\mathbf{B}=\bar{\mathbf{F}}\tilde{\mathbf{F}}^{-1}\mathring{\mathbf{B}}\tilde{\mathbf{F}}^{-T}\bar{\mathbf{F}}^{T}$.

4.2 Constitutive Law for the Heat Flux

We have formulated the differential equations for an arbitrary heat flux vector ${\mathbb{q}}({\mathbb{r}},{\mathbb{u}},\theta)$. We have done so to point out that multiple options of defining the heat flux are available. For instance, the heat flux could be chosen to be consistent with regular formulations where the heat flux is a linear function of the temperature gradient with respect to the stress-free or equilibrium coordinate system. However, the only definition that we have used is one where the heat flux is linear with respect to the temperature gradient in the reference configuration:

[TABLE]

with scalar thermal conductivity $\kappa$.

5 Numerical Examples

We describe two different test cases in this section. The first test case, described in Section 5.1, is meant to be a simple demonstration of the capabilities of the inverse formulation when it is applied to isothermal elasticity. In the second test case, Section 5.2, we employ all the concepts that were introduced in the previous chapters to solve a thermoelastic problem.

5.1 Elastic beam under gravitational forces

In this numerical example, we answer the question of how a clamped elastic beam would have to be shaped, such that it takes on a completely straight shape when subjected to gravity. The test case demonstrates the simplest manner how the inverse formulation could be used, since temperature is ignored and initial stresses are not used.

We consider a simplified two-dimensional beam. The simulation mesh in the reference configuration, along with boundary conditions for the forward simulation, is shown in Figure 5.

We use the Dirichlet boundary condition ${\mathbb{u}}=0$ on the left-hand side boundary for the equilibrium displacement field. This makes sure that the beam stays attached to the wall. Additionally, we use the Dirichlet boundary condition ${\mathbb{r}}=0$ on all boundaries for the initial displacement field. This implies that the smoothing equation will ensure that the initial and reference configurations are identical333The initial displacement field as well as the smoothing equation are not necessary in the forward simulation, but we use them anyway to highlight the differences between forward and inverse setups.. The values for the material properties and the gravitational force are shown in Table 4.

In the equilibrium configuration, the lower right-hand corner of the beam is displaced by $-8.95\times 10^{-3}$ in horizontal direction and $-62.79\times 10^{-3}$ in vertical direction. Figure 6 shows the stresses in the beam in the initial and the equilibrium configurations.

For the inverse problem, we remove most of the Dirichlet boundary conditions on the initial displacement field ${\mathbb{r}}$. We only keep the boundary condition of ${\mathbb{r}}=0$ on the left-hand side boundary. This makes sure that the shape can be attached to a wall. We use Dirichlet boundary conditions for the equilibrium displacement field on all boundaries to completely prescribe the equilibrium shape. The mesh and boundary conditions are shown in Figure 7.

The solution to the inverse problem is shown in Figure 8. The upper right-hand corner of the beam in the initial configuration is displaced by $-9.23\times 10^{-3}$ in horizontal direction and $63.97\times 10^{-3}$ in vertical direction. While the shape may appear to be a mirror image of the solution to the forward problem (Figure 6, bottom), the numbers prove that this is not the case.

5.2 Thermoelastic body under thermal stresses

Our principal test case involves stresses caused by thermal expansion or contraction. The idea is that we start with an inhomogeneously heated body (such as an injection molding part that has been partially cooled from the outside) that is cooled down to an equilibrium temperature. Our example geometry is shown in Figure 9.

We created an artificial temperature field as the initial condition of our simulation. To do so we solved the standard heat conduction equation, configured according to Table 5. Using Dirichlet boundary conditions, a hot body is slowly cooled down. Since we required an inhomogeneous temperature field, this simulation was stopped after five time steps. The result of this simulation is shown on the left-hand side in Figure 10.

Table 6 shows the material parameters used for the thermoelastic forward simulation. One should note that these values do not correspond to any existing material but are chosen for qualitative testing purposes.

Starting from the result of the heat conduction simulation, the actual thermoelastic forward simulation can be run. Figure 10 shows the change that the shape undergoes when the body cools down to a homogeneous temperature.

In addition, Figure 11 shows the residual stresses in the material, when starting from a stress-free state.

The same inhomogeneous temperature field can be used as basis for the inverse simulation. The left-hand side of Figure 12 shows the optimal shape of the body at ejection time, under the assumption of minimal changes in the temperature field, along with the adjusted temperature field.

Figure 13 highlights how our simulation, while ensuring a perfectly shaped result, does not prevent residual stresses in the body in any way.

6 Summary and Conclusion

In this paper, we have shown a new method for the solution of inverse design problems in thermoelasticity, including the formulation of the full, discretized system of equations. Using numerical examples, we have also demonstrated the usability of the method. The examples show that the method works and performs according to intuitive expectations.

One of the key difficulties we dealt with was the prescription of temperature and stress fields on an unknown shape. The proposed method of fitting these fields, when given for a specific shape, into adjusted shapes in a way that causes minimal changes is also demonstrated in the examples and follows the expectations.

The examples show how, in contrast to many shape optimization problems, the objective is achieved exactly. I.e., the desired shape is exactly obtained by the method. This is an optimal result, but one should not forget that it comes at the price of residual stresses in the body. However, since the stress field is calculated as a by-product of the inverse method, subsequent analyses can be performed on the stresses, e.g., if some material tolerances are known.

By choosing an inverse formulation rather than a shape optimization method, we have avoided many issues that come with these methods, such as the need to find a suitable objective function, obtain its derivatives, and choose a low-dimensional parameterization of the unknown shape. Overall, this method has the potential of yielding better results than shape optimization methods in less time. One should note that this is a direct consequence of choosing fully elastic material laws. Any other laws that would require a transient simulation cannot be handled by this method.

At this point, we are confident that this method is ready to be applied to realistic problems in various fields. As we have mentioned before, the method was originally motivated by the cavity shape determination problem in injection molding. Its development marks some considerable progress in this direction, but some work remains to be done. Specifically, efficient methods for the determination of the temperature and stress fields for provided cavity shapes are still under development. Challenges such as the precise modeling of polymer crystallization have to be faced to obtain optimum results.

Appendix A Derivations

A.1 Reformulation of Constitutive Law Based on St. Venant-Kirchhoff

We have defined a constitutive law for thermoelasticity based on the St. Venant-Kirchhoff material. This is given as a relationship between the second Piola-Kirchhoff tensor $\hat{\mathbf{S}}$ and the Green-St. Venant strain tensor $\mathbf{E}$:

[TABLE]

Our differential equations require constitutive laws to be formulated for the first Piola-Kirchhoff stress tensor $\mathbf{P}$ with respect to the reference configuration. The Green-St. Venant strain tensor $\mathbf{E}$ is defined as follows:

[TABLE]

The deformation gradient tensors $\tilde{\mathbf{F}}$ and $\bar{\mathbf{F}}$ can be calculated directly from the fields ${\mathbb{r}}$ and ${\mathbb{u}}$, respectively (see Figure 3). However, since we assume the initial deformation to be provided in the form of the left Cauchy-Green deformation tensor $\mathring{\mathbf{B}}$, the tensor $\mathring{\mathbf{F}}$ is not available, and therefore we cannot calculate $\mathbf{E}$. In order to make Equation (36) useful to us, we need to reformulate it in such a way that the stress tensor $\mathbf{P}$ can be calculated from the deformation tensors $\tilde{\mathbf{F}}$, $\bar{\mathbf{F}}$ and $\mathring{\mathbf{B}}$. At this point it should be noted that the determinant $\mathring{J}:=\mathrm{det}\,\mathring{\mathbf{F}}$ can be calculated from $\mathring{\mathbf{B}}$ as $\mathring{J}=\sqrt{\mathrm{det}\,\mathring{\mathbf{B}}}$.

We will make use of the following identities [27]:

[TABLE]

The first step is to find a relationship between the tensors $\mathbf{P}$ and $\hat{\mathbf{S}}$:

[TABLE]

The following identities, which can be derived easily, will be helpful in the reformulation:

[TABLE]

where we have used $\mathbf{B}=\bar{\mathbf{F}}\tilde{\mathbf{F}}^{-1}\mathring{\mathbf{B}}\tilde{\mathbf{F}}^{-T}\bar{\mathbf{F}}^{T}$.

Now, starting from Equation (LABEL:eq:relprefstot) and using the other relations, we can formulate the finished constitutive law:

[TABLE]

A.2 Reformulation of Constitutive Law Based on Neo-Hooke

The second constitutive law that we have mentioned is based on the Neo-Hooke material. The Cauchy stress according to this law can be related to the left Cauchy-Green deformation tensor $\mathbf{B}$ by

[TABLE]

We need to formulate $\mathbf{P}$ in terms of $J\boldsymbol{\sigma}$:

[TABLE]

The law can now be inserted into the equation:

[TABLE]

Declarations

Competing interests

The authors declare that they have no competing interests.

Funding

The presented investigations were carried out at RWTH Aachen University within the framework of the Collaborative Research Centre SFB1120-236616214 ”Bauteilpräzision durch Beherrschung von Schmelze und Erstarrung in Produktionsprozessen” and funded by the Deutsche Forschungsgemeinschaft e.V. (DFG, German Research Foundation). The sponsorship and support is gratefully acknowledged.