# Optimization of the thickness law for thin axi-symmetric shells

Jean-Marie Boisserie, Roland Glowinski

► To cite this version:

Jean-Marie Boisserie, Roland Glowinski. Optimization of the thickness law for thin axi-symmetric shells. [Research Report] IRIA-RR-263, IRIA. 1977, pp.43. hal-04716438

HAL Id: hal-04716438

https://inria.hal.science/hal-04716438v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Drw = 03639 RR 263

RIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

OPTIMIZATION OF
THE THICKNESS LAW
FOR THIN
AXI-SYMMETRIC SHELLS

Jean-Marie BOISSERIE
Roland GLOWINSKI

Rapport de Recherche N° 263

Novembre 1977
OPTIMIZATION OF THE THICKNESS LAW FOR THIN
AXI-SYMMETRIC SHELLS

Jean-Marie Boisserie, Roland Glowinski

(Dedicated to Professor J.H. Argyris)

Résumé :

On étudie dans cet article la résolution numérique de problèmes de conception optimale rela-
tifs à des coques minces de révolution. Les méthodes de résolution sont basées sur l'utilisation
d'éléments finis à une variable d'espace, du contrôle optimal et de la programmation non linéaire.
On insiste particulièrement sur le calcul du gradient de fonctions coût dont la définition fait inter-
venir une équation d'état. Les résultats d'expériences numériques sont présentées.

Abstract :

In this paper one studies the numerical solution of Optimum Structural Design problems for thin axi-symmetric elastic shells. The methods of solution are based on the use of Finite Elements in one space variable Optimal Control and Non linear Programming. One emphasizes particularly on the calculation of the gradient of cost functionals whose definition involves a state equation. Results from experiments are presented.
## 1. INTRODUCTION

The main goal of this paper is to study the numerical solution of a class of simple *Optimum Structural Design* problems. The structures under consideration are *thin axi-symmetric shells*. To achieve the solution of these structural problems we use a methodology, based on *Optimal Control*, combined with a *finite element approximation* in one space variable, of the *linear elasticity* equations describing the mechanical behaviour of the shells.

Results from numerical experiments are presented.

## 2. NOTATION

In the sequel we use the following notation (illustrated by Fig. 2.1) :

- $\Sigma$ is an axi-symmetric shell,

- $S$: meridian of the middle surface,

- $e$: thickness of the shell measured on the normal to the meridian of the middle surface,

- $\vec{n}$: unit vector of the outward normal to the middle surface,

- $P_i = P - \frac{\vec{n}}{2} e$, $P_e = P + \frac{\vec{n}}{2} e$, $|P_i \cdot P_e| = e$.
The normal component of the stresses being neglected, we respectively denote by $\sigma_s$ and $\sigma_\theta$ the meridional and circumferential stresses. More precisely we consider the vector

$$ \sigma = \{\sigma_1, \sigma_2, \sigma_3, \sigma_4\} $$

defined by

$$ (2.1) \qquad \begin{cases} \sigma_1 = \sigma_s(P_i), \quad \sigma_2 = \sigma_s(P_e), \\ \sigma_3 = \sigma_\theta(P_i), \quad \sigma_4 = \sigma_\theta(P_e). \end{cases} $$

The other notation are obvious from Figure 2.1.

<figure><img src="image_2.png" /><figcaption>Fig. 2.1.</figcaption></figure>

We also have to consider the displacement and strain vectors, respectively defined by,

$$ (2.2) \qquad u = \{u_1, u_2\} $$

where $u_1$ is the meridional displacement and $u_2$ the normal displacement, and
$$ (2.3) \qquad \varepsilon = \{\varepsilon_1, \varepsilon_2, \varepsilon_3, \varepsilon_4\} $$

where (see, e.g., ZIENKIEWICZ [1, Ch. 12])

