Forces and Moments Calculation

Integration Points Throughout the Thickness

The integration is performed using n equally spaced integration points throughout the thickness. The method used assumes a linear variation of stresses between integrations points:(1)

N_{x} = t \sum_{k = 1}^{n} w_{k}^{N} σ_{x k}^{p a}

(2)

M_{x} = t^{2} \sum_{k = 1}^{n} w_{k}^{M} σ_{x k}^{p a}

Table 1 compares the coefficients used to the classical Newton quadrature in case of 3 integration points.

Table 1. w^N for 3 Integration Points
Coefficients	w₁	w₂	w₃
Radioss	0.250	0.500	0.250
Simpson	0.166	0.666	0.166

Table 2. w^M for 3 Integration Points
Coefficients	w₁	w₂	w₃
Radioss	-0.083	0.	0.083
Simpson	-0.083	0.	0.083

Table 3. Gauss Integration Scheme
Number of Points	Position	Weight
1	±0.0000	2.0000
2	±0.5774	1.0000
3	0.0000 ±0.7746	0.8889 0.5556
4	±0.8611 ±0.3400	0.6521 0.3479
5	±0.9062 ±0.5385 0.0000	0.2369 0.4786 0.5689
6	±0.9325 ±0.6612 ±0.2386	0.1713 0.3608 0.4679
7	±0.9491 ±0.7415 ±0.4058 0.0000	0.1295 0.2797 0.3818 0.4180
8	±0.9603 ±0.7967 ±0.5255 ±0.1834	0.1012 0.2224 0.3137 0.3627
9	±0.9681 ±0.8360 ±0.6134 ±0.3243 0.0000	0.0813 0.1806 0.2606 0.3123 0.3302
10	±0.9739 ±0.8650 ±0.6794 ±0.4334 ±0.1489	0.0667 0.1495 0.2191 0.2693 0.2955

Table 4. Lobatto Integratin Scheme
Number of Points	Position	Weight for Membrane w^N	Weight for Bending w^M
1	0.0000	1.0000	0.0000
2	±0.5000	0.5000	±0.0833
3	±0.5000 0.0000	0.2500 0.5000	±0.0833 0.0000
4	±0.5000 ±0.1667	0.1667 0.3333	±0.0648 ±0.0556
5	±0.5000 ±0.2500 0.0000	0.1250 0.2500 0.2500	±0.0521 ±0.0625 0.0000
6	±0.5000 ±0.3000 ±0.1000	0.1000 0.2000 0.2000	±0.0433 ±0.0600 ±0.0200
7	±0.500 ±0.3333 ±0.1667 0.0000	0.0833 0.1667 0.1667 0.1667	±0.0370 ±0.0556 ±0.0278 0.0000
8	±0.5000 ±0.3750 ±0.2500 ±0.1250	0.0714 0.1429 0.1429 0.1429	±0.0323 ±0.0510 ±0.0306 ±0.0102
9	±0.5000 ±0.3750 ±0.2500 ±0.1250 0.0000	0.0625 0.1250 0.1250 0.1250 0.1250	±0.086 ±0.0469 ±0.0313 ±0.0156 0.0000
10	±0.5000 ±0.3889 ±0.2778 ±0.1667 ±0.0555	0.0556 0.1111 0.1111 0.1111 0.1111	±0.0257 ±0.0432 ±0.0309 ±0.0185 ±0.0062

For shell elements, integration points through the thickness are almost Lobatto points.

Global Plasticity Algorithm

In the case of global plasticity, the forces and moments are computed directly. The algorithm is activated by specifying the number of integration points through the thickness as zero. The first step is an obvious elastic calculation:(3)

{\sum^{e l}} = {\sum (t)} + L {\dot{E}} Δ t

The yield criterion used is the uncoupled Iliouchine ¹ form:(4)

F = \frac{{\bar{N}}^{2}}{t^{2}} + \frac{16 {\bar{M}}^{2}}{t^{4}} - σ_{y i e l d}^{2} \leq 0

with(5)

{\bar{N}}^{2} = N_{x}^{2} + N_{y}^{2} - N_{x} N_{y} + 3 N_{x y}^{2}

(6)

{\bar{M}}^{2} = M_{x}^{2} + M_{y}^{2} - M_{x} M_{y} + 3 M_{x y}^{2}

Where,

