Notes

{ AXISYMMETRIC_STRESS.PDE }
{ ******************************************************

This example shows the application of FlexPDE to problems in
axi-symmetric stress.

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

dr(r*Sr)/r - St/r + dz(Trz) + Fr = 0
dr(r*Trz)/r + dz(Sz) + Fz = 0

where Sr, St and Sz are the stresses in the r- theta- and z- directions,
Trz is the shear stress, and Fr and Fz are the body forces in the
r- and z- directions.

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

er = dr(U)
et = U/r
ez = dz(V)
grz = dz(U) + dr(V).

The quantities U,V,er,et,ez,grz,Sr,St,Sz and Trz are related through the
constitutive relations of the material,

Sr = C11*er + C12*et + C13*ez - b*Temp
St = C12*er + C22*et + C23*ez - b*Temp
Sz = C13*er + C23*et + C33*ez - b*Temp
Trz = C44*grz

In isotropic solids we can write the constitutive relations as

C11 = C22 = C33 = G*(1-nu)/(1-2*nu) = C1
C12 = C13 = C23 = G*nu/(1-2*nu) = C2
b = alpha*G*(1+nu)/(1-2*nu)
C44 = G/2

where G = E/(1+nu) is the Modulus of Rigidity
E is Young's Modulus
nu is Poisson's Ratio
and alpha is the thermal expansion coefficient.

from which

Sr = C1*er + C2*(et + ez) - b*Temp
St = C1*et + C2*(er + ez) - b*Temp
Sz = C1*ez + C2*(er + et) - b*Temp
Trz = C44*grz

Combining all these relations, we get the displacement equations:

dr(r*Sr)/r - St/r + dz(Trz) + Fr = 0

dr(r*Trz)/r + dz(Sz) + Fz = 0

These can be written as

div(P) = St/r - Fr
div(Q) = -Fz

where P = [Sr,Trz]
and Q = [Trz,Sz]

The natural (or "load") boundary condition for the U-equation defines the
outward surface-normal component of P, while the natural boundary condition
for the V-equation defines the surface-normal component of Q. Thus, the
natural 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.


The problem analyzed here is a steel doughnut of rectangular cross-section,
supported on the inner surface and loaded downward on the outer surface.

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

title "The spherical shell"

coordinates
ycylinder('R','Z')

select

errlim = 1e-5

cubic


variables
U { declare U and V to be the system variables }
V

definitions
alpha0=3
beta0=80
alpha=alpha0*Pi/180
beta=beta0*Pi/180
rho=sqrt(r^2+z^2)
!theta=arctan(z/rho)
cost=r/rho
sint=z/rho
r1 = 1.0
r2 =1.04
q = 1e4
L = 3.0
r0=(r1+r2)/2

nu = 0.3
En = (20e10)
lambd=nu*En/((1+nu)*(1-2*nu))
mu=En/(2*(1+nu))
C1=lambd+2*mu
C2=lambd




Fr = 0
Fz = 0
mm=1e8




b11 = 5.13e11
b22 = 12e11
b33 = 12e11
b12 = 4.42e11
b13 =4.42e11
b23 = 4.81e11
b55 = 1.85e11



er = dr(U)
et = U/r
ez = dz(V)
grz = dz(U) + dr(V)

Ur=U*cost+V*sint
Ut=-U*sint+V*cost
n1=cost
n2=sint
m1=-sint
m2=cost

Sr = ((cost^2*b11+sint^2*b12)*cost^2+(cost^2*b12+sint^2*b22)*sint^2+4*cost^2*sint^2*b55)*er+(cost^2*b12+sint^2*b23)*et+((cost^2*b11+sint^2*b12)*sint^2+(cost^2*b12+sint^2*b22)*cost^2-4*cost^2*sint^2*b55)*ez+(-(cost^2*b11+sint^2*b12)*cost*sint+(cost^2*b12+sint^2*b22)*cost*sint+2*cost*sint*b55*(cost^2-sint^2))*grz
St = (cost^2*b12+sint^2*b23)*er+b22*et+(sint^2*b12+cost^2*b23)*ez+(-cost*sint*b12+cost*sint*b23)*grz
Sz = ((sint^2*b11+cost^2*b12)*cost^2+(sint^2*b12+cost^2*b22)*sint^2-4*cost^2*sint^2*b55)*er+(sint^2*b12+cost^2*b23)*et+((sint^2*b11+cost^2*b12)*sint^2+(sint^2*b12+cost^2*b22)*cost^2+4*cost^2*sint^2*b55)*ez+(-(sint^2*b11+cost^2*b12)*cost*sint+(sint^2*b12+cost^2*b22)*cost*sint-2*cost*sint*b55*(cost^2-sint^2))*grz
Trz = ((-cost*sint*b11+cost*sint*b12)*cost^2+(-cost*sint*b12+cost*sint*b22)*sint^2+2*cost*sint*b55*(cost^2-sint^2))*er+(-cost*sint*b12+cost*sint*b23)*et+((-cost*sint*b11+cost*sint*b12)*sint^2+(-cost*sint*b12+cost*sint*b22)*cost^2-2*cost*sint*b55*(cost^2-sint^2))*ez+(-(-cost*sint*b11+cost*sint*b12)*cost*sint+(-cost*sint*b12+cost*sint*b22)*cost*sint+(cost^2-sint^2)^2*b55)*grz

rt1=Sr*n1+Trz*n2
rt2=Trz*n1+Sz*n2
SigmaR=rt1*n1+rt2*n2
TauR=-rt1*n2+rt2*n1
tt1=Sr*m1+Trz*m2
tt2=Trz*m1+Sz*m2
SigmaT=-(tt1*m1+tt2*m2)
TauT=-(-tt1*m2+tt2*m1)


initial values
U = 0
V = 0

equations { define the axi-symmetric displacement equations }

U: dr(r*Sr)/r - St/r + dz(Trz) + Fr = 0
V: dr(r*Trz)/r + dz(Sz) + Fz = 0

constraints { prevent rigid-body motion: }

integral(V) = 0


boundaries
region 1
start "BDRY" (r1*cos(alpha),r1*sin(alpha))
load(U) = q* rho^2*sin(Pi/9)*sin(alpha)
load(V) = -q*rho^2*sin(Pi/9)*cos(alpha)


line to (r2*cos(alpha),r2*sin(alpha))
value(U) = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
value(V) = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

ARC ( CENTER = 0,0 ) ANGLE=beta0-alpha0

load(U) = q* rho^2*(cos(Pi/30)-1)*cos(beta) { define a free boundary on inner wall }
load(V) = q* rho^2*(cos(Pi/30)-1)*sin(beta)


line to (r1*cos(beta),r1*sin(beta))


value(U) = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
value(V) = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



ARC ( CENTER = 0,0 ) ANGLE=-beta0+alpha0


monitors
grid(r+mm*U,z+mm*V)

plots
grid(r+mm*U,z+mm*V)


contour(Ur) painted as "Ur-Displacement"
elevation(Ur,0) FROM (r1*cos((alpha+beta)/4),r1*sin((alpha+beta)/4)) TO (r2*cos((alpha+beta)/4),r2*sin((alpha+beta)/4))
{elevation(Ur,0) FROM (r1*cos((alpha+beta)/2),r1*sin((alpha+beta)/2)) TO (r2*cos((alpha+beta)/2),r2*sin((alpha+beta)/2))}
!!!!!!!!!!!!!!!!!!
elevation(Ur,0) FROM (r1*cos((alpha+beta)/2),r1*sin((alpha+beta)/2)) TO (r2*cos((alpha+beta)/2),r2*sin((alpha+beta)/2)) export format "#x#b#1" file="file_bir_Ur.txt"
!!!!!!!!!!!!!!!!!
elevation(Ur,0) FROM (r1*cos((alpha+beta)/1.5),r1*sin((alpha+beta)/1.5)) TO (r2*cos((alpha+beta)/1.5),r2*sin((alpha+beta)/1.5))

contour(Ut) painted as "Ut-Displacement" { show displacement field }
elevation(Ut,0) FROM (r1*cos((alpha+beta)/4),r1*sin((alpha+beta)/4)) TO (r2*cos((alpha+beta)/4),r2*sin((alpha+beta)/4))

!!!!!!!!!!!!!!!!!!
elevation(Ut,0) FROM (r1*cos((alpha+beta)/2),r1*sin((alpha+beta)/2)) TO (r2*cos((alpha+beta)/2),r2*sin((alpha+beta)/2)) export format "#x#b#1" file="file_bir_Ut.txt"
!!!!!!!!!!!!!!!!
elevation(Ut,0) FROM (r1*cos((alpha+beta)/1.5),r1*sin((alpha+beta)/1.5)) TO (r2*cos((alpha+beta)/1.5),r2*sin((alpha+beta)/1.5))

contour(SigmaR) painted as "R-Stress" { show displacement field }
elevation(SigmaR,0) FROM (r1*cos((alpha+beta)/4),r1*sin((alpha+beta)/4)) TO (r2*cos((alpha+beta)/4),r2*sin((alpha+beta)/4))
elevation(SigmaR,0) FROM (r1*cos((alpha+beta)/2),r1*sin((alpha+beta)/2)) TO (r2*cos((alpha+beta)/2),r2*sin((alpha+beta)/2))
!!!!!!!!!!!!!!!!!!
elevation(SigmaR,0) FROM (r1*cos((alpha+beta)/2),r1*sin((alpha+beta)/2)) TO (r2*cos((alpha+beta)/2),r2*sin((alpha+beta)/2)) export format "#x#b#1" file="file_bir_SigmaR.txt"
!!!!!!!!!!!!!!!!!
elevation(SigmaR,0) FROM (r1*cos((alpha+beta)/1.5),r1*sin((alpha+beta)/1.5)) TO (r2*cos((alpha+beta)/1.5),r2*sin((alpha+beta)/1.5))

contour(TauR) painted as "RT-Stress" { show displacement field }
elevation(TauR,0) FROM (r1*cos((alpha+beta)/2),r1*sin((alpha+beta)/2)) TO (r2*cos((alpha+beta)/2),r2*sin((alpha+beta)/2))
!!!!!!!!!!!!!!!!!!
elevation(TauR,0) FROM (r1*cos((alpha+beta)/2),r1*sin((alpha+beta)/2)) TO (r2*cos((alpha+beta)/2),r2*sin((alpha+beta)/2)) export format "#x#b#1" file="file_bir_TauR.txt"
!!!!!!!!!!!!!!!!!

contour(SigmaT) painted as "T-Stress" { show displacement field }
elevation(SigmaT) FROM (r1*cos((alpha+beta)/2),r1*sin((alpha+beta)/2)) TO (r2*cos((alpha+beta)/2),r2*sin((alpha+beta)/2))
!!!!!!!!!!!!!!!!!!
elevation(SigmaT,0) FROM (r1*cos((alpha+beta)/2),r1*sin((alpha+beta)/2)) TO (r2*cos((alpha+beta)/2),r2*sin((alpha+beta)/2)) export format "#x#b#1" file="file_bir_SigmaT.txt"
!!!!!!!!!!!!!!!!!



vector(U,V) painted as "Displacement" { show displacement field }


end 109550200