PHOTON USE
     p



     gr ou z 1
     con thex z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg Linear thermal expansion coefficient
     pause
     con cl;red
     con stif z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg Stiffnes ( Young's module ) contours
     pause
     con cl;red
     vec z 1 sh
     msg Velocity and displacement vectors
     pause
     vec cl;red
     con u1 z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg X-component displacement contours
     pause
     con cl;red
     con v1 z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg Y-component displacement contours
     pause
     con cl;red
     con dil z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg Dilatation contours
     pause
     con cl;red
     con epsx z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg X-strain contours
     pause
     con cl;red
     con epsy z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg Y-strain contours
     pause
     con cl;red
     con epst z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg Thermal strain contours
     pause
     con cl;red
     con strx z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg X-normal stress contours
     pause
     con cl;red
     con stry z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg Y-normal stress contours
     pause
     con cl;red
     con taxy z 1 fil;.001
     gr ou z 1 x 1 3 y 6 m
     gr ou z 1 x 8 m y 6 m
     msg Shear stress contours
     pause
     msg Press E to end
     enduse
    GROUP 1. Run title and other preliminaries

TEXT(2D Stress in solid, XY plane

  DISPLAY
       CASE STUDY: Thermal stresses in composite solid block
  ENDDIS

REAL(LENGTH,WIDTH,ALFA,TISO,STRNY,STRNX,STIFFN,STRSX,STRSY)
REAL(NORTEM,SOUTEM,HEATSOR,POISSN)
BOOLEAN(ISOTHER,TEMLIN,TEMSOR)

TISO=100.0
LENGTH=1.0;WIDTH =1.0

STIFFN=2.e11
STRSY =-1.e8;STRSX=0.0
ALFA  =1.e-05          ; RG(1)=ALFA
STRNY =STRSY           ; RG(2)=STRNY
STRNX =STRSX           ; RG(3)=STRNX
POISSN=0.3333          ; RG(4)=POISSN

NX  =10 ; NY  =10

ISOTHER=t;TEMLIN=f;TEMSOR=f


    GROUP 3. X-direction grid specification
GRDPWR(X,NX,LENGTH,1.0)
    GROUP 4. Y-direction grid specification
GRDPWR(Y,NY,WIDTH,1.0)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1,TEM)
      Provide the storage for thermal expansion coefficient
      ,THEX, stiffnes, STIF, and Poisson ratio to be initialised
      for the presence of different materials

STORE(THEX,STIF,POIS)

STORE(CON1,CON2,CON3,CON4)
    GROUP 8. Terms (in differential equations) & devices
TERMS(U1,Y,N,Y,Y,Y,Y)
TERMS(V1,Y,N,Y,Y,Y,Y)
terms(tem,n,n,y,p,p,p)
    GROUP 9. Properties of the medium (or media)
ENUL=1.
    GROUP 11. Initialization of variable or porosity fields
IF(ISOTHER) THEN
iniadd=f
PATCH(STELP,INIVAL,1,NX,1,NY,1,1,1,1)
INIT (STELP,POIS,0.0,POISSN)

PATCH(STELT,INIVAL,1,NX,1,NY,1,1,1,1)
INIT (STELT,THEX,0.0,ALFA)

PATCH(STELS,INIVAL,1,NX,1,NY,1,1,1,1)
INIT (STELS,STIF,0.0,STIFFN)

fiinit(tem)=TISO
   #
PATCH(ALMT1,INIVAL,1,3,NY/2+1,NY,1,1,1,1)
INIT (ALMT1,THEX,0.0,10.*ALFA)
INIT (ALMT1,STIF,0.0,0.1*STIFFN)
PATCH(ALMT2,INIVAL,8,NX,NY/2+1,NY,1,1,1,1)
INIT (ALMT2,THEX,0.0,10.*ALFA)
INIT (ALMT2,STIF,0.0,0.1*STIFFN)
ENDIF
    GROUP 13. Boundary conditions and special sources
    ** Linear temperature field
IF(TEMLIN) THEN
PATCH(TEMPN,NWALL,1,NX,NY,NY,1,1,1,1)
COVAL(TEMPN,TEM,1.0,0.0)
PATCH(TEMPS,SWALL,1,NX,1,1,1,1,1,1)
COVAL(TEMPS,TEM,1.0,1.0)
ENDIF
    ** Source generated temperature field
IF(TEMSOR) THEN
PATCH(TEMPN,NWALL,1,NX,NY,NY,1,1,1,1)
COVAL(TEMPN,TEM,1.0,0.0)
PATCH(TEMSOR,VOLUME,1,NX,1,NY,1,1,1,1)
COVAL(TEMSOR,TEM,FIXFLU,2.0)
ENDIF

    PLANTBEGIN
PATCH(STELC,INIVAL,1,NX,1,NY,1,1,1,1)
    VAL=(1.-RG(4))/((1.+RG(4))*(1-2.*RG(4)))
INIT (STELC,CON1,0.0,GRND)
    VAL=RG(4)/((1.+RG(4))*(1-2.*RG(4)))
INIT (STELC,CON2,0.0,GRND)
    VAL=1./(1.-2.*RG(4))
INIT (STELC,CON3,0.0,GRND)
    VAL=1./(2.*(1+RG(4)))
INIT (STELC,CON4,0.0,GRND)
              Stress - in -Solid BC
    ---------------- WEST ----------------------------
   ** Zero displacement at the west boundary
PATCH(WESZD,WWALL,1,1,1,NY,1,1,1,1)
COVAL(WESZD,U1,1.,0.0)
    ---------------- EAST ----------------------------
   ** Normal stress at the east solid boundary
PATCH(EASZN,EAST,NX-1,NX-1,1,NY,1,1,1,1)
    VAL=U1[-1,,]$
   +0.333*(P1+4.*(THEX*TEM+RG(3)/STIF))*DXU2D
COVAL(EASZN,U1,FIXVAL,GRND)
   ** Zero shear stress at the east boundary
PATCH(EASZS,EAST,NX,NX,1,NY-1,1,1,1,1)
    VAL=-(U1[-1,+1,]-U1[-1,,])/DYG2D
COVAL(EASZS,V1,FIXFLU,GRND)
   ** Dilatation at the east solid boundary
PATCH(FXDEAS,EAST,NX,NX,2,NY-1,1,1,1,1)
    VAL=2.*(THEX*TEM-RG(3)/STIF)$
                -1.5*(V1-V1[,-1,])/DYV2D
COVAL(FXDEAS,P1,FIXVAL, grnd)

PATCH(FXDEAC,EAST,NX,NX,1,1,1,1,1,1)
    VAL=2.*(THEX*TEM-RG(3)/STIF)$
                -1.5*(V1[,+1,]-V1)/DYV2D[,+1,]
COVAL(FXDEAC,P1,FIXVAL, grnd)

    ---------------- NORTH ----------------------------
   ** Normal stress at the free north boundary
PATCH(NORZN,CELL,1,NX,NY-1,NY-1,1,1,1,1)
    VAL=V1[,-1,]$
    +0.333*(P1+4.*(THEX*TEM+RG(2)/STIF))*DYV2D
COVAL(NORZN,V1,FIXVAL,grnd)
   ** Zero shear stress at the free north boundary
PATCH(NORZS,NORTH,1,NX-1,NY,NY,1,1,1,1)
    VAL=-(V1[+1,-1,]-V1[,-1,])/DXG2D
COVAL(NORZS,U1,FIXFLU,GRND)
   ** Dilatation at the free north boundary
PATCH(FDNW,NORTH,1,1,NY,NY,1,1,1,1)
    VAL=2.*(THEX*TEM-RG(2)/STIF)$
                -1.5*(U1[+1,,]-U1)/DXU2D[+1,,]
COVAL(FDNW,P1,FIXVAL, grnd)
PATCH(FDNOR,NORTH,2,NX-1,NY,NY,1,1,1,1)
    VAL=2.*(THEX*TEM-RG(2)/STIF)$
                -1.5*(U1-U1[-1,,])/DXU2D
COVAL(FDNOR,P1,FIXVAL, grnd)
PATCH(FDNEC,NORTH,NX,NX,NY,NY,1,1,1,1)
    VAL=2.*(THEX*TEM-RG(2)/STIF)$
                -1.5*(U1[-1,,]-U1[-2,,])/DXU2D[-1,,]
COVAL(FDNEC,P1,FIXVAL, grnd)
    ---------------- SOUTH ----------------------------
   ** Zero displacement at the south symmetry
PATCH(SOUZD,SWALL,1,NX,1,1,1,1,1,1)
COVAL(SOUZD,V1,1.,0.0)
    ==================================================
   ** Volumetric mass source for dilatation
PATCH(DILAT,VOLUME,1,nx,1,ny,1,1,1,1)
    VAL=-P1/CON3+THEX*TEM/CON4
COVAL(DILAT,P1,FIXFLU,GRND)
    ===================================================
            OUTPUT INFORMATION
   ** Dilatation
STORE(DIL)
    DIL=-P1/CON3+THEX*TEM/CON4
   ** Stress and strain calculations
STORE(STRX,STRY,TAXY)
STORE(EPSX,EPSY,EPST)
      * Thermal strains
   EPST=THEX*TEM
  REGION(1,NX,1,NY,1,1)
      * Normal strains
   EPSX=(U1-WEST(U1))/DXU2D
  REGION(2,NX-1,1,NY,1,1)
   EPSX=(U1[+1,,]-U1)/DXU2D[+1,,]
  REGION(1,1,1,NY,1,1)
   EPSX=(U1[-1,,]-U1[-2,,])/DXU2D[-1,,]
  REGION(NX,NX,1,NY,1,1)
   EPSY=(V1-SOUTH(V1))/DYV2D
  REGION(1,NX,2,NY-1,1,1)
   EPSY=(V1[,+1,]-V1)/DYV2D[,+1,]
  REGION(1,NX,1,1,1,1)
   EPSY=(V1[,-1,]-V1[,-2,])/DYV2D[,-1,]
  REGION(1,NX,NY,NY,1,1)
      * Normal stresses
   STRX=STIF*(CON1*EPSX+CON2*EPSY-CON3*EPST)
  REGION(1,NX,1,NY,1,1)
   STRY=STIF*(CON1*EPSY+CON2*EPSX-CON3*EPST)
  REGION(1,NX,1,NY,1,1)
      * Shear stresses
STORE(DDYU,DDXV)
              == ddy.u ==
   DDYU=0.25                          * $
    ((U1        -U1[,-1,]  )/DYG2D[,-1,]      + $
     (U1[,+1,]  -U1        )/DYG2D            + $
     (U1[-1,,]  -U1[-1,-1,])/DYG2D[,-1,]    + $
     (U1[-1,+1,]-U1[-1,,]  )/DYG2D )
  REGION(2,NX-1,2,NY-1,1,1)
   DDYU=0.5                            * $
           ((U1-U1[,-1,])/DYG2D[,-1,]          + $
            (U1[,+1,]-U1)/DYG2D)
  REGION(1,1,2,NY-1,1,1)
   DDYU=(U1[,+1,]-U1)/DYG2D
  REGION(1,1,1,1,1,1)
   DDYU=(U1-U1[,-1,])/DYG2D[,-1,]
  REGION(1,1,NY,NY,1,1)
   DDYU=0.5                            * $
           ((U1[-1,,]-U1[-1,-1,])/DYG2D[,-1,]+ $
            (U1[-1,+1,]-U1[-1,,])/DYG2D )
  REGION(NX,NX,2,NY-1,1,1)
   DDYU=(U1[-1,+1,]-U1[-1,,])/DYG2D
  REGION(NX,NX,1,1,1,1)
   DDYU=(U1[-1,,]-U1[-1,-1,])/DYG2D[,-1,]
  REGION(NX,NX,NY,NY,1,1)
   DDYU=0.5                            * $
           ((U1-U1[,-1,])/DYG2D[,-1,]          + $
            (U1[-1,,]-U1[-1,-1,])/DYG2D[,-1,])
  REGION(2,NX-1,NY,NY,1,1)
   DDYU=0.5                            * $
           ((U1[,+1,]-U1)/DYG2D                + $
            (U1[-1,+1,]-U1[-1,,])/DYG2D)
  REGION(2,NX-1,1,1,1,1)
              == ddx.v ==
   DDXV=0.25                          * $
    ((V1[+1,-1,]-V1[,-1,]  )/DXG2D            + $
     (V1[+1,,]  -V1        )/DXG2D            + $
     (V1[,-1,]  -V1[-1,-1,])/DXG2D[-1,,]      + $
     (V1        -V1[-1,,]  )/DXG2D[-1,,] )
  REGION(2,NX-1,2,NY-1,1,1)
   DDXV=(V1[+1,,] -V1)/DXG2D
  REGION(1,1,1,1,1,1)
   DDXV=0.5                           * $
    ((V1[+1,-1,]-V1[,-1,]  )/DXG2D            + $
     (V1[+1,,]  -V1        )/DXG2D )
  REGION(1,1,2,NY-1,1,1)
   DDXV=(V1[+1,-1,]-V1[,-1,]  )/DXG2D
  REGION(1,1,NY,NY,1,1)
   DDXV=(V1-V1[-1,,])/DXG2D[-1,,]
  REGION(NX,NX,1,1,1,1)
   DDXV=0.5                           * $
    ((V1[,-1,]  -V1[-1,-1,])/DXG2D[-1,,]      + $
     (V1        -V1[-1,,]  )/DXG2D[-1,,] )
  REGION(NX,NX,2,NY-1,1,1)
   DDXV=(V1[,-1,] -V1[-1,-1,])/DXG2D[-1,,]
  REGION(NX,NX,NY,NY,1,1)
   DDXV=0.5                           * $
    ((V1[+1,,]  -V1        )/DXG2D            + $
     (V1        -V1[-1,,]  )/DXG2D[-1,,] )
  REGION(2,NX-1,1,1,1,1)
   DDXV=0.5                           * $
    ((V1[+1,-1,]-V1[,-1,]  )/DXG2D            + $
     (V1[,-1,]  -V1[-1,-1,])/DXG2D[-1,,])
  REGION(2,NX-1,NY,NY,1,1)

   TAXY=STIF*CON4*(DDYU+DDXV)
  REGION(1,NX,1,NY,1,1)
    PLANTEND


RELAX(P1,LINRLX,0.005)
RELAX(U1,FALSDT,0.0005)
RELAX(V1,FALSDT,0.0005)

NAMSAT=MOSG

    GROUP 15. Termination of sweeps
LSWEEP=3000
selref=f
RESREF(P1)=1.E-8
RESREF(U1)=1.E-8;RESREF(V1)=1.E-8
tstswp=-1
    GROUP 22. Spot-value print-out
IYMON=NY-1
    GROUP 23. Field print-out and plot control
NXPRIN=1;NYPRIN=1
 LIBREF=458

tstswp=-1
dmpstk=t
DISTIL=T
 EX(P1)=4.671E-03; EX(U1)=3.079E-03; EX(V1)=6.590E-04
 EX(DDXV)=2.347E-03; EX(DDYU)=3.255E-03; EX(EPST)=3.700E-03
 EX(EPSY)=3.955E-03; EX(EPSX)=6.382E-03; EX(TAXY)=1.362E+08
 EX(STRY)=3.664E+08; EX(STRX)=1.830E+08; EX(DIL)=8.618E-03
 EX(CON4)=3.750E-01; EX(CON3)=2.999E+00; EX(CON2)=7.498E-01
 EX(CON1)=1.500E+00; EX(POIS)=3.333E-01; EX(STIF)=1.460E+11
 EX(THEX)=3.700E-05; EX(TEM)=1.000E+02
 LIBREF=458
STOP