Numerical optimization of a new robotized glass blowing process

Abstract

To answer an increasing need for glass product manufacturing in both small and medium series, the first glass-blower robot was recently developed. In the face of this new technology, which particularly interests crystal glass-makers, expertise remains the main decision-making element which intervenes in the choices of the design and implementation of these new processes. Finite element models of this new blowing process were developed. After the analysis of the process and of these stages, an initial sensibility study allowed us to find the essential parameters for the success of the operation. With the results of these sensitivity analyses, an optimizer was developed to adjust virtually the forming process of a linear cylindrical vase by determining the optimal forming parameters. A second optimization allowed us to determine the initial shape of the parison, an essential parameter in the successful forming of a convex cylindrical vase. Finally, the numerical tools were validated during trial campaigns carried out in crystal glass-makers.

Appendices

Appendix 1. Set of material properties for numerical modeling of rheological testing

The set of glass material properties used for the process modeling are indicated in the Table 4. The set of glass material properties is proposed by:

• Sharp, for specific heat [14],

• Primenko and Van-Zee, for conductivity [15, 16],

• Cesar De Sa, for the convection coefficients with the mould [3].

Table 4 Glass material properties

The viscosity and the convection coefficient with the ambient air were determined by two characterization tests developed by the “Modeling of Glass Forming and Tempering “ team of the L.A.M.I.H. of Valenciennes (Fig. 6). The mould material properties (wood) are recapitulated in the Table 5. The initial mould temperature was 150°C.

Fig. 6

Heat transfer coefficient between the glass and the mould as a function of the contact time proposed by Cesar De Sa [3]

Table 5 Material properties of the moulds

Appendix 2. Inverse method and optimization procedure

The identification method used is based on an iterative diagram in three stages. The aim is to determine N parameters p _i in order to reduce the difference between the objective solution ${\mathop f\limits^ \wedge }$ and the numerical solution f ^* ( p _i ). The error function E to be minimized is expressed by:

$$ E{\left( {p_{i} } \right)} = \frac{1} {2}{\sum\limits_{s = 1}^M {{\left[ {r_{s} {\left( {p_{i} } \right)}} \right]}^{2} } }, $$

(8)

with

$$ r_{s} {\left( {p_{i} } \right)} = f^{*}_{s} {\left( {p_{i} } \right)} - {\mathop {f_{s} }\limits^ \wedge },\;{\text{and }}i = 1,\;N, $$

(9)

where M is the number of comparison points.

Constraints can be added to limit the identified parameters to a field of validity by penalizing the error function. These constraints C _j ( p _i ) are expressed by

$$ C_{j} {\left( {p_{i} } \right)} \geqslant 0,\;\;\;j = 1,\;q{\text{,}} $$

(10)

where q represents the number of constraints.

The rewritten error function is then:

$$ E^{ * } {\left( {p_{i} } \right)} = E{\left( {p_{i} } \right)} + {\sum\limits_{j = 1}^q {\xi _{j} {\left( {p_{i} } \right)}} }, $$

(11)

with

$$ \xi _{j} {\left( {p_{i} } \right)} = \frac{{\omega _{j} {\left( {p_{i} } \right)}}} {{C_{j} {\left( {p_{i} } \right)}}}, $$

(12)

where ξ _j ( p _i ) is the penalization function suggested by Caroll [12], and ω _j ( p _i ) a positive scalar of penalization that Schnur [13] proposes to determine by

$$ \omega _{j} {\left( {p_{i} } \right)} = 0.0001 \cdot C_{j} {\left( {p_{i} } \right)} \cdot E{\left( {p_{i} } \right)}. $$

(13)

When the initial parameters are chosen the identification can begin. The resolution uses an algorithm developed by Levenberg–Marquardt [13]. At the k th iteration, the dp _i corrections added to the p _i parameters are expressed by:

