Chapter7-Numerical Methods

Ex3-pg204

In [12]:
import math
#calculate tauA and tauB
import numpy
from numpy.linalg import inv
a=15. ##mm
b=10. ##mm
h=5. ##mm
h1=4.4 ##mm
h2=2.45 ##mm
h3=3. ##mm

x=([[2, 0, 0, 0, 2, -4],[0, 2, 0, 1 ,-4 ,1],[0, 0, 2, -4, 1, 0],[-4, 2 ,0 ,0 ,0 ,1],[1, -4.27, 1, 0, 1.06, 0],[0, 1.25, -7.41, 1.34, 0, 0]])
print(x)
y=([[-2], [-2], [-2], [-2], [-2], [-2]])
print(y)
z=numpy.dot(inv(x),y)
X1=z[1]
X2=z[2]
X3=z[3]
X4=z[0]
X5=z[4]
X6=z[5]
print(X4)

print(X1)
print(X2)
print(X3)
print(X5)
print(X6)
 
dfi=2.075 
d3fi=-0.001
d2fi=-1.383
d4fi=0.002

##tauB=derivative(fi,y)B
tauB=(dfi+(d2fi/2.)-(d3fi/3.)+(d4fi/4.))
print'%s %.2f %s'%('tauB=',tauB,' G*thetab')

dfi=1.536
d2fi=-0.613
d3fi=-0.002
d4fi=0.001
d5fi=0.001
d6fi=-0.002

##tauA=derivative(fi,x)A
tauA=(dfi+(-d2fi/2.)-(d3fi/3.)-(d4fi/4.)+(d5fi/5.)+(d6fi/6.))
print'%s %.2f %s'%('tauA=',tauA,' G*thetaa')

[[2, 0, 0, 0, 2, -4], [0, 2, 0, 1, -4, 1], [0, 0, 2, -4, 1, 0], [-4, 2, 0, 0, 0, 1], [1, -4.27, 1, 0, 1.06, 0], [0, 1.25, -7.41, 1.34, 0, 0]]
[[-2], [-2], [-2], [-2], [-2], [-2]]
[ 2.07117939]
[ 1.76075859]
[ 0.84472513]
[ 1.53616791]
[ 2.45522136]
[ 2.76320038]
tauB= 1.38  G*thetab
tauA= 1.84  G*thetaa

Ex6-pg210

In [13]:
import math
#calculate effective fem at joint B and the change of moment in beam segment AB  and BC
p=15.
P=45.
a=3.
b=1.5
L1=3.
L2=4.5
MfAB=-(p*L1**2)/12.
print'%s %.2f %s'%("in kNm is= ",MfAB,"")
MfBA=(11.25)
print'%s %.2f %s'%("in kNm is= ",MfBA,"")
MfBC=-(P*a*b**2)/L2**2
print'%s %.2f %s'%("in kNm is= ",MfBC,"")
MfCB=(P*b*a**2)/L2**2
print'%s %.2f %s'%("in kNm is= ",MfCB,"")
B=MfBA+MfBC
print'%s %.2f %s'%("effective fem at joint B in kNm is= ",B,"")
AB=0.429*-B    ##                      joint rotates until a change in moment is +3.75
print'%s %.2f %s'%("the change of moment in beam segment AB in kN is=",AB,"")
BC=0.571*-B
print'%s %.2f %s'%("the change of moment in beam segment AB in kN is=",BC,"")

in kNm is=  -11.25 
in kNm is=  11.25 
in kNm is=  -15.00 
in kNm is=  30.00 
effective fem at joint B in kNm is=  -3.75 
the change of moment in beam segment AB in kN is= 1.61 
the change of moment in beam segment AB in kN is= 2.14

Ex7-pg221

In [2]:
#find epsilon
import math
import numpy
from numpy.linalg import inv
p=14. ##MPa
t=0.3 ##cm
E=200. ##GPa
v=0.3
gamma1=77. ##kN/m**3
alpha=12.*10**-6 ## per degree celcius
A=2
T=50. ##degree celcius

