On the numerical solution of a micropolar continuum model Reijo Kouhia1 and Antti H. Niemi2 (1)

Tampere University of Technology Department of Engineering Design [email protected] (2) Aalto University, School of Engineering Department of Civil and Structural Engineering [email protected] NSCM-26 Oslo, 23-25 October 2013

Content • Equilibrium equations • Isotropic, linear elastic constitutive equations • Comprehensible set of material parameters • Conformally invariant microcurvature • Size effect - pure bending of a beam with circular cross-section • Numerical solution • Example

NSCM-26 Oslo 2013

2

Equilibrium equations Force balance

Moment balance

∂σji + ρbi = 0, ∂xj ∂µji + ρci + ijk σjk = 0, ∂xj

σ force stress tensor µ moment stress tensor b body force per unit mass c body moment per unit mass ijk alternating tensor ρ the mass density

NSCM-26 Oslo 2013

3

Constitutive equations Linear material model (γ)

σij = Cijklγkl,

(κ)

µij = Cijklκkl,

the Cosserat’s first strain tensor γij and microcurvature tensor κij are γij = uj,i − kij ϕk ,

and κij = ϕj,i,

ϕ is the microrotation field The material stiffness tensors for an isotropic solid can be expressed as (γ)

Cijkl = λδij δkl + µ(δik δjl + δilδjk ) + µc(δik δjl − δilδjk ), (κ)

Cijkl = αδij δkl + β(δik δjl + δilδjk ) + γ(δik δjl − δilδjk ). NSCM-26 Oslo 2013

4

Comprehensible set of material parameters s the characteristic length in torsion

`t = s

the characteristic length in bending the coupling number and the polar ratio

0 ≤ N ≤ 1,

`b =

β , µ β+γ , 4µ

r µc N= , µ + µc ψ=

β 1 . β + 2α

3 and 0 ≤ ψ ≤ , 2 NSCM-26 Oslo 2013

5

Conformally invariant curvature state Neff et al. 2009, 2010 µij = (α + 32 β)κkk δij + 2βκ0(ij) + 2γκ[ij] γ(ij) = 21 (γij + γji),

γ[ij] = 21 (γij − γji).

The curvature tensor is symmetric and deviatoric γ = 0,

β = µ`2t,

α = − 23 µ`2t.

Four material parameters λ, µ, µc, `t.

NSCM-26 Oslo 2013

6

Size effect - pure bending of a beam Circular cross-section with radius a (a)

(b) 2.2

1.2

2 Ωb (a → 0)

1.25

Ωb

1.15 1.1 1.05 1

1.8 1.6 1.4 1.2

0

2

4

6 a/ℓb

8

10

1 10−4

10−2

100 µc /µ

102

104

(a) From Figure bottom to Ratio top µ /µ =rigidity 0.01, 0.03, . . . ,versus 0.09,standard 0.1. Solid conformally 1. (a) of cbending of 0.02, the polar continuum Cauchy’slines: continuum model Ωb as a function of the radius 2 a of the beam’s cross-section with different values of the coupling invariant curvature case ( β = 4µ` 0);0.02, dashed γ= 0.01β . invariant curvature modulus µc ; from bottom to top µbc,/µγ== 0.01, 0.03, . . . ,lines 0.09, 0.1. Conformally 2 case 0) (β =for 4µℓthe 0) shown by solid lines and the case with parameters b, γ = (b) Ωb(a → conformally invariant curvature state. γ = 0.01β indicated by dashed lines. (b) Limit value of Ωb when a → 0 as a function of the coupling modulus µc for the conformally invariant curvature state. NSCM-26 Oslo 2013

7

Numerical solution The virtual work expression Z V

Z V



T



 T −∇ · σ − ρb · δu + −∇ · µ −  : σ − ρc · δφ dV = 0


 T T σ : (∇δu) − skew(δϕ) + µ : (∇δϕ) dV Z Z − (ρb · δu + ρc · δϕ) dV − (t · δu + m · δϕ) dS = 0 V

S

Standard C0-interpolation for both u and ϕ.

NSCM-26 Oslo 2013

8

modulus µc ; from bottom to top µc /µ = 0.01, 0.02, 0.03, . . . , 0.09, 0.1. Conformally invariant curvature case (β = 4µℓ2b , γ = 0) shown by solid lines and the case with parameters γ = 0.01β indicated by dashed lines. (b) Limit value of Ωb when a → 0 as a function of the coupling modulus µc for the conformally invariant curvature state.

Example - cantilever beam (b)

1

1

0.8

0.8 vtip /vref

vtip /vref

(a)

0.6 0.4 0.2 0 10−8

0.6 0.4 0.2

10−4

100 µc /µ

104

108

0 10−3 10−2 10−1 100 101 102 ℓb /h

103

Figure 2. (a) The effect oftwo the coupling µc to theelements tip displacement the cantilever beam. with The trilinear (a)`b = h/200 . The upper curvesmodulus triquadratic andof the lower ones 3 reference value is the tip deflection according to the Euler-Bernoulli beam model: vref = F L /3EI. The elements. material Solid length lines:in bending rotations freewhere at the edge, dashed lines: rotations are is ℓb are = h/200, h is clamped the cross-section height. The upper two curves supressed.correspond to solutions with triquadratic elements and the lower ones with trilinear elements. Solid lines indicate solutions where the rotations are free at the clamped edge and the dashed line where they are supressed. (b) The effect of internal scale ℓb FE on the FE solution: (b) The effect of internal length scalelength `b on the solution: µcµc== µ. µ.

NSCM-26 Oslo 2013

9

Concluding remarks • Stable and locking free FE-formulation needed for the micropolar model.

NSCM-26 Oslo 2013

10

Thanks!

Edvard Munch, Anxiety, 1894

NSCM-26 Oslo 2013

11

Equilibrium equations in terms of displacements and rotations

(λ + µ − µc)grad divu + (µ + µc)∆u + 2µccurlϕ + ρb = 0

(α + β − γ)grad divϕ + (β + γ)∆ϕ − 4µcϕ + 2µccurlu + ρc = 0

∆u = div gradu

NSCM-26 Oslo 2013

12

Constitutive parameters

e = [γ11, γ22, γ33, γ12, γ21, γ23, γ32, γ31, γ13, κ11, κ22, κ33, κ12, κ21, κ23, κ32, κ31, κ13]T , C = diag





C1 C2 C2 C2 C3 C4 C4 C4 ,     α + 2β α α λ + 2µ λ λ , , α α + 2β α λ λ + 2µ λ C3 =  C1 =  α α α + 2β λ λ λ + 2µ     µ + µc µ − µc β+γ β−γ C2 = , C4 = . µ − µc µ + µc β−γ β+γ µ > 0,

3λ + 2µ > 0,

µc > 0,

β > 0,

3α + 2β > 0,

γ > 0.

NSCM-26 Oslo 2013

13

Stiffness ratio Ωb in pure bending of a beam with circular cross-section with radius a 2

Ωb = 1 +

2

8N 1 − (β − γ)/(β + γ) [(β − γ)/(β + γ) + ν] + 2 1+ν (δa) ζ(δa) + 8N 2(1 − ν) 

 ,

δ 2 = 4µcµ/(β + γ)(µ + µc) x2(xI0(x) − I1(x)) , ζ(x) = xI0(x) − 2I1(x) I0, I1 the modified Bessel functions of the first kind Krishna-Reddy and Venkatasubramanian 1978 NSCM-26 Oslo 2013

14