{ TENSION.PDE }

 { *******************************************************
   This example shows the deformation of a tension bar with a hole.

   The equations of Stress/Strain arise from the balance of forces in a
   material medium, expressed as

        dx(Sx) + dy(Txy) + Fx = 0
        dx(Txy) + dy(Sy) + Fy = 0

   where Sx and Sy are the stresses in the x- and y- directions,
   Txy is the shear stress, and Fx and Fy are the body forces in the
   x- and y- directions.

   The deformation of the material is described by the displacements,
   U and V, from which the strains are defined as

        ex = dx(U)
        ey = dy(V)
        gxy = dy(U) + dx(V).

   The eight quantities U,V,ex,ey,gxy,Sx,Sy and Txy are related through the
   constitutive relations of the material. In general,

        Sx  =  C11*ex + C12*ey + C13*gxy - b*Temp
        Sy  =  C12*ex + C22*ey + C23*gxy - b*Temp
        Txy =  C13*ex + C23*ey + C33*gxy

   In orthotropic solids, we may take C13 = C23 = 0.

   Combining all these relations, we get the displacement equations:

        dx[C11*dx(U)+C12*dy(V)] + dy[C33*(dy(U)+dx(V))] + Fx = dx(b*Temp)
        dy[C12*dx(U)+C22*dy(V)] + dx[C33*(dy(U)+dx(V))] + Fy = dy(b*Temp)

   In the "Plane-Stress" approximation, appropriate for a flat, thin plate
   loaded by surface tractions and body forces IN THE PLANE of the plate,
   we may write

        C11 = G         C12 = G*nu
                        C22 = G
                                        C33 = G*(1-nu)/2

   where G = E/(1-nu**2)
         E is Young's Modulus
   and   nu is Poisson's Ratio.

   The displacement form of the stress equations (for uniform temperature
   and no body forces) is then (dividing out G):

        dx[dx(U)+nu*dy(V)] + 0.5*(1-nu)*dy[dy(U)+dx(V)] = 0
        dy[nu*dx(U)+dy(V)] + 0.5*(1-nu)*dx[dy(U)+dx(V)] = 0

   In order to quantify the load boundary condition mechanism,
   consider the stress equations in their original form:

        dx(Sx) + dy(Txy) = 0
        dx(Txy) + dy(Sy) = 0

   These can be written as

        div(P) = 0
        div(Q) = 0

   where P = [Sx,Txy]
   and   Q = [Txy,Sy]

   The "load" (or "natural") boundary condition for the U-equation defines
   the outward surface-normal component of P, while the load boundary condition
   for the V-equation defines the surface-normal component of Q. Thus, the
   load boundary conditions for the U- and V- equations together define
   the surface load vector.

   On a free boundary, both of these vectors are zero, so a free boundary
   is simply specified by

        load(U) = 0
        load(V) = 0.

   Here we consider a tension strip with a hole, subject to an X-load.

  *************************************************************** }

  Title 'Plane Stress tension strip with a hole'

  select
   NODELIMIT=1000
!   NGRID=30
    regrid=off
    errlim = 1e-3	{ increase accuracy to resolve stresses }
 !   painted	{ paint all contour plots }
 !   fixdt!=0.01

  variables

    U(1e-10)                   { declare U and V to be the system variables }
    V(1e-10)
   Vt(1e-10)
   Ut(1e-10)
   Sxt(1e4)
   Syt(1e4)
   Szt(1e4)
   txyt(1e4)
    Ep(1e-10)
   Stresstt(1e4)
!  vx(1e-11)
!  xm(1e-10) = move(x)
!  vy(1e-11)
!  ym(1e-10) = move(y)
 definitions
!   transfer('plastic_5.txt',Ut0,Vt0,Sxt0,Syt0,Szt0,Txyt0,Ep0)
 t0=0
 t2=10
 ddt=0.1
! t2=3.5
 dte=0.1*ddt
 Temp=15
  Ls=5e-9
! slim=60e-9/t2
 slim=0
  Lw=100e-9-60e-9
  Lh=280e-9-60e-9
  Lgh=0
  Lsp=300e-9-Lw


 stress0=0
 stressE=-1.5E9*1.5

eps=1e-10
!
!xp0=EVAL(Ut,-Ls,y )
!yp0=EVAL(y+Vt,-Ls,y )

xp0=Eval(Ut,0,y)
yp0=Eval(y+Vt,0,y)
xp1=Eval(Ut,Lw,y)
yp1=Eval(y+Vt,Lw,y)
xp01=val(Ut,0,Lh-eps)
yp01=val(y+Vt,0,Lh-eps)
xp11=val(Ut,Lw,Lh-eps)
yp11=val(y+Vt,Lw,Lh-eps)