$N_{x} = \int_{- t / 2}^{t / 2} σ_{x}^{p a} d z$	$M_{x} = \int_{- t / 2}^{t / 2} σ_{x}^{p a} z d z$
$N_{y} = \int_{- t / 2}^{t / 2} σ_{y}^{p a} d z$	$M_{y} = \int_{- t / 2}^{t / 2} σ_{y}^{p a} z d z$
$N_{x y} = \int_{- t / 2}^{t / 2} σ_{x y}^{p a} d z$	$M_{x y} = \int_{- t / 2}^{t / 2} σ_{x y}^{p a} z d z$

An extension of Iliouchine criterion for isotropic hardening is developed here. The yield surface can be expressed as:(7)

f = {\begin{matrix} {N} \\ {M} \end{matrix}}^{t} [F] {\begin{matrix} {N} \\ {M} \end{matrix}} - {(σ_{y}^{0} β)}^{2} = 0

with(8)

[F] = [\begin{matrix} \frac{1}{h^{2}} [A] & \frac{1}{2 \sqrt{3} (s γ h \frac{h^{2}}{6})} [A] \\ s y m & \frac{1}{{(γ \frac{h^{2}}{6})}^{2}} [A] \end{matrix}]

and (9)

s = \frac{| {N}^{t} [A] {M} |}{{N}^{t} [A] {M}}

Where, $β$ and $γ$ are scalar material characteristic constants, function of plastic deformation. They can be identified by the material hardening law in pure traction and pure bending.

In pure traction:(10)

f = \frac{N^{2}}{h^{2}} - {(σ_{y}^{0} β)}^{2} = 0 \to β = \frac{σ_{y} (ε^{p})}{σ_{y}^{0}}

In pure bending: (11)

f = \frac{M^{2}}{{(γ h^{2} / 6)}^{2}} - {(σ_{y}^{0} β)}^{2} = 0 \to γ = \frac{M}{σ_{y} h^{2} / 6}

If no hardening law in pure bending is used, $γ$ is simply computed by $γ = \frac{σ_{y} / E + \frac{3}{2} ε^{p}}{σ_{y} / E + ε^{p}}$ varying between 1.0 and 1.5.

The plasticity flow is written using the normality law:(12)

{\begin{matrix} {d ε_{p}} \\ {d Χ_{p}} \end{matrix}} = d λ \frac{\partial f}{\partial {\begin{matrix} {N} \\ {M} \end{matrix}}} = 2 d λ [F] {\begin{matrix} {N} \\ {M} \end{matrix}}

The equivalent plastic deformation

ε^{P}

is proportional to the plastic work. Its expression is the same as in the case of traction:(13)

σ_{y}^{0} β d ε^{p} = {\begin{matrix} {d ε_{p}} \\ {d χ_{p}} \end{matrix}}^{t} {\begin{matrix} {N} \\ {M} \end{matrix}} = 2 d λ {(σ_{y}^{0} β)}^{2}

This leads to:(14)

d ε^{p} = 2 d λ σ_{y}^{0} β and d β = \frac{H}{σ_{y}^{0}} d ε^{p}

Where, $H$ is the plastic module. The derivative of function $f$ in Equation 7 is discontinuous when ${N}^{t}$ ${A} {M}$ =0. This can be treated when small steps are used by putting s=0 as explained in ².

Then the derivative of

f

with respect to

d ε^{P}

(

\frac{\partial f}{\partial d ε^{P}}

) is carried out. The derived equation is nonlinear in internal efforts and is resolved by Newton-Raphson:(15)

{\begin{matrix} {N} \\ {M} \end{matrix}}_{n + 1} = {\begin{matrix} {N} \\ {M} \end{matrix}}_{*} - 2 d λ [D] [F] {\begin{matrix} {N} \\ {M} \end{matrix}}_{n + 1}

Where,

[D]

is the elastic stiffness matrix and:(16)

{\begin{matrix} {N} \\ {M} \end{matrix}}_{*} = {\begin{matrix} {N} \\ {M} \end{matrix}}_{n} - [D] {\begin{matrix} {d ε} \\ {d χ} \end{matrix}}

¹ Iliouchine A., “Plasticity”, Edition Eyrolles Paris, 1956.

² Crisfield M.A., “Nonlinear finite element analysis of solids and structures”, J. Wiley, Vol. 2, 1997.