# Precision and Flexible Bending Control Strategy Based on Analytical Models and Data Models – Chinese Journal of Mechanical Engineering

Aug 12, 2022

### Flexible Hardening Model for Sheet Metal

To improve the adaptability of the flexible system, the material hardening model which have a wider description capability will be selected preferentially. From references, one type is constructed in the form of a power function, while the other is in the form of an exponential function. Among the power function type are the Hollomon and Ludwik models, whereas the exponential function models include the Voce, Voce+Voce, Voce++, and Hockett–Sherby models. The Ludwik type is an improved variant of the Hollomon model. The power function model is usually an unsaturated one, while the Voce type is saturated, and the other exponential function models are variants of the Voce model. So Ludwik model and Voce model are recommended firstly. So the first recommended elementary function is determined by referring to the Ludwik constitutive equation as presented thus:

sigma = left{ begin{aligned} & Evarepsilon , & quad varepsilon le varepsilon_{s} , hfill \ & sigma_{s} + Kleft( {varepsilon^{n} – varepsilon_{s}^{n} } right){, } & quad varepsilon ge varepsilon_{s} , hfill \ end{aligned} right.

(1)

where σ is true stress, ε is true strain, E is elastic modulus, εs is yield strain, σs is yield strength, K is strength coefficient, and n is hardening index. Referring to the Voce constitutive equation, the second recommended elementary function is determined as follows:

sigma = left{ begin{aligned} & Evarepsilon , & quad varepsilon le varepsilon_{s} , hfill \ & sigma_{s} + Aleft[ {exp ( – Bvarepsilon_{s} ) – exp ( – Bvarepsilon )} right]{, } & quad varepsilon ge varepsilon_{s} , hfill \ end{aligned} right.

(2)

where A and B are Voce model coefficients.

In order to test the hardening models, three representative curves are selected, as shown in Figure 1. High-strength steel DP980 have a higher yield stress, lower uniform elongation, and no evident yield turning point. The aluminum alloy 6061 sheets exhibits an apparent elastic-plastic transition. With the continuation of plastic deformation, the material hardening rate decreases slowly in a small range. For the 304 stainless steel material, there is a visible elastic–plastic transition curve, the hardening rate of the material, which remains larger than those of the other two, changes a little. Other test materials with similar behaviors are not listed there individually.

After fitting the stress and strain data of the material by the above two elementary functions, the fitting results are shown in the Figure 1(b)–(d). DP980 sheet metal is better fitted by Voce model than Ludwik model. The correlations of Voce model and Ludwik model with aluminum alloy 6061 sheet are similar. For the 304 stainless steel, it is same as aluminum alloy 6061 sheet.

When the elementary functions are used to extrapolate the high strain data, the Voce hardening rate changes little and tends to be saturated earlier, while the Ludwik model can continue to harden and be consistent with the cold working hardening behavior of the material. Therefore, when the material hardens slowly, the Ludwik and Voce models can be linearly combined to build a flexible hardening model:

sigma = left{ begin{aligned} & Evarepsilon , & quad varepsilon le varepsilon_{s} , hfill \ & sigma_{s} + alpha Kleft( {varepsilon^{n} – varepsilon_{s}^{n} } right) + (1 – alpha )Aleft[ {exp ( – Bvarepsilon_{s} ) – exp ( – Bvarepsilon )} right]{, } & quad varepsilon ge varepsilon_{s} . hfill \ end{aligned} right.

(3)

When α = 1 and α = 0, the combination function degenerates to elementary Eqs. (1) and (2), respectively. When α(0,1), the combination function can describe hardening behavior of diverse materials precisely. The parameter identification performance of different hardening models for various materials will be analyzed later.

The combination constitutive model has been determined above. It is necessary to give the corresponding analytical model for the hardening model. In the forming process, as the bending stroke of the punch increases, the bending angle of the sheet gradually increases. The mechanical model is shown in Figure 2.

According to the mechanical model established in Figure 2, the static equilibrium equation and the geometric relationship in y direction are obtained as

$$sum {F_{y} } = Ncos theta + mu Nsin theta – frac{P}{2} = 0,$$

(4)

$$N = frac{1}{cos theta + mu sin theta }frac{P}{2}.$$

(5)

Among them, μ is the friction coefficient, θ is the rotation angle at the die fillet, P is the bending force, rd is die radius, rp is punch radius, N is the direct reaction force at the fulcrum, M is bending moment, l is the distance from the origin point in x direction to the tangent point of the sheet and die, L is the distance from the origin point to the center of the round corner on the die. The different mass points on curved section are in different deformation stages. ρx is the curvature radius on the neutral layer at x, ρlim is the limited curvature radius of elastic. When ρx is larger than ρlim, particle which is on the section of the sheet is in the state of complete elastic deformation. When ρx is less than ρlim, particle which is on the section of the sheet is in the state of complete elastic deformation or plastic deformation. The two state calculation methods are as follows:

1. 1)

In the process of elastic deformation, the limited curvature radius of elastic ρlim is obtained by generalized Hooke’s law which is under the principal stress state:

$$rho_{lim } = frac{{E(1 – nu + nu^{2} )^{{{1 mathord{left/ {vphantom {1 2}} right. kern-nulldelimiterspace} 2}}} overline{y}}}{{sigma_{s} (1 – nu^{2} )}},$$

(6)

where ν is Poisson’s ratio, is the ordinate value under the elastic limit.

2. 2)

In the process of elastic-plastic deformation, the deformation zone can be simplified to plane strain state, which means εz=0. According to the assumptions, the stress in thickness direction during the air bending process is neglected, which means σy=0. And the deformation zone conforms to the small deformation assumption, which means the strain distribution on the section is linear. The bending moment of the sheet is calculated by the following formula:

$$M = M_{e} + M_{p} ,$$

(7)

where Me is elastic bending moment, Mp is plastic bending moment.

According to the first elementary function 1 and the second elementary function 2, the combined bending moment calculation model is established as

$$begin{gathered} M_{e} = 2bint_{0}^{{varepsilon_{xs} rho_{lim } }} {frac{E}{{1 – v^{2} }}frac{y}{{rho_{x} }}y{text{d}}y} , hfill \ M_{p} = 2bint_{{varepsilon_{xs} rho_{lim } }}^{t/2} {alpha left[ {frac{2}{sqrt 3 }left( {sigma_{s} – Kvarepsilon_{s}^{n} } right) + frac{2}{sqrt 3 }Kleft( {frac{2}{sqrt 3 }frac{y}{{rho_{lim } }}} right)^{n} } right]} y{text{d}}y + hfill \ quad quad int_{{varepsilon_{xs} rho_{lim } }}^{t/2} {(1 – alpha )left[ {frac{2}{sqrt 3 }left( {sigma_{s} + Aexp ( – Bvarepsilon_{s} )} right) + frac{2}{sqrt 3 }Aleft( { – exp ( – Bfrac{y}{{rho_{x} }})} right)} right]} y{text{d}}y, hfill \ end{gathered}$$

(8)

where y is the ordinate value of sheet cross section, b is the width of cross section , t is the thickness of the sheet metal, when α=1, the combined model degenerates to Eq. (1); when α=0, the combined model degenerates to Eq. (2); when α(0,1), the combined model can describe diverse relationships between bending moment and curvature.

As the bending stroke of the punch increases, the radius of curvature under the punch gradually decreases. When the radius of curvature under the punch is ρmin>rp+t/2, the area closing to the punch of the sheet did not appear, as shown in Figure 3(a); when the radius of curvature under the punch is ρmin=rp+t/2, as shown in Figure 3(b), the area closing to the punch of the sheet appeared.

According to the situations whether the die is attached or not, the curvature of each point on the sheet is calculated by the following formula:

$$rho_{x} = left{ begin{gathered} r_{p} + {t mathord{left/ {vphantom {t {2, , rho_{min } (M) le r_{p} + t/2,}}} right. kern-nulldelimiterspace} {2, , rho_{min } (M) le r_{p} + t/2,}} hfill \ rho (M), , rho_{min } (M) > r_{p} + t/2. hfill \ end{gathered} right.$$

(9)

Next, these control equations are used to establish the relationship between bending force and bending stroke during loading for each deformation state. Through a trial analysis, when the die span is larger, there is no monotonous relationship between the bending force and the stroke, in other words, when the stroke is increasing, the bending force which is corresponding to the next stroke increment may be larger or smaller than the current force. In order to simplify solution, a return back algorithm is used for calculation. At first, it is assumed that the bending force is larger than the current force, then the algorithm will find the corresponding stroke. If the stroke can not be found in the forward search direction, the search direction will be adjusted backward. Until the convergence condition is meet, the algorithm exits the search program. This calculation process is shown in Figure 4.

According to the reference, the elastic modulus is exponentially described as a function of the total strain, as shown in the following equation:

$$E_{u} = E_{0} – (E_{0} – E_{a} )left[ {1 – exp ( – xi varepsilon )} right],$$

(10)

where Eu is elastic modulus during unloading, E0 is initial elastic modulus, Ea is saturated modulus.

It can be concluded that under the bending moment, the relationship between the springback curvature and the bending moment is shown as follows:

$$M = 2int_{0}^{t/2} {E_{u} } frac{{y^{2} }}{{rho_{x}^{u} }}b{text{d}}y,$$

(11)

where ({rho_{x}^{u}}) is the radius of curvature at a point x after unloading the bending moment M. The amount of variation in curvature after springback can be further expressed as

$$Delta K = frac{M}{{E_{a} I + (E_{0} – E_{a} )I^{prime}}} = frac{M}{W},$$

(12)

$$frac{1}{{rho_{x}^{u} }} = frac{1}{{rho_{x} }} – Delta K,$$

(13)

where (I^{prime} = 2bleft( {frac{{t^{2} B}}{4A} – frac{tB}{{A^{2} }} + frac{2B}{{A^{3} }} – frac{2}{{A^{3} }}} right)), (A = – frac{xi }{{rho_{x} }}), (B = exp frac{At}{2}), W is bending stiffness.

### Sequential Identification Strategy for Nominal Parameters

The force and stroke sensors collect data from the bending machine synchronously, The original data are recorded as (hoi, Poi), i=1, 2,…, me. The curves of the data are shown in Figure 6. The interpolation algorithm can be used to resample the experimental data at equal intervals to solve the problem that the experimental data and the predicated data from the mechanical prediction model are unequal on the abscissa. The interpolated experimental data points are recorded as (hei, Pei), i=1, 2,…, mint. The sampling interval is Δh, and the calculation formula is presented below. The force stroke data points given by the analytical model are recorded as (hai, Pai), i=1, 2,…, mint, and mint is the number of experimental data set after resampling.

$$P_{j} = P_{i}^{e} + left( {jDelta h – h_{i}^{e} } right)frac{{P_{i + 1}^{e} – P_{i}^{e} }}{{h_{i + 1}^{e} – h_{i}^{e} }}{,}; , j = 1, 2,…, m_{{text{int}}} ,;i in left{ {1,;2,…,;m_{e} } right}left| {left( {h_{i}^{e} le jDelta h le h_{i + 1}^{e} } right),} right.$$

(14)

$$h_{j} = jDelta h, , j = 1,;2,…,;m_{{text{int}}} .$$

(15)

Based on the above-described model, these experimental data points can be divided into four segments according to different purposes. The data points before yield bending are used to optimize the elastic modulus, and the transition points between elastic and plastic bending are used to optimize the yield strength, while the data points after yield bending are used to optimize the hardening parameters, and the unloading data points are used to optimize the variable modulus parameters. The window method [32] is used to automatically identify the each segment boundary points of the experimental and predicated force–stroke curves. Finally, automatic curve segmentation is completed.

The boundary points for each segment on the curve can be stably extracted by moving the window. The thresholds of the abscissa and ordinate are set to Δhc and ΔFc, respectively. In order to keep the threshold universal, the curve data are normalized, the number of window points is determined by the following formula:

$$n_{w}^{i} = left{ {omega_{i} :omega_{i} in left{ {2,;4,…,;m_{{text{int}}} – i} right}left| {left( begin{gathered} sumlimits_{j = 1}^{{omega_{i} }} {abs(h_{i + j}^{norm} – h_{i + j – 1}^{norm} ) ge Delta h_{c} } vee hfill \ left( {F_{i + j}^{norm} – F_{i}^{norm} } right) ge Delta F_{c} ,;j = 2,;4,…,;omega_{i} hfill \ end{gathered} right)} right.} right},;i = 1,;2,…,$$

(16)

where (i) is the index of the starting point (hinorm, Finorm) of the window data. When the length of the horizontal coordinate or the height of the vertical coordinate of the window whose end point (hnormi+j, Fnormi+j) meets the conditions, niw+1 is the adaptive window size which i is the index starting point. the points between [i, i+niw/2] and [i+niw/2, i+niw] are fitted linearly respectively, ois is defined as the starting point on the straight line, oim is defined as the middle point, oif is defined as the end point, the rotation angle is calculated from crossing the two vectors, it is expressed as follows:

$$alpha_{i} = – frac{{overrightarrow {{o_{s}^{i} o_{m}^{i} }} times overrightarrow {{o_{m}^{i} o_{f}^{i} }} }}{{left| {overrightarrow {{o_{s}^{i} o_{m}^{i} }} times overrightarrow {{o_{m}^{i} o_{f}^{i} }} } right|}}arccos left( {frac{{overrightarrow {{o_{s}^{i} o_{m}^{i} }} cdot overrightarrow {{o_{m}^{i} o_{f}^{i} }} }}{{left| {overrightarrow {{o_{s}^{i} o_{m}^{i} }} } right|left| {overrightarrow {{o_{m}^{i} o_{f}^{i} }} } right|}}} right),;i = 1,;2,…$$

(17)

The objective functions of the elastic, yield, plastic, and elastic unloading segments are defined as follows. When the predicated and experimental data are sufficiently close, in other words, each of the objective functions takes the minimum value, the selected nominal parameters are reasonable. Now four sub-optimization problems are built. They are optimized one by one sequentially, this is a sequential identification strategy for nominal parameters.

$$G_{e} (X_{k}^{e} ) = frac{1}{2}sumlimits_{{i in Omega_{e} }}^{{}} {frac{{left( {P_{i}^{e} – P_{i}^{a} } right)^{2} }}{{i_{e} }}} ,;i = 1,;2,…,;i_{e} ,$$

(18)

$$G_{y} (X_{k}^{y} ) = frac{1}{2}sumlimits_{{i in Omega_{y} }}^{{}} {frac{{left( {P_{i}^{e} – P_{i}^{a} } right)^{2} }}{{i_{y} }}} ,;i = 1,;2,…,;i_{y} ,$$

(19)

$$G_{p} (X_{k}^{p} ) = frac{1}{2}sumlimits_{{i in Omega_{p} }}^{{}} {frac{{left( {P_{i}^{e} – P_{i}^{a} } right)^{2} }}{{i_{p} }}} ,;i = 1,;2,…,;i_{p} ,$$

(20)

$$G_{u} (X_{k}^{u} ) = frac{1}{2}sumlimits_{{i in Omega_{u} }}^{{}} {frac{{left( {P_{i}^{e} – P_{i}^{a} } right)^{2} }}{{i_{u} }}} ,;i = 1,;2,…,;i_{u} .$$

(21)

The nominal parameter ranges for each sub-optimization problem are set as

$$left{ begin{gathered} a_{k}^{e} le X_{k}^{e} le b_{k}^{e} , , k = 1,;2, cdots ,;l_{e} , hfill \ a_{k}^{y} le X_{k}^{y} le b_{k}^{y} , , k = 1,;2, cdots ,;l_{y} , hfill \ a_{k}^{p} le X_{k}^{p} le b_{k}^{p} , , k = 1,;2, cdots ,;l_{p} , hfill \ a_{k}^{u} le X_{k}^{u} le b_{k}^{u} , , k = 1,;2, cdots ,;l_{u} , hfill \ end{gathered} right.$$

(22)

where le, lv, lp, lu are the number of nominal parameters that need to be optimized, and variables Xek, Xyk, Xpk, Xuk are the degree of freedom. For each sub-optimization problem, a quadratic approximation function is used for each objective function. At the kth iteration, there is an approximate function shown as follows:

$$Q_{k} left( {{varvec{Y}}_{j} } right) = G_{s} left( {{varvec{Y}}_{j} } right), , j = 1,;2,;3, cdots ,;m,$$

(23)

where m[(ls+2), (ls+1)(ls+2)/2], s{e, y, p, u}. The optimization problems of above four objective functions are transformed into the minimizing problems of approximate functions. A optimization problem is shown as follows:

begin{aligned} & min , {Q_{k}} ({varvec{X}}_{k} + {varvec{d}}), \ & {text{s.t.}} {varvec{a}} le {varvec{X}}_{k} + {varvec{d}} le {varvec{b}}, | {varvec{d}} |_{2} le Delta_{k} , hfill \ & {varvec{X}}_{k} in { {varvec{Y}}_{j} :j = 1,;2, ldots ,m} left| min left{ G_{s} ( {varvec{Y}}_{j} ):j = 1,;2, ldots , m right} right., end{aligned}

(24)

where d is the vector step of each iteration and Δk is the confidence radius at the kth iteration.