{ St2LS_1.pde - program for FlexPDE}
 { Plane stress tension of the plate with a hole}
 { The quarter of the plate is considered by virtue of the symmetry problem}

  Title 'Plane Stress tension strip with a hole'

  select
    errlim = 0.5e-4	{error tolerance for stress determinations}

  variables
    U           {displacements U and V are the variables}
    V

  definitions
    a=5	{length of a plate part}
    b=2	{width of a plate part}
    R=0.25	{radius of a hole}
    p=1e3  {value of the stretching loading (kG/sm^2)}
    nu = 0.3	{Poisson's ratio}
    E  = 2e6 	{Young's modulus  (kG/sm^2)}

    mu  = E/2/(1+nu)	{shear modulus}
    lambda1 =  2*mu*nu/(1-2*nu)   {Lame's coefficient}
	{lambdaz is the Lame's coefficient for plane stress condition}
    lambdaz = 2*mu*lambda1/(lambda1+2*mu)
    Sxx = dx(U)		{axial strain Sxx}
    Syy = dy(V)		{axial strain Syy}
    Sxy = (dy(U)+dx(V))/2	{shear strain Sxy}
    Txx= lambdaz*(Sxx+Syy)+2*mu*Sxx        {stress Txx}
    Tyy= lambdaz*(Sxx+Syy)+2*mu*Syy       {stress Tyy}
    Txy=2*mu*Sxy			 {stress Txy}
    mcoeff=0.5    {scale factor for grid-plot output}
    Scal=mcoeff*globalmax(magnitude(x,y))/globalmax(magnitude(u,v))


{for tangential stresses computation, can be omitted for lab1 individual tasks}
    RR=sqrt(x*x+y*y)     theta=arcsin(y/RR)  {RR,  theta are the polar coordinates}
    {components of the stress tensor in polar coordinate system}
    Trr=Txx*(cos(theta))^2+Tyy*(sin(theta))^2+Txy*sin(2*theta)
    Tthetatheta=Txx*(sin(theta))^2+Tyy*(cos(theta))^2-Txy*sin(2*theta)

  equations      {differential equilibrium equations for a plane elastic problem}
    U:  dx(Txx) +dy(Txy)  = 0
    V:  dx(Txy)+dy(Tyy)  = 0


  boundaries       {region and boundary conditions}
    region 1
      start (R,0)
      load(U)=0       {symmetry conditions along the axis Y=0}
      value(V)=0
      line to (a,0)

      load(U)=p      {pressure}
      load(V) = 0
      line to (a,b)

      load(U)=0        {free boundary, no stress}
      load(V)=0
      line to (0,b)

      value(U)=0     {symmetry conditions along the axis X=0}
      load(V) = 0
      line to (0,R)

      load(U)=0        {free boundary, no stress}
      load(V)=0
      arc(center=0,0)  to close


  monitors
      grid(x+Scal*U,y+Scal*V)    {show deformed grid of solution in progress}

  plots               {resulting plots}
     grid(x+Scal*U,y+Scal*V)  as "Deformed shape"  {show deformed grid}
report(Scal)
    contour(U) painted as "Displacement Ux" {gradient fill of colors}
     contour(U) as "Displacement Ux" {isolines}
      contour(V) painted as "Displacement Uy"
     vector(U,V) as "Displacement vector"
contour(magnitude(u,v)) painted as"Magnitude of displacement vector"
       elevation(U) from (0,(b-R)/2) to (a,(b-R)/2) as "Graph for displacement Ux"
       contour(Txx) painted as "Stress Txx"
      contour(Txy) painted as "Stress Txy"
     contour(Tyy) painted as "Stress Tyy"
      contour(Sxx) painted as "Strain Sxx"
      contour(Sxy) painted as "Strain Sxy"
     contour(Syy) painted as "Strain Syy"

    {Pictures of tangential stresses distribution, for stress T_thetatheta}
     contour(Tthetatheta) painted as "Stress T_thetatheta"
     contour(Tthetatheta)  as "Stress T_thetatheta"
     elevation(Tthetatheta) from (0,R) to (0,b) as "T_thetatheta along line x=0"
  end