D=numpy.matrix([[3.33, 0.99, 0],[0.99 ,3.3, 0],[0, 0, 1.16]])
print(D)
##[D*]=(t*[D])/4*A
Dj=(10**6)*D

D1=Dj/8. 

print(D1)

##solution a: stiffness matrix
xi=0.
x1=0.
xj=4.
x2=4.
xm=0.
x3=0.
yi=-1.
y1=-1.
yj=-1.
y2=-1.
ym=1.
y3=1.

ai=0.-4.
a1=0.-4.
print(ai,a1)
aj=0-0.
a2=0-0.
print(aj,a2)
am=4.-0.
a3=4.-0.
print(am,a3)

bi=-1-1
b1=-1-1
print(bi,b1)
bj=1.+1.
b2=1.+1.
print(bj,b2)
bm=-1.+1.
b3=-1+1.
print(bm,b3)

k11=(10**6/8.)*(3.3*4+1.16*16)
print'%s %.2f %s'%('k11=',k11,'')
k12=(10**6/8)*(3.3*2*-2+0)
print'%s %.2f %s'%('k12=',k12,'')
k13=(10**6/8)*(0+1.16*4*-4)
print'%s %.2f %s'%('k13=',k13,'')
k22=(10**6/8.)*(3.3*4+0)
print'%s %.2f %s'%('k22=',k22,'')
k23=0.
print'%s %.2f %s'%('k23=',k23,'')
k32=0.
print'%s %.2f %s'%('k32=',k32,'')
k21=(10**6/8)*(3.3*2*-2+0)
print'%s %.2f %s'%('k21=',k21,'')
k31=(10**6/8)*(0+1.16*4*-4)
print'%s %.2f %s'%('k31=',k31,'')
k33=(10**6/8.)*(0+1.16*16)
print'%s %.2f %s'%('k33=',k33,'')

kuu=numpy.matrix([[k11 ,k12, k13],[k21 ,k22 ,k23],[k31, k32, k33]])
print(kuu)
kuv=10**6*numpy.matrix([[2.15, -1.16, -0.99],[-0.99, 0, 0.99],[-1.16 ,1.16, 0]])
print(kuv)
kvv=10**6*numpy.matrix([[7.18 ,-0.58, -6.6],[-0.58, 0.58 ,0],[-6.6, 0, 6.6]])
print(kvv)
kvu=numpy.matrix([[2.15, -0.99, -1.16],[-1.16, 0, 1.16],[-0.99, 0.99, 0]])
print(kvu)



##solution b:
Fx=0
Fy=0.077 ##N/cm**2
Qbe=[0,0,0,-0.0308,-0.0308,-0.0308]##N
print(Qbe)


Qp3=numpy.matrix([[0,-420,-420,0,-840,-840]])
print('Qp3')

epsilon=alpha*T
print'%s %.2f %s'%('epsilon=',epsilon,'')
##Qte=[B']*[D]*epsilon*At
Qte=(1/8.)*numpy.matrix([[-2, 0, -4],[2, 0, 0],[0, 0, 4],[0, -4, -2],[0 ,0, 2],[0, 4 ,0]])*((200*10**5)/0.91)*numpy.matrix([[1, 0.3, 0,],[0.3 ,1 ,0],[0, 0 ,0.35]])*numpy.matrix([[0.0006],[0.0006],[0]])*(1.2)
print(Qte)

Qe=numpy.matrix([[-5142.85,4742.85,-400,-10285.71,-840.03,9445.67]])
print (Qe)

[[ 3.33  0.99  0.  ]
 [ 0.99  3.3   0.  ]
 [ 0.    0.    1.16]]
[[ 416250.  123750.       0.]
 [ 123750.  412500.       0.]
 [      0.       0.  145000.]]
(-4.0, -4.0)
(0.0, 0.0)
(4.0, 4.0)
(-2, -2)
(2.0, 2.0)
(0.0, 0.0)
k11= 3970000.00 
k12= -1650000.00 
k13= -2320000.00 
k22= 1650000.00 
k23= 0.00 
k32= 0.00 
k21= -1650000.00 
k31= -2320000.00 
k33= 2320000.00 
[[ 3970000. -1650000. -2320000.]
 [-1650000.  1650000.        0.]
 [-2320000.        0.  2320000.]]