$$ {\left[ {{\left( {{\mathbf{J}}^{k} } \right)}^{T} \cdot {\mathbf{J}}^{k} + \lambda ^{k} \cdot {\mathbf{I}} + H^{k} } \right]} \cdot {\text{d}}p^{k}_{i} = - {\left( {{\mathbf{J}}^{k} } \right)}^{T} \cdot r^{k} + G^{k} , $$

(14)

where λ is the Levenberg–Marquardt parameter. The Jacobian matrix J of the error function E ( p _i ), G and H the first and second derivatives, respectively, of the penalization function ξ _j ( p) are expressed by:

$$ {\mathbf{J}}_{{si}} = \frac{{\partial r_{s} }} {{\partial p_{i} }} = \frac{{\partial f^{ * }_{s} }} {{\partial p_{i} }},\;\;\;s = 1,\;M,\;\;{\text{and}}\;\;i = 1,\;N; $$

(15)

$$ G_{i} = {\sum\limits_{j = 1}^q {\frac{{\partial \xi _{j} }} {{\partial p_{i} }}} },\;\;\;i = 1,\;N{\text{;}} $$

(16)

$$ H_{{il}} = {\sum\limits_{j = 1}^q {\frac{{\partial ^{2} \xi _{j} }} {{\partial p_{i} \cdot \partial p_{l} }}} },\;\;\;i,\;l = 1,\;N{\text{;}} $$

(17)

where J _si (the terms of the Jacobian matrix) are approximated by finite differences. A coefficient of sensitivity S _p is added to the algorithm for the calculation of the partial derivative of each parameter:

$$ \frac{{\partial f^{ * } {\left( {p_{i} } \right)}}} {{\partial p_{i} }} = \frac{{f^{ * } {\left( {p_{i} + \Delta p_{i} } \right)} - f^{ * } {\left( {p_{i} } \right)}}} {{\Delta p_{i} }},\;\;\;{\text{with}}\;\Delta p_{i} = p_{i} S_{p} . $$

(18)

The Levenberg–Marquardt parameter λ determines simultaneously the direction and the size of the dp correction. Marquardt showed that if this parameter tends to the infinite, the direction will take the greatest slope and the step of correction will be smaller [17].

The determination of the solution ( p 1,..., pn) that minimizes the error function E ( pi) is made with the algorithm

• 1. Determine the initial values of the parameters to identify p _i ⁽⁰⁾, the sensitivity S _p and the Levenberg–Marquardt parameter λ ⁽⁰⁾.

• 2. Calculate the numerical solution of the problem f ^* ( p _i ⁽⁰⁾) and determine the error function E ⁽⁰⁾.

• 3. Determine the initial constraint function ξ _j ⁽⁰⁾ and the error function E *⁽⁰⁾.

• 4. For an iteration k ≥1:

  (a):

  Calculate the Jacobian matrix J ^{( k)} by solving an approximation by finite differences.

  (b):

  Calculate the first derivative H ^{( k)} and the second G ^{( k)} of the constraint function.

  (c):

  Determine the parameter corrections dp _i ^{( k)} by the resolution of Eq. (14) p _i ^{( k +1)}= p _i ^{( k)}+dp _i ^{( k)}.

  (d):

  Calculate the numerical solution of the problem f^* (p_i ^{(k +1)}), the constraint function ξ _j ^{( k +1)} and the error function E ^{*( k +1)}.

  (e):

  Check the convergence of the solution E ^{*( k +1)}< E ^{*( k)}. If the convergence is checked and if the error is lower than the stop criteria, the identification is finished. If the convergence is checked and if the error is higher than the stop criteria, the identification goes on (stage f). If the convergence is not checked, multiply the Levenberg–Marquardt parameter by 10 ( λ ^{( k)}=10 λ ^{( k)}) and begin again at stage (c).

  (f):

  Reduce the weight factors ω _j ( p _i )^{( k)} and λ ^{( k)}, increment k to 1, begin again at stage (a).