plasticity

<< Click to Display Table of Contents >>

Navigation:  Sample Problems > Applications > Stress >

plasticity

Previous pageReturn to chapter overviewNext page

{ PLASTICITY.PDE

 

This is a moving-mesh model of large displacement plasticity.

It uses a pseudo-time evolution to relax into its final configuration.

The "time" behavior is not intended to model the actual temporal behavior of such a loaded bar.

The Elasticity is at first strongly diffused to avoid erratic hopping between elastic and plastic states.

Once the shape is established, we decrease the diffusion to focus in on the transition surface.

}

 

TITLE     'Beam, Shear-Loaded at End'

 

SELECT

regrid=off         ! Don't allow mesh refinement.  It confuses the elasticity transition

prefer_stability   ! allow Newton steps as necessary for convergence

order=2           ! use quadratic FEM interpolation

 

VARIABLES

  u(1e-4)

  v(1e-3)

  xm= move(x)

  ym= move(y)

  E(10e9)

 

DEFINITIONS

  L=1.0

  h=0.2

  mu=0.3

  s_pl=1e9       { Stress limit for elasticity }

  E0   =200e9

 

  C=E/(1-mu^2)     G=E/(2*(1+mu))

  uv=vector(u,v)     uvm=magnitude( uv)

  ex=dx(u)     ey=dy(v)     exy=dx(v)+ dy(u)

 

  !>> Mises for rest elasticity determines where the elasticity switches.

  C0=E0/(1-mu^2)     G0=E0/(2*(1+mu))

  sx0=(C0*(ex+ mu*ey))     sy0=(C0*(mu*ex+ ey))     sxy0=(G0*exy)

  ez0=-mu*(sx0+ sy0)/E0

  mises0=sqrt( 0.5*( (sx0-sy0)^2+ sx0^2+ sy0^2)+ 3*sxy0^2)

!  E1 = SWAGE(mises0-s_pl, E0, E0/20, s_pl/10)

  E1 = if(mises0 < s_pl) then E0 else E0/20

 

  sx=C*(ex+ mu*ey)     sy=C*(mu*ex+ ey)     sxy=G*exy

  ez=-mu*(sx+ sy)/E

  p_ang=0.5* arctan( 2*sxy/(sx-sy) )                           { Radians }

  sxp0=(sx+sy)/2+ (sx-sy)/2*cos(2*p_ang)+ sxy*sin(2*p_ang)

  syp0=(sx+sy)/2- (sx-sy)/2*cos(2*p_ang)- sxy*sin(2*p_ang)

  { Test for highest algebraic value: }

  p_angl=     if sxp0>syp0 then p_ang     else p_ang+ pi/2

  sxp=     if sxp0>syp0 then sxp0     else syp0               { sx' }

  syp=     if sxp0>syp0 then syp0     else sxp0               { sy' }

  p_angle=p_angl*180/pi                                       { Degrees }

  mises=sqrt( 0.5*( (sx-sy)^2+ sx^2+ sy^2)+ 3*sxy^2)

  energy_d=(1/2)*( sx*ex+ sy*ey+ sxy*exy)               { Energy density }

 

 

  ! report the plastic distribution

  plastic=   if mises0>1e9 then 1 else 0

 

  ! start out with a heavily smeared distribution of E,

  !   then gradually narrow the transition by reducing the smear coefficient

  smear = (h/(t+2))^2

 

  endload = -6e7*min(t,1)

  endpoint = POINT(L,0)

 

INITIAL VALUES

 E=E0

 

EQUATIONS

  ! instantaneous evolution of displacement

  u:        dx( sx)+ dy( sxy)=0

  v:        dx( sxy)+ dy( sy)=0

  ! mesh motion follows displacement

  xm:       dt(xm) = dt(u)

  ym:       dt(ym) = dt(v)

  ! relax the effective elasticity towards the computed value

  E:        dt(E) = smear*div(grad(E)) + 10*(E1 - E)

 

BOUNDARIES

region 'steel'

  mesh_spacing = h/20*(1+3*x)   ! denser mesh near mount

 

  start 'outer' (0,-h/2)

 

  LABEL 'bottom'               ! label a piece of the boundary

  load(u)=0     load(v)=0 line to (L,-h/2)

  ENDLABEL 'bottom'

 

  load(u)=0     load(v)=endload line to(endpoint) to (L,h/2)

 

  LABEL 'top'                   ! label another piece of the boundary

  load(u)=0     load(v)=0     line to (0,h/2)

  ENDLABEL 'top'

 

  value(u)=0     value(v)=0 velocity(xm)=0 velocity(ym)=0

  mesh_spacing=h/40

  line to close

 

FRONT (mises0-s_pl, s_pl/10)

resolve(E)

 

TIME 0 by 1e-5 to 5

 

MONITORS

for cycle=10

  contour(E) painted as "Effective Elasticity"

  contour(E1) painted as "Computed Elasticity"

  contour( mises0) painted

  contour(plastic) painted    

 

PLOTS

  for cycle=50

  grid( x, y)

  contour( u) contour(dt(u))

  contour( v) contour(dt(v))

  contour(E)  painted as "Effective Elasticity"

  contour(E1) painted as "Computed Elasticity"

  contour( mises0) painted

  elevation(mises0) on 'top'

  contour( plastic) painted

  contour( sxy) painted

  elevation(y, h/2+v) on 'top'   as "Validate Displacement (V::Y)" ! verify that mesh and displacement agree

  elevation(y, -h/2+v) on 'bottom'   as "Validate Displacement (V::Y)"

  history( sxy) at  (endpoint)

  history( sxy) at  (0,h/2)

  history(endload)

  history(dt(v)) at (endpoint)

  history(deltat) LOG

  history(smear)

END