x00=VAL(Ut,0,Lh )
y00=VAL(y+Vt,0,Lh )
rad00=IF(yp0<y00-eps) THEN arctan((x00-xp0)/(y00-yp0)) ELSE arctan((x00-xp01)/(y00-yp01))



x10=VAL(Ut,Lw,Lh )
y10=VAL(y+Vt,Lw,Lh )

 rad10=IF(yp1<y10-eps) THEN arctan((x10-xp1)/(y10-yp1)) ELSE arctan((x10-xp11)/(y10-yp11))



rad0X=IF(yp0<y10) THEN arctan((Lsp+xp0-x10)/(y10-yp0)) ELSE pi/2
rad1X=IF(yp1<y00) THEN arctan((Lsp-xp1+x00)/(y00-yp1)) ELSE  pi/2


 rad1Y=SAVE(Eval((rad1X-rad10*ustep(rad10))/pi,Lw+Ls,y))
 rad1Z=SAVE(Eval(-rad10*ustep(-rad10)/pi,Lw+Ls,y))

 rad0Y=SAVE(Eval((rad0X+rad00*ustep(-rad00))/pi,-Ls,y))
 rad0Z=SAVE(Eval(rad00*ustep(rad00)/pi,-Ls,y))




   ratio=0.5
  ratioL=save((1+ratio)*rad0Z+(1-ratio)*rad0Y)
  ratioR=save((1-ratio)*rad1Z+(1+ratio)*rad1Y)
slimL=SAVE(slim*ratioL)
slimR=SAVE(slim*ratioR)

