Resolve

title "RESOLVE shear stress in bent bar"

{

This is a test problem from Timoshenko: Theory of Elasticity, p41

The RESOLVE statement has been added to force regridder to resolve the

shear stress.

}

select

   elevationgrid=500

variables

   U   { X-displacement }

   V    { Y-displacement }

definitions

   L = 1               { Bar length }

   hL = L/2

   W = 0.1             { Bar thickness }

   hW = W/2

   eps = 0.01*L

   I = 2*hW^3/3        { Moment of inertia }

   nu = 0.3            { Poisson's Ratio }

   E  = 2.0e11         { Young's Modulus for Steel (N/M^2) }

                       { plane stress coefficients }

   G  = E/(1-nu**2)

   C11 = G

   C12 = G*nu

   C22 = G

   C33 = G*(1-nu)/2

   amplitude=1e-6      { a guess for grid-plot scaling }

   mag=0.1/amplitude

   force = 250         { total loading force in Newtons (~10 pound force) }

   dist = 0.5*force*(hW^2-y^2)/I      { Distributed load }

   Sx = [C11*dx(U) + C12*dy(V)]        { Stresses }

   Sy = [C12*dx(U) + C22*dy(V)]

   Txy = C33*[dy(U) + dx(V)]

   Vexact = [force/(6*E*I)]*[(L-x)^2*(2*L+x) + 3*nu*x*y^2]

   Uexact = [force/(6*E*I)]*[3*y*(L^2-x^2) +(2+nu)*y^3 -6*(1+nu)*hW^2*y]

   Sxexact = -force*x*y/I

   Txyexact = -0.5*force*(hW^2-y^2)/I

   small = 1e-5

initial values

   U = 0

   V = 0

equations             { define the displacement equations }

   U:        dx[C11*dx(U) + C12*dy(V)] + dy[C33*(dy(U) + dx(V))] = 0

   V:        dx[C33*(dy(U) + dx(V))]   + dy[C12*dx(U) + C22*dy(V)] = 0

{  force regridder to resolve the shear stress.

    Avoid the ends, where the stress is extreme. }

resolve (Txy, 100*(x/L)*(1-x/L))

boundaries

   region 1

     start (0,-hW)

     load(U)=0         { free boundary on bottom, no normal stress }

     load(V)=0

     line to (L,-hW)

     value(U) = Uexact      { clamp the right end }

     line to (L,0) point value(V) = 0

     line to (L,hW)

     load(U)=0         { free boundary on top, no normal stress }

     load(V)=0

     line to (0,hW)

     load(U) = 0

     load(V) = dist    { apply distributed load to Y-displacement equation }

     line to close

plots

   grid(x+mag*U,y+mag*V)   as "deformation"   { show final deformed grid }

   elevation(V,Vexact) from(0,0) to (L,0) as "Center Y-Displacement(M)"

   elevation(V,Vexact) from(0,hW) to (L,hW) as "Top Y-Displacement(M)"

   elevation(U,Uexact) from(0,hW) to (L,hW) as "Top X-Displacement(M)"

   elevation(Sx,Sxexact) from(0,hW) to (L,hW) as "Top X-Stress"

   elevation(Sx,Sxexact) from(0,0) to (L,0) as "Center X-Stress"

   elevation(Txy,Txyexact) from(0,hW) to (L,hW) as "Top Shear Stress"

   elevation(Txy,Txyexact) from(0,0) to (L,0) as "Center Shear Stress"

   elevation(Txy,Txyexact) from(hL,-hW) to (hL,hW) as "Center Shear Stress"

end