[[ 2150000. -1160000.  -990000.]
 [ -990000.        0.   990000.]
 [-1160000.  1160000.        0.]]
[[ 7180000.  -580000. -6600000.]
 [ -580000.   580000.        0.]
 [-6600000.        0.  6600000.]]
[[ 2.15 -0.99 -1.16]
 [-1.16  0.    1.16]
 [-0.99  0.99  0.  ]]
[0, 0, 0, -0.0308, -0.0308, -0.0308]
Qp3
epsilon= 0.00 
[[ -5142.85714286]
 [  5142.85714286]
 [     0.        ]
 [-10285.71428571]
 [     0.        ]
 [ 10285.71428571]]
[[ -5142.85   4742.85   -400.   -10285.71   -840.03   9445.67]]

Ex8-pg223

In [3]:
import math
#calculate epsilon and sigma 
import numpy
from numpy.linalg import inv
t=0.3 ##cm
E=200 ##GPa
v=0.3
i=2.
j=4.
m=3.
L=5000. ##N

a1=0-4
a2=0-4
print(a1,a2)
aj=4-0
a4=4-0
print(aj,a4)
am=4-4
a3=4-4
print(am,a3)

bi=1-1
b2=1-1
print(bi,b2)
bj=1+1
b4=1+1
print(bj,b4)
bm=-1-1
b3=-1-1
print(bm,b3)

k22=(10**6/8.)*(3.3*0+1.16*16)
print'%s %.2f %s'%('k22=\n',k22,'')
k44=(10**6/8.)*(3.3*4*+1.16*16)
print'%s %.2f %s'%('k44=\n',k44,'')
k24=(10**6/8.)*(3.3*0+1.16*4*-4)
print'%s %.2f %s'%('k24=\n',k24,'')
k42=(10**6/8)*(3.3*0+1.16*4*-4)
print'%s %.2f %s'%('k42=\n',k42,'')
k23=0
print'%s %.2f %s'%('k23=\n',k23,'')
k32=0
print'%s %.2f %s'%('k32=\n',k32,'')
k43=(10**6/8)*(3.3*2*-2+1.16*0)
print'%s %.2f %s'%('k43=\n',k43,'')
k34=(10**6/8)*(3.3*2*-2+1.16*0)
print'%s %.2f %s'%('k34=\n',k34,'')
k33=(10**6/8)*(3.3*4+1.16*0)
print'%s %.2f %s'%('k33=\n',k33,'')


kuu=numpy.matrix([[k22, k23, k24],[k32, k33, k34,],[k42, k43, k44]])
print(kuu)
kuv=10**6*numpy.matrix([[0 ,1.16 ,-1.16],[0.99 ,0 ,-0.99],[-0.99, -1.16, 2.15]])
print(kuv)
kvv=10**6*numpy.matrix([[6.6, 0 ,-6.6],[0, 0.58, -0.58],[-6.6 ,-0.58, 7.18]])
print(kvv)
kvu=10**6*numpy.matrix([[0 ,0.99, -0.99],[1.16, 0, -1.16],[-1.16 ,-0.99, 2.15]])
print(kvu)





k1=numpy.matrix([[3.97, -1.65, -2.32, 0],[-1.65, 1.65 ,0 ,0],[-2.32, 0, 2.32, 0],[0, 0 ,0 ,0]])
print(k1)
k2=numpy.matrix([[2.15, -1.16, -0.99, 0],[-0.99, 0, 0.99, 0],[-1.16, 1.16, 0, 0],[0, 0 ,0 ,0]])
print(k2)
k3=numpy.matrix([[2.15, -0.99 ,-1.16, 0],[-1.16, 0 ,1.16, 0],[-0.99 ,0.99, 0 ,0],[0 ,0, 0, 0]])
print(k3)
k4=numpy.matrix([[7.18 ,-0.58, -6.6, 0],[-0.58, 0.58, 0 ,0],[-6.6, 0, 6.6, 0],[0, 0 ,0, 0]])
print(k4)




