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 k = 1 n w k N σ x k p a
(2)
M x = t 2 k = 1 n w k M σ x k p a
表 1 compares the coefficients used to the classical Newton quadrature in case of 3 integration points.
表 1. wN for 3 Integration Points
Coefficients w1 w2 w3
Radioss 0.250 0.500 0.250
Simpson 0.166 0.666 0.166
表 2. wM for 3 Integration Points
Coefficients w1 w2 w3
Radioss -0.083 0. 0.083
Simpson -0.083 0. 0.083
表 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

表 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)
{ e l } = { ( t ) } + L { E ˙ } Δ t
The yield criterion used is the uncoupled Iliouchine 1 form:(4)
F = N ¯ 2 t 2 + 16 M ¯ 2 t 4 σ y i e l d 2 0
with(5)
N ¯ 2 = N x 2 + N y 2 N x N y + 3 N x y 2
(6)
M ¯ 2 = M x 2 + M y 2 M x M y + 3 M x y 2
Where,
N x = t / 2 t / 2 σ x p a d z M x = t / 2 t / 2 σ x p a z d z
N y = t / 2 t / 2 σ y p a d z M y = t / 2 t / 2 σ y p a z d z
N x y = t / 2 t / 2 σ x y p a d z M x y = 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 = { { N } { M } } t [ F ] { { N } { M } } ( σ y 0 β ) 2 = 0
with(8)
[ F ] = [ 1 h 2 [ A ] 1 2 3 ( s γ h h 2 6 ) [ A ] s y m 1 ( γ h 2 6 ) 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 = N 2 h 2 ( σ y 0 β ) 2 = 0 β = σ y ( ε p ) σ y 0
In pure bending: (11)
f = M 2 ( γ h 2 / 6 ) 2 ( σ y 0 β ) 2 = 0 γ = M σ y h 2 / 6

If no hardening law in pure bending is used, γ is simply computed by γ = σ y / E + 3 2 ε 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 } } = 2 d λ [ 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)
σ y 0 β d ε p = { { d ε p } { d χ p } } t { { N } { M } } = 2 d λ ( σ y 0 β ) 2
This leads to:(14)
d ε p = 2 d λ σ y 0 β and d β = H σ y 0 d ε p

Where, H MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyyaaaa@36DC@ is the plastic module. The derivative of function f in 式 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 ( f d ε 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 } } 2 d λ [ 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.