$$ (2.4) \qquad \left\{ \begin{array}{l} \varepsilon_1 = \frac{du_1}{ds}, \quad \varepsilon_2 = \frac{u_1 \sin \alpha + u_2 \cos \alpha}{r}, \\ \\ \varepsilon_3 = - \frac{d^2 u_2}{ds^2}, \quad \varepsilon_4 = - \frac{\sin \alpha}{r} \frac{du_2}{ds}. \end{array} \right. $$

### 3. FORMULATION OF THE OPTIMIZATION PROBLEMS

We use in this section the terminology of Optimal Control (see, e.g., LIONS [2]).

#### 3.1. Cost functions

The cost function is chosen so that we have the stresses as little as possible; then an obvious choice for the cost function is

$$ (3.1) \qquad j_1 = \frac{1}{2} \int_S \left( \sum_{k=1}^{4} \sigma_k^2 \right) r(s) \, ds. $$

A more sophisticated possibility is to take a cost function directly related to the second invariant of the stress tensor. Such a choice is more satisfactory from a mechanical point of view. This second invariant $I_2(\bar{\sigma})$ being defined by

$$ (3.2) \qquad \frac{3}{2} I_2(\bar{\sigma}) = \sigma_s^2 + \sigma_\theta^2 - \sigma_s \sigma_\theta $$

we define a second cost function (actually a family of cost functions) by

$$ (3.3) \qquad j_2 = \frac{1}{2} \int_S \left( \beta_i(s) I_2(\bar{\sigma}_i) + \beta_e(s) I_2(\bar{\sigma}_e) \right) r(s) \, ds. $$

In (3.3), $\beta_i$ and $\beta_e$ are weighting functions, and $\bar{\sigma}_i$, $\bar{\sigma}_e$ are respectively related to the stresses on the interior (inner) and exterior (outer) surfaces of the shell.
### 3.2. The state equations

The state equations are the equations describing the elastic behaviour of the shell under consideration. Using the classical Love-Kirchoff-Koiter shell theory it is convenient to introduce the following energy functional of the shell (in which vectors are identified to the corresponding column matrices)

$$ (3.4) \quad \frac{1}{2} \int_S \varepsilon^t(v)[D]\varepsilon(v)r(s)ds - \int_S F.vr(s)ds. $$

By differentiation of the functional (3.4) we obtain the variational formulation of the elasticity equations

$$ (3.5) \quad \left\{ \begin{array}{l} \displaystyle \int_S \varepsilon^t(u)[D]\varepsilon(v)r(s)ds = \int_S F.vr(s)ds \quad \forall v \in V_o, \\ u \in V. \end{array} \right. $$

In (3.4), (3.5)

- $V_o$ is the space of the test functions, i.e. the space of the displacements, compatible with the boundary conditions that the shell has to satisfy.

- $V$ is the space in which the solution is sought; the boundary conditions may occur in the definition of $V$.

- $F$ is the vector of the external forces acting on the shell.

In (3.4), (3.5) we have $[D]$ defined by

$$ (3.6) \quad [D] = \frac{Ee}{1-v^2} [d] $$

where $E$ is the Young modulus and $v$ the Poisson ratio and where

$$ (3.7) \quad [d] = \begin{bmatrix} 1 & v & 0 & 0 \\ v & 1 & 0 & 0 \\ 0 & 0 & \frac{e^2}{12} & v \frac{e^2}{12} \\ 0 & 0 & v \frac{e^2}{12} & \frac{e^2}{12} \end{bmatrix}. $$
If F is given and e is known, the displacement u can be obtained
by solving (3.5). Then σ is derived from u by (2.4) and

$$
(3.8) \quad \sigma = [MD][D]\epsilon,
$$

where

$$
(3.9) \quad [MD] = \begin{bmatrix} \frac{1}{e} & 0 & \frac{6}{e^2} & 0 \\ \frac{1}{e} & 0 & -\frac{6}{e^2} & 0 \\ 0 & \frac{1}{e} & 0 & \frac{6}{e^2} \\ 0 & \frac{1}{e} & 0 & -\frac{6}{e^2} \end{bmatrix} .
$$

It is of fundamental importance to observe that F being given,
$j_1$ and $j_2$ (cf. (3.1), (3.3)) are, in fact, functions of e, only,
through (3.5)-(3.9). Using the above notation we have

$$
(3.10) \quad j_1(e) = J_1(u(e),e)
$$

where

$$
(3.11) \quad J_1(v,e) = \frac{1}{2} \int_S \varepsilon^{\mathrm{t}}(v) G_1(e) \varepsilon(v) r(s) \mathrm{d}s
$$

with

$$
(3.12) \quad G_1(e) = [D(e)]^t [MD(e)]^t [MD(e)][D(e)] .
$$

Similarly we have

$$
(3.13) \quad j_2(e) = J_2(u(e),e)
$$

where

$$
(3.14) \quad J_2(v,e) = \frac{1}{2} \int_S \epsilon^{\mathrm{t}}(v) G_2(e) \epsilon(v) r(s) \mathrm{d}s
$$

with

$$
(3.15) \quad G_2(e) = [D(e)]^t [MD(e)]^t [C(s)] [MD(e)] [D(e)] .
$$

In (3.15), we have
$$ (3.16) \quad [C(s)] = \begin{bmatrix} \beta_i(s) & 0 & -\frac{1}{2}\beta_i(s) & 0 \\ 0 & \beta_e(s) & 0 & -\frac{1}{2}\beta_e(s) \\ -\frac{1}{2}\beta_i(s) & 0 & \beta_i(s) & 0 \\ 0 & -\frac{1}{2}\beta_e(s) & 0 & \beta_e(s) \end{bmatrix} $$

where $\beta_i$ and $\beta_e$ are the weighting functions introduced in (3.3).

### 3.3. Set of feasible thickness

In the following we shall require that the thickness e belong to $E_f$ where $E_f$ is the convex set defined by

$$ (3.17) \quad E_f = \{e \mid 0 < e_m \le e \le e_M, \left|\frac{de}{ds}\right| \le C, 2\pi \int_S e(s)r(s)ds = W\}. $$

The relation $2\pi \int_S e(s)r(s)ds = W$ implies that the weight of the structure is a priori given.

It is clear that many other possibilities exist for the choice of $E_f$.

### 3.4. Formulation of the Optimization problems

Using the above notation, the minimization problem that we consider is

$$ (3.18) \quad \begin{array}{ll} \text{Min} & j_i(e), \\ \text{e} \in E_f \end{array} $$

where for i=1 (resp. i=2), $j_i$ is defined by (3.1) (resp. (3.3)).

In (3.1) (resp. (3.3)), $j_1$ (resp. $j_2$) is a function of e via (3.5)-(3.12) (resp. (3.5)-(3.9), (3.13)-(3.16)).

### 4. GENERALITIES ABOUT THE NUMERICAL SOLUTION OF (3.18)

Concerning the numerical solution of (3.18) we have to consider the two following fundamental steps :
### 4.1. Approximation

We approximate the continuous problem (3.18) by a finite dimensional problem. To do this we approximate $V, V_o, E_f$ by respectively $V_h, V_{oh}, E_f^h$, then the state equation (3.5) by a linear system in finite dimension. Then the functional $J_i$ (defined by (3.11) if $i=1$ and by (3.14) if $i=2$) is approximated by $J_{ih}$ which induces an approximation $j_{ih}$ of $j_i$. Using the above approximations we define the approximate optimization problems by

$$ (4.1) \quad \begin{array}{ll} \displaystyle \ Min_{e_h \in E_f^h} & j_{ih}(e_h), \\ \text{for } i=1,2. \end{array} $$

### 4.2. Iterative solution of the approximate problem

The approximate problems (4.1) are Non Linear Programming problems ; therefore, their numerical solution will require the use of an iterative method. We may use a method of gradient type, which has the following general description :

$$ (4.2) \quad e_h^o \in E_f^h \text{ is given} $$

then for $n \ge 0$,

$$ (4.3) \quad \left\{ \begin{array}{l} \text{assuming that } e_h^n \text{ is known, we compute } \nabla j_h(e_h^n), \text{ and then} \\ \text{by some process to be defined, we obtain } e_h^{n+1} \text{ from } e_h^n \text{ and} \\ \qquad \nabla j_h(e_h^n), \end{array} \right. $$

where $\nabla j_h$ denotes the gradient of $j_h$.

In (4.3) we have used $j_h$ instead of $j_{ih}$ or $j_{2h}$, since the methods to be described later can be applied to other cost functionals than $j_1$ or $j_2$.
The two above steps in the numerical solution of (4.1) will be
detailed in the following sections.

5. - <u>APPROXIMATION</u>

5.1. <u>Generalities</u>

The displacement u = {u₁,u₂} and the thickness e are approximated
by finite elements (in one space variable). More precisely we
introduce a subdivision {P₀,P₁...Pₙ} on the meridian S of the shell.
We denote by sⱼ the curvilinear abscissa of Pⱼ and we suppose that
sⱼ < sⱼ₊₁ ∀j ≥ 0; P₀ and Pₙ coincide with the extremities of S.
As usual we define h by

$$
(5.1) \quad h = \max_j |s_{j+1} - s_j| .
$$

Then we approximate u_l by u_{1h}, where u_{1h} is a function of the curvi-
linear abscissa s, <u>globally continuous</u> over S and <u>piecewise linear</u>
over [s_j, s_{j+1}] ∀j ≥ 0 (see Figure 5.1).

<figure><img src="image_3.png" /><figcaption>Fig. 5.1</figcaption></figure>
Hence the approximation of $u_1$ is of *Lagrange type*, the natural *degrees of freedom* being the values taken by $u_{1h}$ at the $P_j$, i.e. $u_{1h}(s_j)$, $j=0,...,N$.

To approximate $u_2$ we use a function $u_{2h}$, *globally continuous* over $S$ just as its first *derivative*; furthermore we take $u_{2h}$ *piecewise cubic* over $[s_j,s_{j+1}]$ $\forall j \ge 0$ (see Fig. 5.2).

<figure><img src="image_4.png" /><figcaption>Fig. 5.2.</figcaption></figure>

Hence the approximation of $u_2$ is of *Hermite type* with *two degrees of freedom* by node $P_j$ which are $u_{2h}(s_j)$ and $\frac{du_{2h}}{ds}(s_j)$.

Finally we approximate $e$ by $e_h$, *piecewise constant* over $[s_j,s_{j+1}]$ $\forall j \ge 0$ (see Figure 5.3).
<figure><img src="image_5.png" /><figcaption>Fig. 5.3</figcaption></figure>

**Remark 5.1** : We could also have used an approximation $e_h$ of $e$, globally continuous and piecewise linear, i.e. of the same type than for $u_1$.

## 5.2. Approximation of the state equation

Once the displacement $u$ and the thickness $e$ have been approximated as indicated in Sec. 5.1, one has to approximate the state equation (3.5). It will be done by using a variational formulation in finite dimension.

Let $W_{1h}$ (resp. $W_{2h}$) be the space of the functions of $s$ globally continuous over $S$ (resp. globally continuous just as their first derivative) and piecewise linear (resp. piecewise cubic) over $[s_j,s_{j+1}] \forall j \ge 0$. Therefore from Sec. 5.1. we have $u_h \in W_{1h} \times W_{2h}$.
Actually since we have to take into account the boundary conditions that $u$ has to satisfy, $u_h$ belongs to a space $V_h$ such that $V_h \subset W_{1h} \times W_{2h}$. Similarly the space $V_o$ of the continuous test-functions is approximated by $V_{oh} \subset W_{1h} \times W_{2h}$ in which one takes into account the fact that the approximate test-functions $v_h$ have to be compatible with the boundary conditions satisfied by $u_h$.

It is then natural to approximate the state equation (3.5) by

$$ (5.2) \qquad \left\{ \begin{array}{l} \displaystyle \int_S \varepsilon^t(u_h)[D(e_h)]\varepsilon(v_h)r(s)ds = \int_S F.v_h r(s)ds \quad \forall v_h \in V_{oh}, \\ \\ u_h \in V_h. \end{array} \right. $$

Actually the direct use of (5.2) leads to fairly heavy calculations in the derivation of the *linear system equivalent* to (5.2). To avoid, or at least to simplify these calculations we have used systematically *numerical integration* to compute the integral occurring in (5.2).

Finally $u_h$ is the solution of a *finite dimensional linear variational problem*, of the following type :

$$ (5.3) \qquad \left\{ \begin{array}{l} a_h(e_h; u_h, v_h) = L_h(v_h) \quad \forall v_h \in V_{oh}, \\ \\ u_h \in V_h. \end{array} \right. $$

In (5.3) :

• $a_h(e_h; \cdot, \cdot)$ is, if $e_h \in E_f^h$ is *fixed*, a *bilinear form*, *symmetric and positive definite* over $V_{oh} \times V_{oh}$,

• $L_h(v_h)$ is an approximation of $\int_S F \cdot v r(s) ds$.

**Remark 5.2** : $a_h(e_h; \cdot, \cdot)$ is positive definite over $V_{oh} \times V_{oh}$ if $e_h$ restricted to $]s_j, s_{j+1}]$ is strictly positive. ■
The method to obtain the linear system equivalent to (5.3) is
classical but owing to the practical importance of this step
we shall consider it with some details : using the above notation
let us define $\gamma_h$ by

$$
(5.4) \quad \gamma_h = W_{1h} \times W_{2h}.
$$

Let now $V_{\gamma h}$ be a complementary subspace of $V_{oh}$ in $\gamma_h$, then

$$
(5.5) \quad \mathcal{V}_h = \mathcal{V}_{oh} \oplus \mathcal{V}_h.
$$

Let $\beta_{oh}$ and $\beta_{\gamma h}$ be respectively basis of $V_{oh}$ and $V_{\gamma h}$, with

$$
(5.6) \quad \beta_{oh} = \{w_j\}_{j=1}^{N_{oh}}, \quad N_{oh} = \dim (V_{oh}),
$$

$$
(5.7) \quad \beta_{\gamma h} = \{\bar{\omega}_k\}_{k=1}^{N_{\gamma h}}, \quad N_{\gamma h} = \dim (V_{\gamma h}).
$$

It follows then from (5.5)-(5.7) that $\forall v_h \in \gamma_h$ we have the
unique decomposition

$$
(5.8) \quad v_h = \sum_{j=1}^{N_{oh}} \mu_j w_j + \sum_{k=1}^{N_{\gamma h}} \beta_k \bar{w}_k .
$$

in (5.8) the $\beta_k$ are determined by the <u>boundary conditions</u>
satisfied by $v_h$.

Let $u_h$ be the approximate solution (i.e. the solution of (5.3));
from (5.8) we have for $u_h$ the following expansion

$$
(5.9) \quad u_h = \sum_{j=1}^{N_{oh}} \lambda_j w_j + \sum_{k=1}^{N_{\gamma h}} \alpha_k \bar{w}_k.
$$

In (5.9), the $\alpha_k$ are known from the <u>boundary conditions</u> and the
$\lambda_j$ are the true unknowns of the problem.

Using the bilinearity of $a_h(e_h; \cdot, \cdot)$ it is clear that (5.3) is equi-
valent to the following linear system in $\lambda_j$, $j=1, \dots, N_{oh}$:

$$
(5.10) \quad \left\{
\begin{array}{l}
\displaystyle \sum_{j=1}^{N_{\text{oh}}} a_h(e_h; w_j, w_i) \lambda_j = L_h(w_i) - \sum_{k=1}^{N_{\gamma h}} a_h(e_h; \bar{w}_k, w_i) \alpha_k, \\
\qquad 1 \le i \le N_{\text{oh}}.
\end{array}
\right.
$$
From the properties of $a_h(e_h; \cdot, \cdot)$ the matrix

$$ (5.11) \quad A_h(e_h) = (a_h(e_h; w_j, w_i))_{1 \le i,j \le N_{oh}} $$

is symmetric and positive definite (if $e_h$ is conveniently chosen, see Remark 5.2), hence the linear system (5.10) has a unique solution and it may be solved by a Method of Cholesky.

**Remark 5.3** : A basis of $\nu_h$, implying the sparsity of $A_h(e_h)$ is easily found from the canonical basis of $W_{1h}$ and $W_{2h}$. These basis may be found in many text books on Splines and Finite Elements (see, e.g., ODEN-REDDY [3, Chap. 6, Sec. 6.3], STRANG-FIX [4, Ch. 1, Sec. 1.7]). For the sake of completeness we have shown on Fig. 5.4, 5.5 the graphs of typical basis functions in, respectively, $W_{1h}$ (piecewise linear) and $W_{2h}$ (piecewise cubic) (the modifications for the basis functions related to the boundary nodes are obvious).

<figure><img src="image_6.png" /><figcaption>Fig. 5.4.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 5.5.</figcaption></figure>

### 5.3. Approximation of the cost functions

To approximate the cost functions $j_i$ (see Sec. 3.1, 3.2) we use similar methods than those we have considered in Sec. 5.2 to approximate the state equation. We have in particular used numerical integration to derive $j_{ih}$ from $J_i(v,e)$ (see (3.10)-(3.12) and (3.13)-(3.16)).

### 5.4. Approximation of $E_f$

We recall (see Sec. 3.3) that the set of feasible thickness $E_f$ is defined by
$$ E_f = \{e | 0 < e_m \le e \le e_M, \left| \frac{de}{ds} \right| \le C, 2\pi \int_S e(s)r(s)ds = W \}. $$

It is easy to approximate $E_f$ by a convex set of piecewise constant functions over $S$, relatively to the subdivision $\{s_0, s_1, \dots, s_N\}$.
Let us denote by $e_{j+1/2}$, $j=0, \dots, N-1$, the value of $e_h$ over $]s_j, s_{j+1}]$. Then we approximate $E_f$ by the closed convex set $E_f^h$ defined by :

$$ e_h \in E_f^h \text{ iff} $$

$$ (5.12) \qquad e_m \le e_{j+1/2} \le e_M \quad \forall j=0,1,\dots,N-1, $$

$$ (5.13) \qquad \frac{\left| e_{j+1/2} - e_{j-1/2} \right|}{\left( \frac{s_{j+1} - s_{j-1}}{2} \right)} \le C \quad \forall j=1, \dots, N=1, $$

$$ (5.14) \qquad 2\pi \sum_{j=0}^{N-1} (s_{j+1} - s_j) r_{j+1/2} \quad e_{j+1/2} = W $$

where, in (5.14), $r_{j+1/2}$ may be defined either by

$$ r_{j+1/2} = r\left(\frac{s_j + s_{j+1}}{2}\right) $$

or by

$$ r_{j+1/2} = \frac{1}{2} (r(s_j) + r(s_{j+1})) \quad (\text{trapezoidal rule}). $$

We observe that from (5.12)-(5.14), $E_f^h$ is defined by linear inequality or equality constraints only; this is important from a computational point of view as it will be seen in the following sections.
### 5.5. Approximation of the Optimization problems

Using the results and notation of the above sections it is clear that the structural optimization problems (3.18) are approximated, $\forall i=1,2$, by

$$ (5.15) \quad \Min_{e_h \in E_f^h} j_{ih}(e_h). $$

Since $E_f^h$ is a closed, bounded subset of a finite dimensional vector space it is compact. Moreover, $\forall i=1,2$, $j_{ih}$ is a continuous function of $e_h$ over $E_f^h$. These properties imply that (5.15) has at least a solution. The problem of uniqueness is much more difficult, and in fact it is still open.

**Remark 5.4** : Using more sophisticated compactness arguments, based on Ascoli-Arzela Theorem (see, e.g. YOSIDA [5]) it can be proved that the continuous problems (3.18) have always a solution.

## 6. - ITERATIVE SOLUTION OF THE APPROXIMATE PROBLEMS

### 6.1. Orientation

The considerations to follow are fairly classical in Numerical Optimal Control and in Non Linear Programming. In Sec. 6.2. assuming that the cost function is differentiable, we give necessary conditions of optimality involving the differential of the cost function. In Sec. 6.3 we describe an algorithm of Frank-Wolfe type which makes use of the differential of the functional to minimize. Finally in Sec. 6.4, we give a method for computing the differential of a cost function whose definition involves a state equation.
## 6.2. Generalities. Necessary conditions of optimality.

We observe that the optimization problems considered in the previous sections have the following formulation

$$ (6.1) \qquad \begin{array}{ll} \underset{e \in E_f}{\text{Min}} & j(e) \\ & e \in E_f \end{array} $$

(the subscripts i and/or h have been omitted for the sake of simplicity).

In (6.1), $E_f$ is a closed convex subset of an Hilbert space $\mathcal{H}$ (possibly finite dimensional) whose scalar product is denoted by $(\cdot, \cdot)$. Moreover we assume that $j(\cdot)$ is differentiable (it is the case for the specific problems we have considered in the above sections).

Let $\bar{e}$ be a solution of (6.1), then

$$ (6.2) \qquad \left\{ \begin{array}{l} \bar{e} \in E_f \\ j(e) \ge j(\bar{e}) \quad \forall e \in E_f. \end{array} \right. $$

Then from the convexity of $E_f$ and from (6.2) we have

$$ (6.3) \qquad \left\{ \begin{array}{l} j(\bar{e}+t(e-\bar{e})) = j(te+(1-t)\bar{e}) \ge j(\bar{e}) \\ \forall t \in ]0,1[, \quad \forall e \in E_f. \end{array} \right. $$

It follows from (6.3) that

$$ (6.4) \qquad \frac{j(\bar{e}+t(e-\bar{e}))-j(\bar{e})}{t} \ge 0 \quad \forall t \in ]0,1[, \quad \forall e \in E_f. $$

Therefore at the limit in (6.4) when $t \to 0$, and since $j(\cdot)$ is differentiable, we have
$$ (6.5) \quad (\nabla j(\bar{e}), e-\bar{e}) \ge 0 \quad \forall e \in E_f; $$

in (6.5), $\nabla j(e)$ denotes the differential (one says also the gradient) of $j(\cdot)$ at $e$. If $\delta = \mathbf{R}^N$ and if we use as scalar product

$$ (6.6) \quad (x,y) = \sum_{i=1}^{N} x_i y_i \quad \forall x,y \in \mathbf{R}^N, $$

then

$$ (6.7) \quad \nabla j(e) = \{ \frac{\partial j}{\partial e_1}(e), \dots, \frac{\partial j}{\partial e_k}(e), \dots, \frac{\partial j}{\partial e_N}(e) \}, $$

i.e. we have recovered classical concepts.

The relations

$$ (6.8) \quad \left\{ \begin{array}{l} (\nabla j(\bar{e}), e-\bar{e}) \ge 0 \quad \forall e \in E_f, \\ \bar{e} \in E_f \end{array} \right. $$

which have to be satisfied by any solution of (6.1), are necessary conditions of optimality for the solutions of (6.1). These conditions are not sufficient in general (if $j(\cdot)$ is convex, conditions (6.8) are necessary and sufficient, i.e. (6.2) and (6.8) are equivalent).

**6.3. Solution of the approximated problems by a method of Frank-Wolfe**

We suppose in this section that (6.1) is an optimization problem in finite dimension and that $E_f$ is defined by linear equality or inequality relations; it follows from Sec. 5 that the approximate structural optimization problems of the previous sections obey these conditions.
To solve the above problems we have used an *iterative* method of
*Frank-Wolfe type* (see, e.g., CEA [6]) in which we take advantage
of the fact that $E_f$ is defined by *linear relations*.

The algorithm is as follows

(6.9) $e^o \in E_f$ is given,

then for $n \ge 0$, assuming that $e^n$ is known one obtains $e^{n+1}$ by

(6.10) One uses $e^n$ to compute $\nabla j(e^n)$.

Then one computes $\bar{e}^{n+1}$ by solving the linear program

$$
(6.11) \quad \left\{
\begin{array}{l}
(\nabla j(e^n), \bar{e}^{n+1} - e^n) \le (\nabla j(e^n), e - e^n) \quad \forall e \in E_f \cap K_d^n, \\
\bar{e}^{n+1} \in E_f \cap K_d^n,
\end{array}
\right.
$$

where

$$
K_d^n = \{e \mid |e - e^n| \le d\} ;
$$

in practice, instead of a constant d, one uses a decreasing sequence
${d_n}$ to define $K_d^n$.

Finally one obtains $e^{n+1}$ by solving the *one-dimensional* minimization problem

$$
(6.12) \quad \left\{
\begin{array}{l}
j(e^{n+1}) \le j(e) \quad \forall e \in [e^n, \bar{e}^{n+1}], \\
e^{n+1} \in [e^n, \bar{e}^{n+1}],
\end{array}
\right.
$$

for which very efficient solvers exist (methods of *Dichotomy*, *Fibonacci*,
etc..., see WILDE-BEIGHTLER [7] for a very complete description of such
methods).
To conclude this subsection we observe that (6.11) can be solved
by linear programming methods and that the use of algorithm (6.9)-
(6.12) requires the knowledge of $\nabla j$ at the points $e^n$ for $n \ge 0$. The
computations of $\nabla j$ is considered in the following Sec. 6.4.

6.4. Computation of $\nabla j$.

6.4.1. Generalities. Reduction to a canonical formulation.

The iterative method described in the above Sec. 6.3, requires at
each iteration the knowledge of $\nabla j(e^n)$. Since in the structural
optimization problems considered in this paper the cost functions
are defined through the solution of a state equation, obtaining
the gradient $\nabla j(e)$ of $j(\cdot)$ at $e$ is not completely trivial operation.
In fact, in our opinion, it is the most crucial point in the numerical
solution of most optimum design problems. For the sake of generality
(in view of other possible applications) we don't particularize the
cost functions and the convex $E_f$ under consideration.

Like in Sec. 6.2. we consider the following class of minimization
problems

$$
(6.13) \quad \begin{array}{ll} \text{Min} & j(e) \\ & e \in E_f \end{array}
$$

where $E_f$ is a closed, convex, non empty subset of an Hilbert space
$\varepsilon$ (possibly finite dimensional) whose scalar product is denoted by
$(\cdot, \cdot)_\varepsilon$.

We assume furthermore that $j(e)$ is actually defined by

$$
(6.14) \quad j(e) = J(u(e), e).
$$
In (6.14):

J(•, •) is a functional from $\mathcal{E} \times H$ to $\mathbb{R}$, where $H$ is also an Hilbert space (possibly finite dimensional) whose scalar product is denoted by $(\cdot, \cdot)_H$.

u(e) is the solution of the linear variational problem

$$ (6.15) \qquad \left\{ \begin{array}{l} (a(e;u(e),v) = L(v) \quad \forall v \in H_o, \\ u(e) \in H; \end{array} \right. $$

in (6.15), $H_o$ is a closed subspace of $H$ and $L(\cdot)$ is a linear functional, continuous from $H$ to $\mathbb{R}$. We suppose moreover (in order to have a well-posed problem (6.15)) that for every given, the functional $a(e;\cdot, \cdot)$ is bilinear, continuous, possibly not symmetric, and is strongly elliptic over $H_o \times H_o$. $\forall e \in E_f$ (this last property means that $\forall e \in E_f$, there exists $\alpha(e) > 0$ such that

$$ (6.16) \qquad a(e;v,v) \geq \alpha(e) \|v\|_H^2 \quad \forall v \in H_o). $$

Finally we assume that $u(e)$ can be written

$$ (6.17) \qquad u(e) = u(e) - u_o + u_o $$

where

$$ (6.18) \qquad u_o \text{ is known (and usually independent of } e). $$

$$ (6.19) \qquad u(e) - u_o \in H_o; $$

in our applications $u_o$ is known from the boundary conditions. Since $u_o$ is known, $u(e)-u_o$ may be viewed as the true unknown and we observe that from (6.15), (6.17), (6.19), $u(e)-u_o$ is the solution of
$$ (6.20) \quad \left\{ \begin{array}{l} a(e;u(e)-u_o,v) = L(v)-a(e;u_o,v) \quad \forall v \in H_o \\ u(e)-u_o \in H_o. \end{array} \right. $$

From the properties of $a(e;\cdot,\cdot)$ and $L(\cdot)$ it can be proved that $\forall e \in E_f$, the linear variational problem (6.20) has a unique solution (this existence and uniqueness result follows from the Lax-Milgram Theorem, proved in, e.g., YOSIDA, loc. cit. ; in this proof the property of strong ellipticity (6.16) plays a fundamental role).

**Remark 6.1** : From a practical point of view, in finite element codes one never has to know $u_o$, since if the boundary conditions are correctly taken into account, the splitting (6.17) is implicitly and automatically done. ■

The formulation (6.20) suggests what we call a *canonical formulation* for our optimization problems ; more precisely we consider the following class of problems, defined by (6.13), (6.14) and the *state equation*

$$ (6.21) \quad \left\{ \begin{array}{l} a(e;u(e),v) = L(e;v) \quad \forall v \in H_o \\ u(e) \in H_o \end{array} \right. $$

where $H_o$, $a(e;\cdot,\cdot)$ are as before and $L(e;\cdot)$ is, for $e$ given, linear and continuous from $H_o$ to $\mathbb{R}$.

It is clear that modulo a translation by $u_o$ the original problem (6.13)-(6.15) can be reduced to our canonical formulation (6.13), (6.14), (6.21).
### 6.4.2. Calculation of $\nabla j(e)$ in the canonical case. Applications.
We suppose that the above functionals $J(\cdot, \cdot)$, $a(\cdot; \cdot, \cdot)$, $L(\cdot; \cdot)$ are differentiable in all their arguments. We observe then that the differential $\nabla j(e)$ of $j(\cdot)$ at $e$ verifies

$$ (6.22) \quad \lim_{\substack{t \to 0 \\ t \neq 0}} \frac{j(e+t\delta e) - j(e)}{t} = (\nabla j(e), \delta e)_{\xi} (= \delta j(e)). $$

In (6.22) we don't have to assume that $\delta e$ is infinitesimal.
By differentiation we have

$$ (6.23) \quad \delta j(e) = \delta J(u(e),e) = \left(\frac{\partial J}{\partial u}(u(e),e), \delta u(e)\right)_H + \left(\frac{\partial J}{\partial e}(u(e),e), \delta e\right)_\xi $$

In (6.23), $\frac{\partial J}{\partial u}$ and $\frac{\partial J}{\partial e}$ are the partial gradients of $J$ and $\delta u(e)$ is the solution of

$$ (6.24) \quad \left\{ \begin{array}{l} a(e; \delta u(e), v) = \left( \frac{\partial L}{\partial e} (e; v), \delta e \right)_{\xi} - \left( \frac{\partial a}{\partial e} (e; u(e), v), \delta e \right)_{\xi} \quad \forall v \in H_o, \\ \\ \delta u(e) \in H_o; \end{array} \right. $$

(6.24) is obtained by differentiation of (6.21) and has a unique solution if the above conditions on $a(\cdot; \cdot, \cdot)$ and $L(\cdot; \cdot)$ are satisfied. It is of fundamental importance to observe that $(\frac{\partial a}{\partial e} (e; u(e), \cdot), \cdot)$ and $(\frac{\partial L}{\partial e} (e; \cdot), \cdot)$ are bilinear over $H_o \times \xi$.

The main drawback of (6.23), (6.24) is that it contains $\delta u(e)$ which is a linear function of $\delta e$ through (6.24). In order to eliminate $\delta u(e)$ we use the classical method of the adjoint equation; more precisely we introduce the following linear variational problem

$$ (6.25) \quad \left\{ \begin{array}{l} a(e;v,p(e)) = \left(\frac{\partial J}{\partial u}(u(e),e),v\right)_H \quad \forall v \in H_o, \\ \\ p(e) \in H_o. \end{array} \right. $$
From the above hypotheses on $a(e;\cdot,\cdot)$ and $J(\cdot,\cdot)$, (6.25) has a unique solution.

**Remark 6.2** : We have to observe that if $a(e;\cdot,\cdot)$ is not symmetric, the bilinear form occurring in (6.25) is not $a(e;\cdot,\cdot)$, but its adjoint bilinear form $a^*(e;\cdot,\cdot)$ which is defined by

$$ (6.26) \qquad a^*(e;u,v) = a(e;v,u) \quad \forall u,v \in H. \blacksquare $$

Let us see now that the introduction of $p(e)$, via (6.25), gives an elegant technique for eliminating $\delta u(e)$ in (6.23), modulo the solution of (6.25) (we emphasize on the fact that this technique is classical in *Optimal Control*, see e.g. [2]).

Taking $v = \delta u(e)$ in (6.25) we obtain from (6.25), (6.24) that

$$ (6.27) \qquad \left\{ \begin{aligned} \left(\frac{\partial J}{\partial u}(u(e),e), \delta u(e)\right)_H &= a(e;\delta u(e),p(e)) = \\ &= \left(\frac{\partial L}{\partial e}(e;p(e)), \delta e\right)_\xi - \left(\frac{\partial a}{\partial e}(e;u(e),p(e)), \delta e\right)_\xi \end{aligned} \right. $$

which, combined with (6.23), implies

$$ \delta j(e) = (\nabla j(e), \delta e)_\xi = \left( \left( \frac{\partial J}{\partial e}(u(e),e) + \frac{\partial L}{\partial e}(e;p(e)) \right) - \frac{\partial a}{\partial e}(e;u(e),p(e)), \delta e \right)_\xi $$

Finally we have proved that

$$ (6.28) \qquad \nabla j(e) = \frac{\partial J}{\partial e}(u(e),e) + \frac{\partial L}{\partial e}(e;p(e)) - \frac{\partial a}{\partial e}(e;u(e),p(e)). $$

**Application** : Let us consider again the problem (6.13)-(6.15) ; we assume that, in (6.18), $u_0$ is independent of $e$. We can easily prove that in this particular case we have
$$
(6.29) \quad \nabla j(e) = \frac{\partial J}{\partial e}(u(e),e) - \frac{\partial a}{\partial e}(e;u(e),p(e)),
$$

where $u(e)$ is the solution of (6.15) and $p(e)$ the solution of (6.25).

6.4.3. The finite dimensional case.

This section is of practical importance since effective calculations
are done in finite dimension. If $\mathcal{E}, H, H_o$ are finite dimensional with
for example $\mathcal{E} = \mathbb{R}^M, H_o = \mathbb{R}^N$, the canonical situation (6.13), (6.14),
(6.21) reduces to

$$
(6.30) \quad \text{Min} \ j(e) \quad \varepsilon \in E_f
$$

with

$$
(6.31) \quad j(e) = J(u(e),e)
$$

and $u(e)$ is solution of the linear system (the state equation)

$$
(6.32) \quad A(e)u(e) = f(e).
$$

In (6.32), $f(e) \in \mathbb{R}^N$ and $A(e)$ is a $N \times N$ matrix, positive definite
if $e \in E_f$ ; we denote by $a_{ij}(e)$ the general coefficient of $A(e)$.
Assuming that $J,A,f$ are differentiable functions of $e$, we have

$$
(6.32) \quad \delta j(e) = \left( \frac{\partial J}{\partial u}(u(e),e), \delta u(e) \right)_{\mathbb{R}^N} + \left( \frac{\partial J}{\partial e}(u(e),e), \delta e \right)_{\mathbb{R}^M},
$$

$$
(6.33) \quad A(e)\delta u(e) + \delta A(e)u(e) = \delta f(e) = \frac{df}{de}(e)\delta e.
$$

In (6.33), $\frac{df}{de}(e)$ is the $N \times M$ matrix whose general coefficient
$f'_{ik}(e)$ is defined by

$$
f'_{ik}(e) = \frac{\partial f_i}{\partial e_k}(e), \quad 1 \le i \le N, \ 1 \le k \le M,
$$
and where $\delta A(e)$ is a $N \times N$ matrix. It is in fact convenient to
observe that

$$
\delta A(e) u(e) = \left( \frac{dA}{de}(e) u(e) \right) \delta e
$$

where $\frac{dA}{de}(e)$ is the *third-order tensor* whose general coefficient
$a'_{ijk}(e)$ is defined by

$$
a'_{ijk}(e) = \frac{\partial a_{ik}}{\partial e_j}(e), \quad 1 \le i, k \le N, \ 1 \le j \le M;
$$

this implies that $\frac{dA}{de}(e)u(e)$ is a $N \times M$ matrix whose general coefficient
is (with $1 \le i \le N$, $1 \le j \le M$)

$$
\left( \frac{dA}{de}(e) u(e) \right)_{ij} = \sum_{k=1}^{N} \frac{\partial a_{ik}}{\partial e_j} (e) u_k(e).
$$

To eliminate $\delta u(e)$ in (6.32) we introduce the *adjoint equation*

$$
(6.34) \quad A^t(e)p(e) = \frac{\partial J}{\partial u}(u(e),u);
$$

it follows from (6.33), (6.34) that

$$
(6.35) \quad \left\{
\begin{aligned}
\left( \frac{\partial J}{\partial u}(u(e),e), \delta u(e) \right)_{\mathbb{R}^N} &= \left( A^t(e) p(e), \delta u(e) \right)_{\mathbb{R}^N} = \\
&= \left( A(e) \delta u(e), p(e) \right)_{\mathbb{R}^N} = \left( \frac{df}{de}(e) \delta e, p(e) \right)_{\mathbb{R}^N} - \left( \delta A(e) u(e), p(e) \right)_{\mathbb{R}^N}.
\end{aligned}
\right.
$$

From the above remark on $\delta A(e)u(e)$, it follows from (6.35) that

$$
(6.36) \quad \left\{
\begin{aligned}
& \left( \frac{\partial J}{\partial u}(u(e),e), \delta u(e) \right)_{\mathbb{R}^N} = \left( \left( \frac{df}{de}(e) \right)^t p(e), \delta e \right)_{\mathbb{R}^M} - \\
& \quad - \left( \left( \frac{dA}{de}(e)u(e) \right)^t p(e), \delta e \right)_{\mathbb{R}^M} .
\end{aligned}
\right.
$$

Finally from (6.32), (6.36) we obtain
$$ (6.37) \quad \nabla j(e) = \frac{\partial J}{\partial e}(u(e),e) + \left(\frac{df}{de}(e)\right)^t p(e) - \left(\frac{dA}{de}(e)u(e)\right)^t p(e). $$

**Remark 6.3**: To compute $\nabla j(e)$ we have to solve the two linear systems (6.32) and (6.34). From a practical point of view if one uses a **Gauss** or **Cholesky** method for the solution of these two systems, we clearly have only **one factorization** to do. Thus solving (6.32), (6.34) is almost as cheaper as solving (6.32) alone.

## 7. APPLICATIONS. (I) A FIRST EXAMPLE

In this section, and in Sec. 8, 9 we apply the methods of Sec. 4, 5, 6 to the numerical solution of three problems which are particular cases of those we have considered in Sec. 2, 3.

### 7.1. Description of the problem.

#### 7.1.1. Geometrical description :

The structure under consideration is an **axi-symmetric reservoir** containing water, up to a given level. The geometry of the meridian of the middle surface is shown on Fig. 7.1

<figure><img src="image_8.png" /><figcaption>Fig. 7.1.</figcaption></figure>
The meridian of the middle-surface is an arc of circle whose radius R = 40 meters ; the center is O. We have $\phi_o = \frac{3\pi}{10} = 54$ degrees, $\phi_1 = \frac{\pi}{4} = 45$ degrees. If L is the length of the meridian we have $L = \frac{3\pi}{10} R \approx 37.70$ meters.

### 7.1.2. Mechanical description

We assume that the shell is perfectly clamped at A ; this implies for $u(A)$ that

$$ (7.1) \qquad \left\{ \begin{array}{l} u_1 (A) = u_2 (A) = 0, \\ \frac{du_2}{ds} (A) = 0. \end{array} \right. $$

We assume moreover that the shell is free at B.

Concerning the forces acting on the shell, we neglect gravity and consider only the hydro-static pressure of the water filling the reservoir ; from these assumptions one obtains easily the function $F$ occurring in (3.4),(3.5). About the constitutive material of the structure we take $E=1$ for the Young modulus and $\nu = 0.2$ for the Poisson ratio.

**Remark 7.1** : The gravity can be easily introduced in $F$.

### 7.1.3. Definition of $E_f$. Choice of the cost function.

The convex $E_f$ of feasible thickness is defined by

$$ (7.2) \qquad 2.40 \text{ meters} \le e(s) \le 6.50 \text{ meters,} $$

$$ (7.3) \qquad \left| \frac{de}{ds} \right| \le \frac{1}{2}, $$

$$ (7.4) \qquad \int_{0}^{L} e(s)r(s)ds = \int_{0}^{L} e^{\phi}(s)r(s)ds $$
where

$$ (7.5) \qquad e^0(s) = 4 + 10 \left(\frac{s}{L}\right)^2 \left(1 - \frac{s}{L}\right)^2 . $$

The meridian S is oriented from B to A with $s(B) = 0$. As cost function we have used for this example the functional $j_1(\cdot)$ defined by (3.1), (3.10)-(3.12).

## 7.2. Solution of the approximate problem.

### 7.2.1. Discretization parameters.

We divide AB in 24 sub-arcs. Then using the approximations and notation of Sec. 5 we have

$$ (7.6) \qquad \begin{cases} \dim W_{1h} = 25, \\ \dim W_{2h} = 50. \end{cases} $$

We can easily take into account the boundary conditions (7.1).

### 7.2.2. Initialization of the Frank-Wolfe algorithm

As starting vector in algorithm (6.9)-(6.12) (see Sec. 6.3) we have taken an approximation of the function $e^0$ defined by (7.5).

## 7.3. Convergence of the iterative method.

We may consider that the convergence was effective after 28 iterations (corresponding to 173 seconds of IBM 370-165).

On the table below we have shown the variation of the cost functional with respect to the number of iterations.
<table><caption>Table 7.1</caption><thead><tr><td>Iteration</td><td>Cost Function</td></tr></thead><tbody><tr><td>0</td><td>0.2775</td></tr><tr><td>1</td><td>0.2055</td></tr><tr><td>2</td><td>0.1852</td></tr><tr><td>3</td><td>0.1815</td></tr><tr><td>5</td><td>0.1769</td></tr><tr><td>10</td><td>0.1735</td></tr><tr><td>15</td><td>0.1722</td></tr><tr><td>20</td><td>0.1717</td></tr><tr><td>25</td><td>0.1714</td></tr></tbody></table>

### 7.4. Mechanical analysis of the results.

The maximum value of the compression stresses has been reduced from 2410 tons/m² to 1490 tons/m². The maximum value of the traction stresses has been reduced from 957 tons/m² to 625 tons/m². The maximum value of the radial displacement has been reduced of about 50%.

About the aspect of the optimal thickness law, it agrees in fact with mechanical intuition, since thickness is reinforced at the clamping point A (s=L) and at the application-point of the resultant of hydrostatic pressure forces (s ≈ 21 meters). Since the weight of the structure is given, local reinforcements of thickness of the shell must be compensated by reduction elsewhere; this is very clear on Fig. 7.2, 7.3 where are shown $e^0(s)$, $e^{15}(s)$, $e^{28}(s)$.

On Fig. 7.4, 7.5, 7.6 and 7.7 we have shown the internal meridional and circumferential stresses at iterations n=0,15,28.
<figure><img src="image_9.png" /><figcaption>Figure 7.2</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7.3</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 7.4</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 7.5</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 7.6</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 7.7</figcaption></figure>

## 8. APPLICATIONS. (II) A SECOND EXAMPLE

The structure that we consider in this example is the _pressure vessel_ shown on Fig. 8.1; cylinder and spherical caps have a radius of 5 meters.

The _Young modulus_ is $E=1$ and the _Poisson ratio_ is $\nu=0.29$.
<figure><img src="image_15.png" /><figcaption>Figure 8.1</figcaption></figure>

The **load** consists of an **internal pressure** of 5 MN/m².

The cost function is $j_2(\cdot)$ (cf. (3.3)) with $\beta_i(s) = \beta_e(s) = 1$ for all $s$. The convex $E_f$ is defined by

$$ (8.1) \qquad 0.15 \text{ meter} \le e(s) \le 0.40 \text{ meter} \quad \forall s, $$

$$ (8.2) \qquad \left| \frac{de}{ds} \right| \le \frac{1}{2} \quad \forall s, $$

$$ (8.3) \qquad \int_S e(s)r(s)ds = \int_S e^\circ(s)r(s)ds $$

with

$$ (8.4) \qquad e^\circ(s) = 0.25 \text{ meter } \forall s. $$
To approximate the structural problem by the finite element method of Sec. 5, we have used a subdivision of the meridian S in 34 subintervals.

After 30 iterations of the Frank-Wolfe algorithm (6.9)-(6.12) (started from $e^°$ defined by (8.4)) the maximal value of the stresses is reduced from $105 \ \text{MN/m}^2$ to $87 \ \text{MN/m}^2$; moreover, as seen on Fig. 8.2, 8.3, the stress discontinuities at the interface caps/cylinder are substantially reduced.

The computational time corresponding to the above run is 3 minutes on IBM 370-168.

<figure><img src="image_16.png" /><figcaption>Figure 8.2</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 8.3</figcaption></figure>
## 9. APPLICATIONS. (III) A THIRD EXAMPLE

The structure considered in this example is the vessel shown on Fig. 9.1. The cylindrical part is clamped at its top.

<figure><img src="image_18.png" /><figcaption>Figure 9.1</figcaption></figure>

The vessel is full with a liquid of density 0.85. Thermal stresses are neglected. The material constants are still $E=1$, $v=0.29$. The cost function is $j_2(\cdot)$ with $\beta_i(s) = \frac{1}{2} \forall s$, $\beta_e(s) = 1 \forall s$.
The convex $E_f$ is defined by

$$ (9.1) \qquad 0.032 \text{ meter} \le e(s) \le 0.048 \text{ meter} $$

$$ (9.2) \qquad \left| \frac{de}{ds} \right| \text{ is not bounded} $$

$$ (9.3) \qquad \int_S e(s)r(s)ds = \int_S e^\circ(s)r(s)ds $$

with

$$ (9.4) \qquad e^\circ(s) = 0.04 \text{ meter } \forall s. $$

After 15 iterations (1 minute on IBM 370/168) the maximal value of the stresses has been reduced from 68 MN/m² to 55 MN/m² (the starting thickness is $e^°$ defined by (9.4)). Of course the *cylindrical part* of the structure becomes *thinner* and the *bottom* becomes *thicker*. On Figures 9.2, 9.3, 9.4, 9.5, 9.6 we have shown thickness and stress distributions at iterations n=0 and 15.

<figure><img src="image_19.png" /><figcaption>Figure 9.2</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 9.3</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 9.4</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 9.5</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 9.6</figcaption></figure>
## 10. COMMENTS AND CONCLUSION

One may object to the above approach to Structural Optimization, the fact that the structures under consideration have a *given weight*. Indeed in many applications one tries to minimize the weight, with the requirement that some functionals of stresses and/or displacements stay below given yields. In this type of problems the cost function is simple, but in the set of relations defining $E_f$ appear u and/or $\sigma$; hence we have to take into account *state constraints*. It is well known in Optimal Control that state constraints usually bring many complications in the numerical solution of corresponding optimization problems. It is to overcome these difficulties that we have worked at given weight with a cost functional containing *stress terms* (we can also easily handle *displacement terms*). Indeed if the weight is not given it is always possible to take it as a parameter and to solve a sequence of constant weight problems.

There exists a large number of papers, reports, etc... on Structural Optimum Design and it is not possible to mention all of them. We would like to mention PRAGER [8], GALLAGHER-ZIENKIEWICZ [9], FLEURY-FRAEJS DE VEUBEKE [10] and the bibliographies therein.

The optimization of *axi-symmetric elastic membranes* is considered in KOSTEM [11].
Compared to the above works we have systematically used the
point of view of *Optimal Control* which has the great advantage
of putting a *hierarchy* between the *control variable* - in fact the
*design variable* - (here the *thickness e*) and the *state variable*
(here the *displacement u*). The Optimal Control approach which
is very general (it can handle *time dependent* and/or *non linear*
*state equations*) leads in a very natural way to a great variety
of *descent algorithms* in which the *differential* of the cost
functional is obtained via the solution of an *adjoint state
equation*.

Among possible generalizations of the methods of this paper,
we would like to mention that we have in view Optimum Structural
Design problems, involving non axi-symmetric shells, in which
the approximation is done by the methods in ARGYRIS-HAASE-
MALEJANNAKIS [12], CIARLET [13].

Shorter versions of the present paper may be found in BOISSERIE-
GLOWINSKI [14], [15].

REFERENCES

[1] O.C. ZIENKIEWICZ, *The finite element method in engineering science*, Mc Graw Hill, London, (1971).

[2] J.L. LIONS, *Contrôle Optimal des systèmes gouvernés par des équations aux dérivées partielles*, Dunod-Gauthier Villars, Paris (1968).

[3] J.T. ODEN and J.N. REDDY, *An introduction to the mathematical theory of finite elements*, Wiley-Interscience, New-York (1976).
[4] G. STRANG and G. FIX, *An analysis of the finite element method*, Prentice Hall, Englewood Cliffs, N.J. (1973).

[5] K. YOSIDA, *Functional Analysis*, 3rd Edn., Springer-Verlag, New-York (1971).

[6] J. CEA, *Optimisation : Théorie et Algorithmes*, Dunod, Paris (1971).

[7] D.J. WILDE and C.S. BEIGHTLER, *Foundations of Optimization*, Prentice Hall, Englewood Cliffs, N.J. (1967).

[8] W. PRAGER, Methods of structural optimization, in *Computing Methods in Applied Science and Engineering*, R. Glowinski and J.L. Lions Ed., Lecture Notes in Computer Science, Vol. 10, Springer-Verlag, Berlin, (1974), pp. 1-19.

[9] R.M. GALLAGHER, O.C. ZIENKIEWICZ (Ed.), *Optimum Structural Design : Theory and Applications*, J. Wiley and Sons, London (1973).

[10] C. FLEURY and B. FRAEJS DE VEUBEKE, Structural Optimization, in *Optimization Techniques, I.F.I.P. Technical Conference*, G.I. Marchuk Ed., Lecture Notes in Computer Science, Vol. 27, Springer-Verlag, (1975), pp. 314-326.

[11] G.N. KOSTEM, On the variational optimization of rotational membranes, *Proceedings of the International Conference on Variational Methods in Engineering*. Southampton, September 1972.

[12] J.H. ARGYRIS, M. HAASE and G.A. MALEJANNAKIS, Natural geometry of surfaces with specific reference to the matrix displacement analysis of shells, I, II and III. *Proc. Kon. Ned. Akad. Wetensch.*, B. 76 (1973), 361-410.
[13] P.G. CIARLET, *Numerical Analysis of the finite element method*, Presses de l'Université de Montréal (1976).

[14] J.M. BOISSERIE and R. GLOWINSKI, *Optimization de la loi d'épaisseur pour une coque mince de révolution*, in *Computing Methods in Applied Sciences and Engineering*, R. Glowinski and J.L. Lions Ed., Lecture Notes in Economics and Mathematical Systems, Vol. 134, Springer-Verlag, Berlin, (1976), pp. 258-275.

[15] J.M. BOISSERIE and R. GLOWINSKI, *Optimization of thickness law for a pressure vessel*, in *Proceedings of the 3rd International Conference on Pressure Vessel Technology, Part 1 Analysis, Design and Inspection*, A.S.M.E., New-York, (1977), pp. 275-281.