k5=numpy.matrix([[0, 0 ,0 ,0],[0 ,2.32, 0, -2.32],[0, 0, 1.65, -1.65],[0 ,-2.32, -1.65, 3.97]])
print(k5)
k6=numpy.matrix([[0, 0, 0 ,0],[0, 0 ,1.16, -1.16],[0, 0.99, 0 ,-0.99],[0, -0.99, -1.16, 2.15]])
print(k6)
k7=numpy.matrix([[0, 0 ,0 ,0],[0, 0, 0.99, -0.99],[0 ,1.16, 0, -1.16],[0, -1.16, -0.99, 2.15]])
print(k7)
k8=numpy.matrix([[0 ,0, 0 ,0],[0 ,6.6 ,0 ,-6.6],[0, 0 ,0.58 ,-0.58],[0 ,-6.6 ,-0.58, 7.18]])
print(k8)





Qy4=((3.*(-5000.))/4.*1.)*((1./2.)*(1.+1.)+0.33*(-0.25*(1.-1.+1.)-0.75))
print'%s %.2f %s'%('Qy4=',Qy4,'')                                 ## textbook ans is wrong
Qy2=((3*(-5000))/4*1)*((1/2)*(1+1)-0.33*(1+0.75*(1-1+1)-0.75))
print'%s %.2f %s'%('Qy2=',Qy2,'')                                 ## textbook ans is wrong

Q=numpy.matrix([[0, 0, 0, 0 ,0 ,Qy4 ,0, Qy2]])
print(Q)
u1=0
u3=0
v1=0
v3=0

Z=numpy.matrix([[3.97, -2.32, 0, -1.16],[-2.32, 3.97, -0.99, 2.15],[0, -0.99, 7.18, -6.6],[-1.16, 2.15, -6.6, 7.18]])
print(Z)
z=numpy.linalg.inv(Z)
print(z)
X=z*numpy.matrix([[0],[0],[-2512.5],[-2512.5]]) 
print(X)
X1= X*10**-6
print("u2 u4 v2 v4 is= ",X1,"")

Y=numpy.matrix([[-2, 2, 0, 0, 0, 0],[0, 0, 0, -4, 0, 4],[-4 ,0, 4, -2, 2, 0]])
print(Y)
W=Y*numpy.matrix([[0],[-0.0012],[0],[0],[-0.0068],[0]])
print(W)
W1=W*(1./8.)
print("W1")

y=numpy.matrix([[1, 0.3, 0],[0.3, 1, 0],[0, 0, 0.35]])*W1
print(y)
u=(200.*10**9/0.91)
print(u)
U=u*y
print(U,"sigmax sigmay tauxy in Pa is= ")

(-4, -4)
(4, 4)
(0, 0)
(0, 0)
(2, 2)
(-2, -2)
k22=
 2320000.00 
k44=
 30624000.00 
k24=
 -2320000.00 
k42=
 -2320000.00 
k23=
 0.00 
k32=
 0.00 
k43=
 -1650000.00 
k34=
 -1650000.00 
k33=
 1650000.00 
[[  2320000.         0.  -2320000.]
 [        0.   1650000.  -1650000.]
 [ -2320000.  -1650000.  30624000.]]
[[       0.  1160000. -1160000.]
 [  990000.        0.  -990000.]
 [ -990000. -1160000.  2150000.]]
[[ 6600000.        0. -6600000.]
 [       0.   580000.  -580000.]
 [-6600000.  -580000.  7180000.]]
[[       0.   990000.  -990000.]
 [ 1160000.        0. -1160000.]
 [-1160000.  -990000.  2150000.]]
[[ 3.97 -1.65 -2.32  0.  ]
 [-1.65  1.65  0.    0.  ]
 [-2.32  0.    2.32  0.  ]
 [ 0.    0.    0.    0.  ]]
[[ 2.15 -1.16 -0.99  0.  ]
 [-0.99  0.    0.99  0.  ]
 [-1.16  1.16  0.    0.  ]
 [ 0.    0.    0.    0.  ]]
[[ 2.15 -0.99 -1.16  0.  ]
 [-1.16  0.    1.16  0.  ]
 [-0.99  0.99  0.    0.  ]
 [ 0.    0.    0.    0.  ]]
