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 accuracy of the computed values to be assessed by
being divided by the exact values, and printed as U1/T, etc ,
which, if equal to 1.0 signify correctness. [for casenos 1-11
only at present.]
A menu is provided which enables:
* 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.
In order to ensure accuracy, the grid distribution in beam-
thickness is such as to provide one very thin cell near each beam
surface.
ENDDIS
PHOTON USE
p;;;;
set prop off
cl
msg x-direction-displacement field
gr ou z 1
cont U1 z 1 fil;.0001
vec z 1 y 2 22 col 0
pause
cl
msg y-direction-displacement field
gr ou z 1
cont V1 z 1 fil;.0001
vec z 1 col 0
pause;cl
msg v1/v1_exact field
gr ou z 1
cont V1/t z 1 fi;.0001
pause;cl
msg u1/u1_exact field
gr ou z 1
cont u1/t z 1 fil;.0001
pause;cl
msg p1/p1_exact field
gr ou z 1
cont p1/t z 1 fil;.0001
ENDUSE
************************************************************
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
POISSON=0.0 ! Poisson's ratio
INTEGER(CASENO) ! start of menu
mesg(caseno 1 : default settings
mesg(caseno 3 : Poisson's ratio = 0.0
mesg(caseno 5 : zwlast = 1.0 m
mesg(caseno 7 : lx = 10.0 * ly
mesg(caseno 9 : not bending but tension
mesg(caseno 11 : point downward load at right-hand end
mesg(caseno 13 : point load but no rotation at RH end
mesg(caseno 15 : uniform downward load; RH end fixed
mesg(caseno 17 : no load; x-direction temperature gradient
mesg(caseno 19 : no load; y-direction temperature gradient
caseno=1
mesgm(caseno = :caseno: Enter another if not OK
readvdu(caseno, int, 1)
caseno
libref=caseno ! useful because this appears in graphical monitor
TEXT(bent beam in plane strain; s220
************************************************************
Group 2. Time dependence
STEADY = T
************************************************************
INTEGER(NXBODY,NYBODY) ! grid-fineness menu
NXBODY = 20
NYBODY = 20
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
mesga(Note that your choices will not be reflected in the
mesg(PHOTON USE file)
mesgm(insert other nxbody
readvdu(nxbody,int,:nxbody:)
nxbody
mesgm(insert other nybody, not less than 3
readvdu(nybody,int,:nybody:)
nybody
endif
mesga(grid intervals are uniform except for thin cells
mesgb(near beam surface. 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)
powery
endif
Group 3. X-Direction Grid Spacing
if(caseno.eq.7) 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) 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,DILT)
IF(CASENO.EQ.17.OR.CASENO.EQ.19) THEN
SOLVE(TEM1)
STORE(EPST,DVO1)
ENDIF
STORE(VISL,DVO1,DRH1)
************************************************************
GROUP 8. ITERATION NUMBERS ETC
RESFAC=1.e-8
RESREF(V1)=0.0
RESREF(U1)=0.0
RESREF(P1)=0.0
endit(v1)=0.0
liter(v1)= 100
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.POISSON.EQ.0.0) THEN
IPRPS=161
POISSON=0.0 ! 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(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) 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; s111
fx=0
ENDIF
IF(CASENO.EQ.19) 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; s111
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) then
form=:fx:*:LY: ! simple tension
TEXT(beam in simple tension; s220
endif
if(caseno.gt.10) then
(SOURCE of V1 at FORC01 is COVAL(FIXFLU,-0.1*:fx:)) ! downward force
form=0.0
TEXT(cantilevered beam with point load; s111
ENDIF
if(caseno.eq.15) 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) 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(rightend,cell,nx-1,nx-1,2,ny-1,1,1,1,1) ! better end condition for u
if(nybody.gt.1) then
(source of u1 at rightend is coval(1.e10,u1t))
endif
patch(rightenv,ewall,nx-1,nx-1,1,ny-1,1,1,1,1) ! better end condition for v
(source of v1 at rightenv is coval(1.0,v1t))
patch(rightenv,cell,nx-1,nx-1,1,ny-1,1,1,1,1) ! better end condition for v
(source of v1 at rightenv is coval(1.0e10,v1t))
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
PATCH(AXESVV,cell,2,2,1,ny-1,1,1,1,1) ! V1 fixed at bottom on left
PATCH(AXESVV,cell,2,3,1,ny-1,1,1,1,1) ! V1 fixed at bottom on left
COVAL(AXESVV,V1,FIXVAL,0.0)
(SOURCE OF V1 AT AXESVV IS COVAL(1.E10,V1T))
! PLANE-STRAIN, EPSZ = 0
SPEDAT(BOUNDARY,ZCONST,R,1.0e20)
if(caseno.eq.13) then
TEXT(cantilever; no end rotation: s12
endif
if(caseno.eq.15) then
TEXT(uniform load; no end rotation: s111
endif
title
************************************************************
GROUP 15. TERMINATE SWEEPS
LSWEEP = 1600
if(caseno.ne.7) then
spedat(rlxfac,rlxu1d,r,0.5)
spedat(rlxfac,rlxv1d,r,0.1)
else
spedat(rlxfac,rlxu1d,r,1.5)
spedat(rlxfac,rlxv1d,r,1.5)
endif
lsweep
ISG21=LSWEEP
************************************************************
GROUP 17. RELAXATION
#CONPROM
RELAX(P1 ,LINRLX, 1.000000E+00)
************************************************************
GROUP 19. DATA TRANSMITTED TO GROUND
STRA = T
PARSOL = F
ISG52 = 3 ! probe & res
************************************************************
GROUP 23.FIELD PRINT-OUT & PLOT CONTROL
output(drh1,n,n,n,n,n,n)
output(visl,n,n,n,n,n,n)
patch(vprof,profil,1,nx,1,1,1,1,1,1)
coval(vprof,v1,0.0,0.0)
TSTSWP = - 1 ! graphic-mode
NXPRIN=1; ! IXPRL=NX; IXPRF=NX-5
NYPRIN=1; ! IYPRL=NY; IYPRF=NY-5
IXMON = NX-2
IYMON = ny/2
IZMON = 1
#maxmin
#endpause
#$s003
inform7begin
real(RA,RB,RC,RD)
RA = FX/YOUNG*(1-POISSON**2) ! reference is required for the origin
RB = -FX/YOUNG/2*(1-POISSON**2) ! of these formulae
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:) ! YV was incorrect
FR4=(Yg-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) then
fr3=:RA:*(XU-0.01*:LX:)*:LY:
fr1=(:rb:*:ly:*(YV-0.01*:ly:)*:rc:
fr2=*2.0)
fr4=:ly:
endif
Exact solutions not yet provided for cases 11 and above
(STORED VAR V1T IS :FR1::FR2:)
coval(vprof,v1t,0.0,0.0)
(STORED VAR DILT IS :RD:*:FR4: with imat>100)
(STORED VAR U1T IS :FR3: with imat>100)
(STORED VAR u1/t IS u1/u1t with imat>100)
(STORED VAR v1/t IS v1/v1t)
(STORED VAR p1/t IS p1/dilt with imat>100)
inform7end
STOP