stressL= IF(Eval(Stresstt,-Ls,y)<StressE) THEN    0 else StressE*ratioL
stressR=IF(Eval(Stresstt,Lw+Ls,y)<StressE) THEN     0 else StressE*ratioR
stressT=IF(t>1) THEN     0 else StressE




 StressV=0




 sv=34e6
 st=10E6

 Stress=stressV



   H=st
   S0=sv
    nu = 0.34            { define Poisson's Ratio }
    E  = 2.7E9           { Young's Modulus x 10**-11 }
    k=0.18
!   Tp=50

!   Tp=50

   G=0.5*E/(1+nu)
    C11 = 2*G*(1-nu)/(1-2*nu)
    C12 = 2*G*nu/(1-2*nu)
     C14=0
    C22 = 2*G*(1-nu)/(1-2*nu)
    C24=0
    C44 = G
    C13=2*G*nu/(1-2*nu)
    C23=2*G*nu/(1-2*nu)
    C43=0
   Ex=dx(U)
   Ey=dy(V)
   Gxy=dy(U) + dx(V)

       Sx1  =  C11*Ex+ C12*Ey + C14*Gxy-Stress
       Sy1  =  C12*Ex + C22*Ey + C24*Gxy-Stress
       Sz1 =  C13*Ex + C23*Ey + C43*Gxy-Stress
       Txy1 =  C14*Ex+ C24*Ey  + C44*Gxy

       Sxt1=Sxt+Sx1*dte
       Syt1=Syt+Sy1*dte
       Szt1=Szt+Sz1*dte
       Txyt1=Txyt+Txy1*dte
    SR=sqrt(0.5*(Sxt-Syt)^2+0.5*(Syt-Szt)^2+0.5*(Szt-Sxt)^2+3*Txyt^2)
    S1=sqrt(0.5*(Sxt1-Syt1)^2+0.5*(Syt1-Szt1)^2+0.5*(Szt1-Sxt1)^2+3*Txyt1^2)


    S=max(S1,1e-5)

    p=(Sxt+Syt+Szt)/3
    Sxx=Sxt-p
    Syy=Syt-p
    Szz=Szt-p

    Hep=S0+H*Ep
     eps0=3*G/S/(H+3*G)
     dEp=(Sxx*Ex+Syy*Ey+Txyt*Gxy)*eps0
     code=1-ustep(S-S0)*ustep(dEp-1E-6)
     E0=(1-code)*9*G^2/S^2/(H+3*G)
  !  E0=0
    D11=SAVE(C11-E0*Sxx^2)
    D12=SAVE(C12-E0*Sxx*Syy)
    D14=SAVE(C14-E0*Sxx*Txyt)
    D22=SAVE(C22-E0*Syy^2)
    D24=SAVE(C24-E0*Syy*Txyt)
    D44=SAVE(C44-E0*Txyt^2)
    D13=SAVE(C13-E0*Sxx*Szz)
    D23=SAVE(C23-E0*Syy*Szz)
    D43=SAVE(C43-E0*Txyt*Szz)




        Sx  =  D11*Ex+ D12*Ey + D14*Gxy-Stress
        Sy  =  D12*Ex + D22*Ey + D24*Gxy-Stress
        Sz =  D13*Ex + D23*Ey + D43*Gxy-Stress
        Txy =  D14*Ex+ D24*Ey  + D44*Gxy


  initial values

     U =0              { tell SPDE the approximate variable range }
    V = 0

    Ut=0
    Vt=0
   Sxt=0
   Syt=0
   Szt=0
   Txyt=0
   Ep=0
   Stresstt=0
!  vx = 0
!  vy = 0

   equations
{ the heat equation }
!    Tp:	dx[k*dx(Tp)] + dy[k*dy(Tp)] = 0
{ define the Plane-Stress displacement equations }
  U:      dx(Sx) + dy(Txy) = 0
  V:      dx(Txy) + dy(Sy) = 0
  Ut: dt(Ut)=U
  Vt:  dt(Vt)=V
  Sxt:  dt(Sxt)=Sx
  Syt:  dt(Syt)=Sy
  Szt:  dt(Szt)=Sz
  Txyt:  dt(Txyt)=Txy
  Ep: dt(Ep)=dEp
  Stresstt:dt(Stresstt)=Stress
!vx:  dxx(vx)=0
!Xm:  dt(xm) = vx
!vy:  dyy(vy)=0
!Ym:  dt(ym) = vy

  boundaries
    region 1 {Resist}
      Stress=stressV
      S0=sv

      start  (0,0)
      value(U)=0        { fixed displacement on left edge }
      value(V)=0
!      nobc(vx) nobc(xm) value(vy) = 0 dt(ym) = 0
      line to (Lw,0)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      value(vx) = -slim dt(xm)=-slim nobc(vy) nobc(ym)
      line to (Lw,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
  !     nobc(vx) nobc(xm) value(vy) = -slim dt(ym) = -slim
      line to (0,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
   !    value(vx) = slim dt(xm)=slim nobc(vy) nobc(ym)
      line to close



    region 2 {Left}
      Stress=stressL
      S0=sv



      start (0,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
    !   nobc(vx) nobc(xm) value(vy) = -slim dt(ym) = -slim
      line to (-Ls,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
    !   value(vx) = slim dt(xm)=slim nobc(vy) nobc(ym)
      line to (-Ls,0)
      value(U)=0        { fixed displacement on left edge }
      value(V)=0
     ! nobc(vx) nobc(xm) value(vy) = 0 dt(ym) = 0
      Line to  (0,0)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
     !  value(vx) = slim dt(xm)=slim nobc(vy) nobc(ym)
      line to close

    region 3 {right}

      Stress=stressR
      S0=sv
      start (Lw,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
     !  nobc(vx) nobc(xm) value(vy) = -slim dt(ym) = -slim
      line to (Lw+Ls,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
     !  value(vx) = -slim dt(xm)=-slim nobc(vy) nobc(ym)
      line to (Lw+Ls,0)
      value(U)=0        { fixed displacement on left edge }
      value(V)=0
     ! nobc(vx) nobc(xm) value(vy) = 0 dt(ym) = 0
      line to (Lw,0)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
     !  value(vx) = -slim dt(xm)=-slim nobc(vy) nobc(ym)
      line to close

    region 4 {top}
      Stress=stressT
      start (-2*Ls,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
      line to (-2*Ls,Lh+Ls)
      line to (Lw+2*Ls,Lh+Ls)
      line to (Lw+2*Ls,Lh)
      line to close



    region 5 {Left}
      Stress=0
      S0=sv



      start (-Ls,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      nobc(vx) nobc(xm) value(vy) = -slim dt(ym) = -slim
      line to (-2*Ls,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      value(vx) = Eval(slimL,-Ls,y) dt(xm)=Eval(slimL,-Ls,y)  nobc(vy) nobc(ym)
      line to (-2*Ls,0)
      value(U)=0        { fixed displacement on left edge }
      value(V)=0
 !     nobc(vx) nobc(xm) value(vy) = 0 dt(ym) = 0
      Line to  (-Ls,0)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      value(vx) = Eval(slimL,-Ls,y) dt(xm)=Eval(slimL,-Ls,y)  nobc(vy) nobc(ym)
      line to close

    region 6 {right}

      Stress=0
      S0=sv
      start  (Lw+Ls,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      nobc(vx) nobc(xm) value(vy) = -slim dt(ym) = -slim
      line to (Lw+2*Ls,Lh)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      value(vx) = -Eval(slimR,Lw+Ls,y) dt(xm)=-Eval(slimR,Lw+Ls,y) nobc(vy) nobc(ym)
      line to (Lw+2*Ls,0)
      value(U)=0        { fixed displacement on left edge }
      value(V)=0
 !     nobc(vx) nobc(xm) value(vy) = 0 dt(ym) = 0
      line to (Lw+Ls,0)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      value(vx) = -Eval(slimR,Lw+Ls,y) dt(xm)=-Eval(slimR,Lw+Ls,y) nobc(vy) nobc(ym)
      line to close

    region 7 {top}
      Stress=0
      start (-2*Ls,Lh+Ls)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !     value(vx) = Eval(slimL,-Ls,y) dt(xm)=Eval(slimL,-Ls,y)  nobc(vy) nobc(ym)
!       nobc(vx)  nobc(xm) nobc(vy) nobc(ym)
      line to (-2*Ls,Lh+2*Ls)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      nobc(vx) nobc(xm) value(vy) = -slim dt(ym) = -slim
      line to (Lw+2*Ls,Lh+2*Ls)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      nobc(vx)  nobc(xm) nobc(vy) nobc(ym)
 !       value(vx) = -Eval(slimR,Lw+Ls,y) dt(xm)=-Eval(slimR,Lw+Ls,y) nobc(vy) nobc(ym)
      line to (Lw+2*Ls,Lh+Ls)
      load(U)=0      { free boundary, no normal stress }
      load(V)=0
 !      nobc(vx) nobc(xm) value(vy) = -slim dt(ym) = -slim
      line to close



Feature

START"left" (-Ls,0) LINE TO (-Ls,Lh)
START"right" (Lw+Ls,0) LINE TO (Lw+Ls,Lh)
START"resist" (0,0) LINE TO (0,Lh) Line to (Lw,Lh) Line to (Lw,0)
START"outer" (-Ls,0) LINE TO (-Ls,Lh+Ls) Line to (Lw+Ls,Lh+Ls) Line to (Lw+Ls,0)


time

 from 0 to t2  by ddt


  monitors
  for t=0 {by te*0.01 to te*0.1} by ddt to endtime

    grid(x,y)
     grid(x+U*t2,y+V*t2)   ! zoom(0,0, Lw,Lh)  { show final deformed grid }
     grid(x+Ut,y+Vt)
!     elevation(vx,vy,-ratioR*slim,ratioL*slim)  on 'resist'
     elevation( 0,stresstt, stressE)  on 'left'
     elevation( 0,stresstt, stressE)  on'right'
!     elevation( 0,slimL,slim)  on 'left'
!     elevation( 0,slimR, slim)  on'right'
     elevation(rad10,rad00,rad0X,rad1X)  from (-Ls,0) to (-Ls,Lh)
     elevation(rad0Z,rad0Y,rad1Z,rad1Y)  from (-Ls,0) to (-Ls,Lh)

  plots                 { hardcopy at to close: }
  for t=0 {by te*0.01 to te*0.1} by ddt*5 to endtime
    grid(x,y)
     grid(x+U*t2,y+V*t2)   ! zoom(0,0, Lw,Lh)  { show final deformed grid }
     grid(x+Ut,y+Vt)
     grid(x+Ut,y+Vt)   on region 1
      contour(Stresstt) as "stress" painted
     vector(Ut,Vt) as "Displacement"        on region 1
     contour(Ut) as "X-Displacement" on region 1 painted
     contour(Vt) as "Y-Displacement" on region 1 painted
     contour(Ep) as "Ep" on region 1 painted
     contour(Hep) as "H(ep)" on region 1 painted
     contour(S) as "stress" on region 1 painted
     contour(dEp) as "dEp" on region 1 painted
     contour(code) as "1:elastic,0:plastic" on region 1 painted
     contour(SR-Hep) as "error:should be <0" on region 1 painted
     contour(C11) as "C11" on region 1 painted
     contour(D11) as "D11" on region 1 painted

{
     elevation(x+U,y+V)  from (0,0) to (0,Lh) export format "#x#b#y#b#1#b#2" file="LeftA2.txt"
     elevation(x+U,y+V)  from (0,Lh) to (Lw,Lh) export format "#x#b#y#b#1#b#2" file="TopA2.txt"
     elevation(x+U,y+V)  from (Lw,Lh) to (Lw,0) export format "#x#b#y#b#1#b#2" file="rightA2.txt"
}
!     elevation(x+Ut,y+Vt)  on 'left' export format "#x#b#y#b#1#b#2" file="leftA3.txt"
     elevation(x+Ut,y+Vt)  on 'resist' export format "#x#b#y#b#1#b#2" file="resistI20.txt"
     elevation(x+Ut,y+Vt)  on 'outer' export format "#x#b#y#b#1#b#2" file="outerI20.txt"
!     elevation(x+Ut,y+Vt)  on 'right' export format "#x#b#y#b#1#b#2" file="rightA3.txt"


     elevation( 0,stresstt, stressE)  on 'left'
     elevation( 0,stresstt, stressE)  on'right'
     elevation(rad10,rad00,rad0X,rad1X)  on'left'
     elevation(rad0Z,rad0Y,rad1Z,rad1Y)  on'right'


  transfer(Ut,Vt,Sxt,Syt,Szt,Txyt,Ep,Stresst) file'Is20.txt'

histories
!    history(rate)  range=(0,te)





  end     15788318