[[ 7.18 -0.58 -6.6   0.  ]
 [-0.58  0.58  0.    0.  ]
 [-6.6   0.    6.6   0.  ]
 [ 0.    0.    0.    0.  ]]
[[ 0.    0.    0.    0.  ]
 [ 0.    2.32  0.   -2.32]
 [ 0.    0.    1.65 -1.65]
 [ 0.   -2.32 -1.65  3.97]]
[[ 0.    0.    0.    0.  ]
 [ 0.    0.    1.16 -1.16]
 [ 0.    0.99  0.   -0.99]
 [ 0.   -0.99 -1.16  2.15]]
[[ 0.    0.    0.    0.  ]
 [ 0.    0.    0.99 -0.99]
 [ 0.    1.16  0.   -1.16]
 [ 0.   -1.16 -0.99  2.15]]
[[ 0.    0.    0.    0.  ]
 [ 0.    6.6   0.   -6.6 ]
 [ 0.    0.    0.58 -0.58]
 [ 0.   -6.6  -0.58  7.18]]
Qy4= -2512.50 
Qy2= 1237.50 
[[    0.      0.      0.      0.      0.  -2512.5     0.   1237.5]]
[[ 3.97 -2.32  0.   -1.16]
 [-2.32  3.97 -0.99  2.15]
 [ 0.   -0.99  7.18 -6.6 ]
 [-1.16  2.15 -6.6   7.18]]
[[ 0.42911692  0.17986034  0.25168713  0.24682604]
 [ 0.17986034  0.4831756  -0.25583582 -0.3507947 ]
 [ 0.25168713 -0.25583582  1.36613468  1.37304916]
 [ 0.24682604 -0.3507947   1.37304916  1.54633026]]
[[-1252.51435698]
 [ 1524.15919903]
 [-6882.19940602]
 [-7334.94080944]]
('u2 u4 v2 v4 is= ', matrix([[-0.00125251],
        [ 0.00152416],
        [-0.0068822 ],
        [-0.00733494]]), '')
[[-2  2  0  0  0  0]
 [ 0  0  0 -4  0  4]
 [-4  0  4 -2  2  0]]
[[-0.0024]
 [ 0.    ]
 [-0.0136]]
W1
[[ -3.00000000e-04]
 [ -9.00000000e-05]
 [ -5.95000000e-04]]
2.1978021978e+11
(matrix([[ -6.59340659e+07],
        [ -1.97802198e+07],
        [ -1.30769231e+08]]), 'sigmax sigmay tauxy in Pa is= ')

Ex9-pg226

In [3]:
import math
#calculate sigma 
L=76.2 ##mm
h=50.8 ##mm
t=25.4 ##mm
p=6895. ##kPa
E=207. ##GPa
v=0.15

##solution a: exact solution
##p=Mh/I
##sigmax=-(y/h)*p
sigmay=0.
tauxy=0.
##derivative(u,x)=-(yp/Eh)
##derivative(v,y)=(v*y*p)/(Eh)
##derivative(u,y)+derivative(v,x)=0
##u=-(p/E*h)*x*y                        ## for u(0,0)=v(0,0)=0 and u(L,0)=0
##v=-(p/2*E*h)*(x**2+v*y**2)
##sigmax=-(1/0.0508)*(y*p)
sigmaxmax=6895. ##kPa
##u(0.0762,-0.0254)=25.4*10**-6 ##m
##v(0.0762,0)=1.905*10**-6 ##m

##solution b:
Qx10=((0.0254*0.0254)/6.)*((2.*sigmaxmax)+3447.5)
print'%s %.2f %s'%("in mN is= ",Qx10,"")
Qx11=((0.0254*0.0254)/6.)*(2*3447.5+sigmaxmax)+((0.0254*0.0254)/6)*(2*3447.5+0.)
print'%s %.2f %s'%("in mN is= ",Qx11,"")
Qx12=((0.0254*0.0254)/6.)*(0+3447.5)
print'%s %.2f %s'%("in mN is= ",Qx12,"")

in mN is=  1.85 
in mN is=  2.22 
in mN is=  0.37