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)
Nx=tnk=1wNkσpaxkNx=tnk=1wNkσpaxk
(2)
Mx=t2nk=1wMkσpaxkMx=t2nk=1wMkσpaxk
Table 1 compares the coefficients used to the classical Newton quadrature in case of 3 integration points.
Table 1. wN for 3 Integration Points
Coefficients w1 w2 w3
Radioss 0.250 0.500 0.250
Simpson 0.166 0.666 0.166
Table 2. wM for 3 Integration Points
Coefficients w1 w2 w3
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

wN

Weight for Bending

wM

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)
{el}={(t)}+L{˙E}Δt{el}={(t)}+L{˙E}Δt
The yield criterion used is the uncoupled Iliouchine 1 form:(4)
F=ˉN2t2+16ˉM2t4σ2yield0F=¯¯¯¯N2t2+16¯¯¯¯M2t4σ2yield0
with(5)
ˉN2=N2x+N2yNxNy+3N2xy¯¯¯¯N2=N2x+N2yNxNy+3N2xy
(6)
ˉM2=M2x+M2yMxMy+3M2xy¯¯¯¯M2=M2x+M2yMxMy+3M2xy
Where,
Nx=t/2t/2σpaxdzNx=t/2t/2σpaxdz Mx=t/2t/2σpaxzdzMx=t/2t/2σpaxzdz
Ny=t/2t/2σpaydzNy=t/2t/2σpaydz My=t/2t/2σpayzdzMy=t/2t/2σpayzdz
Nxy=t/2t/2σpaxydzNxy=t/2t/2σpaxydz Mxy=t/2t/2σpaxyzdzMxy=t/2t/2σpaxyzdz
An extension of Iliouchine criterion for isotropic hardening is developed here. The yield surface can be expressed as:(7)
f={{N}{M}}t[F]{{N}{M}}(σ0yβ)2=0
with(8)
[F]=[1h2[A]123(sγhh26)[A]sym1(γh26)2[A]]
and (9)
s=|{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=N2h2(σ0yβ)2=0β=σy(εp)σ0y
In pure bending: (11)
f=M2(γh2/6)2(σ0yβ)2=0γ=Mσyh2/6

If no hardening law in pure bending is used, γ is simply computed by γ=σy/E+32εpσy/E+εp varying between 1.0 and 1.5.

The plasticity flow is written using the normality law:(12)
{{dεp}{dΧp}}=dλf{{N}{M}}=2dλ[F]{{N}{M}}
The equivalent plastic deformation εP is proportional to the plastic work. Its expression is the same as in the case of traction:(13)
σ0yβdεp={{dεp}{dχp}}t{{N}{M}}=2dλ(σ0yβ)2
This leads to:(14)
dεp=2dλσ0yβanddβ=Hσ0ydε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 2.

Then the derivative of f with respect to dεP ( fdεP ) is carried out. The derived equation is nonlinear in internal efforts and is resolved by Newton-Raphson:(15)
{{N}{M}}n+1={{N}{M}}2dλ[D][F]{{N}{M}}n+1
Where, [D] is the elastic stiffness matrix and:(16)
{{N}{M}}*={{N}{M}}n[D]{{dε}{dχ}}
1
Iliouchine A., “Plasticity”, Edition Eyrolles Paris, 1956.
2
Crisfield M.A., “Nonlinear finite element analysis of solids and structures”, J. Wiley, Vol. 2, 1997.