TALK=f;RUN(1,1)
DISPLAY
This is a 2D (xy) Plane-Strain problem, in which a beam is bent by
a torqe applied at its right-hand end.
The problem has an exact analytical solution which is used to enable
the percentage errors on the x- and y-direction displacements to be
computed and printed.
A menu is provided which enables:
* the standard and new (simultaneous) solvers to be compared
* the Poisson's ratio to be varied
* the z-dimension to be varied
* the beam to be made much more slender
* the end torque to be replaced by a tension force
* point-wise and distributed forces to be applied
* temperature gradients to be applied in x and y directions
* the grid-fineness and distribution to be varied.
ENDDIS
PHOTON USE
p;;;;
set prop off
cl
msg x-direction-displacement field
gr ou z 1
cont U1 z 1 x 2 22 y 1 22 fil;.0001
vec z 1 x 2 22 y 2 22 col 0
pause
cl
msg y-direction-displacement field
gr ou z 1
cont V1 z 1 x 2 22 y 1 22 fil;.0001
vec z 1 x 2 22 y 2 22 col 0
ENDUSE
prt0begin
print0 1.0e-20
prt0end
************************************************************
Group 1. Run Title and Number
************************************************************
Declarations and settings
REAL(FX,LX,LY,POISSON,YOUNG)
FX=40.0e6 ! H/m^2 = 40 N/mm^2 STRX(Right) = FX/(LY/2) * (Y - LY/2) = Sigma0 * (Y - LY/2)
LX=120.e-3
LY=40.e-3
FX=FX/(LY/2) ! force gradient
YOUNG = 1/0.5E-11 ! Young's modulus
YOUNG
POISSON=0.3 ! Poisson's ratio
INTEGER(CASENO) ! start of menu
mesg(caseno 1 or 2 : default settings
mesg(caseno 3 or 4 : Poisson's ratio = 0.0
mesg(caseno 5 or 6 : zwlast = 1.0 m
mesg(caseno 7 or 8 : lx = 10.0 * ly
mesg(caseno 9 or 10 : not bending but tension
mesg(caseno 11 or 12 : point downward load at right-hand end
mesg(caseno 13 or 14 : point load but no rotation at RH end
mesg(caseno 15 or 16 : uniform downward load; RH end fixed
mesg(caseno 17 or 18 : no load; x-direction temperature gradient
mesg(caseno 19 or 20 : no load; y-direction temperature gradient
mesg(caseno odd = StressModel, even = ModModel361
caseno=1
mesgm(caseno = :caseno: Enter another if not OK
readvdu(caseno, int, 1)
caseno
libref=caseno ! useful because this appears in graphical monitor
Integer(SolvMod)
! SolvMod = 1 for Stone Solver (caseno = even)
! 2 for simultaneous Solver (caseno = odd )
!
solvmod=1 + (caseno+1)/2-caseno/2
solvmod
TEXT(bent beam in plane strain; s220
************************************************************
Group 2. Time dependence
STEADY = T
************************************************************
INTEGER(NXBODY,NYBODY) ! grid-fineness menu
NXBODY = 21
NYBODY = 21
REAL(POWERX,POWERY)
mesgm(in beam, nx and ny are :nxbody: and :nybody: OK? Y/n
readvdu(ans,char,y)
if(:ans:.eq.n) then
mesgm(insert other nxbody
readvdu(nxbody,int,:nxbody:)
nxbody
mesgm(insert other nybody
readvdu(nybody,int,:nybody:)
nybody
endif
mesgm(grid intervals are uniform. OK? Y/n
readvdu(ans,char,y)
if(:ans:.eq.n) then
mesga(To use power-law x-grid insert powerx
mesg(powerx > 1 gives close spacing at low x
mesg(powerx < 1 gives close spacing at high x
mesg(insert powerx
readvdu(powerx,real,1.0)
powerx
mesga(To use symmetric power-law y-grid insert powery
mesg(powery > 1 gives close spacing at low and high y
mesg(powery < 1 gives close spacing at intermediate y
mesg(insert powery
readvdu(powery,real,1.0)
powerx
endif
Group 3. X-Direction Grid Spacing
if(caseno.eq.7) then
lx=ly*10.0
endif
if(caseno.eq.8) then
lx=ly*10.0
endif
lx
CARTES = T
NREGX=3 ! 3 regions
IREGX=1;GRDPWR(X,1,0.01*LX,1.0) ! single inner fluid cell
IREGX=2;GRDPWR(X,NXBODY,LX,POWERX)
IREGX=3;GRDPWR(X,1,0.01*LX,1.0) ! single outer fluid cell
************************************************************
Group 4. Y-Direction Grid Spacing
NREGY=3 ! 3 regions
IREGY=1;GRDPWR(Y,1,0.01*LY,1.0) ! single inner fluid cell
IREGY=2;GRDPWR(Y,-NYBODY,LY,POWERY)
IREGY=3;GRDPWR(Y,1,0.01*LY,1.0) ! single outer fluid cell
************************************************************
Group 5. Z-Direction Grid Spacing
NZ=1
ZWLAST = 0.001
if(caseno.eq.5.or.caseno.eq.6) then
ZWLAST = 1.0
endif
zwlast
************************************************************
Group 7. Variables: STOREd,SOLVEd,NAMEd
ONEPHS = T
SOLVE(P1,V1,U1)
STORE(PRPS)
STORE(STRX,STRY,STRZ,STXY)
STORE(EPSY,EPSX,EPSZ)
STORE(U1T,V1T,UERR,VERR,DILT)
IF(CASENO.GT.16.AND.CASENO.LT.21) THEN
SOLVE(TEM1)
STORE(EPST,DVO1)
ENDIF
************************************************************
GROUP 8. ITERATION NUMBERS ETC
RESFAC=1.e-8
RESREF(V1)=-1
RESREF(U1)=-1
RESREF(P1)=-1
IF(SOLVMOD.EQ.1) THEN ! combined with the LSWEEP settings below
LITER(V1) = 50 ! these set up a 'fair' contest between the
LITER(U1) = 50 ! two solvers
ELSE
LITER(V1) = 50
LITER(U1) = 50
ENDIF
LITER(P1) = 4
************************************************************
GROUP 9. PROPERTIES
CSG10='Q1' ! materials with differing POISSON ratios
MATFLG=T;NMAT=2
160 7800.0 0.3 473.0 43.0 1.0e-5 0.5e-11
161 7800.0 0.0 473.0 43.0 1.0e-5 0.5e-11
INTEGER(IPRPS)
IPRPS=160 ! used in initial-condition setting below
IF(CASENO.EQ.3.OR.CASENO.EQ.4) THEN
IPRPS=161
POISSON=0.3 ! for use in exact-solution calculation
ENDIF
************************************************************
GROUP 11. INITIAL VALUES
fiinit(p1)=0.0
fiinit(u1)=0.0
fiinit(v1)=0.0
fiinit(U1T) = 0.0
fiinit(V1T)= 0.0
fiinit(UERR)= 0.0
fiinit(VERR)= 0.0
fiinit(DILT)= 0.0
FIINIT(PRPS)=0
PATCH(BODY,INIVAL,2,NX-1,2,NY-1,1,1,1,1)
INIT(BODY,PRPS,FIXVAL,IPRPS)
************************************************************
GROUP 13. BOUNDARY & SPECIAL SOURCES
IF(CASENO.EQ.17.OR.CASENO.EQ.18) THEN ! x-wise temperature gradient
PATCH(LEFTT,CELL,2,2,2,NY-1,1,1,1,1)
COVAL(LEFTT,TEM1,FIXVAL,-10)
PATCH(RIGHTT,CELL,NX-1,NX-1,2,NY-1,1,1,1,1)
COVAL(RIGHTT,TEM1,FIXVAL,10)
PATCH(RIGHT,NORTH,NX-1,NX-1,1,1,1,1,1,1)
COVAL(RIGHT,V1,FIXVAL,0.0)
TEXT(x-wise temperature gradient; s220
fx=0
ENDIF
IF(CASENO.EQ.19.OR.CASENO.EQ.20) THEN ! y-wise temperature gradient
PATCH(bottomT,CELL,2,nx-1,2,2,1,1,1,1)
COVAL(bottomT,TEM1,FIXVAL,-10)
PATCH(topT,CELL,2,NX-1,NY-1,ny-1,1,1,1,1)
COVAL(topT,TEM1,FIXVAL,10)
PATCH(RIGHT,NORTH,NX-1,NX-1,1,1,1,1,1,1)
COVAL(RIGHT,V1,FIXVAL,0.0)
TEXT(y-wise temperature gradient; s220
fx=0
ENDIF
PATCH(AXESUU,WEST,1,1,2,NY-1,1,1,1,1) ! LEFT : x-fixed
COVAL(AXESUU,U1,1.e6,0.0)
char(form)
form=:FX:*(YG-0.51*:LY:) ! 0.51 because thickness is LY+2*0.01*LY
PATCH(FORC01,EAST,NX-1,NX-1,2,NY-1,1,1,1,1) ! Right hand end
if(caseno.eq.9.or.caseno.eq.10) then
form=:fx:*:LY: ! simple tension
TEXT(beam in simple tension; s220
endif
if(caseno.gt.10.and.caseno.lt.15) then
(SOURCE of V1 at FORC01 is COVAL(FIXFLU,-0.1*:fx:)) ! downward force
form=0.0
TEXT(cantilevered beam with point load; s220
ENDIF
if(caseno.eq.15.or.caseno.eq.16) then
patch(topforce,north,2,nx-1,ny-1,ny-1,1,1,1,1)
coval(topforce,v1,fixflu,-0.1*:fx:/:ly:) ! uniform downward force
patch(right,north,nx-1,nx-1,1,1,1,1,1,1)
coval(right,v1,fixval,0.0)
endif
if(caseno.gt.12.and.caseno.lt.17) then
(SOURCE of U1 at FORC01 is COVAL(FIXVAL,0.0)) ! U1 fixed to zero
else
(SOURCE of U1 at FORC01 is COVAL(FIXFLU,form)) ! bending or tension
endif
PATCH(AXESVV,cell,2,2,NY/2+1,NY/2+1,1,1,1,1) ! V1 - fixed
PATCH(AXESVV,cell,2,2,1,1,1,1,1,1) ! V1 fixed at bottom on left
COVAL(AXESVV,V1,FIXVAL,0.0)
! PLANE-STRAIN, EPSZ = 0
SPEDAT(BOUNDARY,ZCONST,R,1.0e20)
if(caseno.eq.13.or.caseno.eq.14) then
TEXT(cantilever; no end rotation: s220
endif
if(caseno.eq.15.or.caseno.eq.16) then
TEXT(uniform load; no end rotation: s220
endif
title
************************************************************
GROUP 15. TERMINATE SWEEPS
IF(SOLVMOD.EQ.1) THEN
LG(40) = F
spedat(rlxfac,rlxu1d,r,0.5)
spedat(rlxfac,rlxv1d,r,0.5)
LSWEEP = 1000
IF(CASENO.EQ.8) THEN
LSWEEP=2000
ENDIF
ELSE
LG(40) = T
LSWEEP = 400
IF(CASENO.EQ.7) THEN
LSWEEP= 600 ! 6
ENDIF
ENDIF
lsweep
ISG21=LSWEEP
************************************************************
GROUP 17. RELAXATION
#CONPROM
RELAX(P1 ,LINRLX, 1.000000E+00)
spedat(rlxfac,rlxu1d,r,1)
spedat(rlxfac,rlxv1d,r,1)
************************************************************
GROUP 19. DATA TRANSMITTED TO GROUND
STRA = T
PARSOL = F
ISG52 = 3 ! probe & res
************************************************************
GROUP 23.FIELD PRINT-OUT & PLOT CONTROL
TSTSWP = - 1 ! graphic-mode
NXPRIN=1; IXPRL=NX-1; IXPRF=NX-5
NYPRIN=1; IYPRL=NY-1; IYPRF=NY-5
IXMON = NX-2
IYMON = ny/2
IZMON = 1
inform7begin
real(RA,RB,RC,RD)
RA = FX/YOUNG*(1-POISSON**2)
RB = -FX/YOUNG/2*(1-POISSON**2)
RC = POISSON/(1-POISSON)
RD = FX/YOUNG*(1+POISSON)*(1-2*POISSON)
RA
RB
RC
RD
**** CALCULATE analytical solution ***
char(FR1,FR2,FR3,FR4)
FR4=(YV-0.51*:LY:)
FR1=:RB:*((YV-0.51*:LY:)^2*:RC:
FR2= +(XG-0.01*:LX:)^2)
FR3= :RA:*(XU-0.01*:LX:)*(YG-0.51*:LY:)
if(caseno.eq.9.or.caseno.eq.10) then
fr3=:RA:*(XU-0.01*:LX:)*:LY:
fr1=(:rb:*:ly:*(YV-0.01*:ly:)*:rc: ! this appears to be incorrect ( correct VIA)
fr2=*2.0)
fr4=:ly:
endif
Exact solutions not yet provided for cases 11 and above
(STORED VAR V1T IS :FR1::FR2: with imat>100)
(STORED VAR DILT IS :RD:*:FR4: with imat>100)
(STORED VAR U1T IS :FR3: with imat>100)
(STORED VAR UERR IS ABS(U1-U1T)/(ABS(U1T)+1.e-10)*100 with imat>100)
(STORED VAR VERR IS ABS(V1-V1T)/(ABS(V1T)+1.e-10)*100 with imat>100)
inform7end
STOP