Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

All published worksheets from http://sagenb.org

Views: 168770
Image: ubuntu2004

Analiza statyczna płyty cienkiej z wykorzystaniem Metody Elementów Skończonych

# -*- coding: utf-8 -*-
"""
Created on Mon Dec 21 17:33:38 2009

@author: krystian.rosinski [at] gmail [dot] com
"""

import numpy as np

E = 31000000
Ni = 0.20

aa=12
bb=16  
h1=0.2 
h2=0.3  
la=12  
lb=16  
q1=4  
q2=6   

a=aa/la     
b=bb/lb     
le=la*lb     
nr_el=np.arange(0,le).reshape(lb,la)

l_wez=(la+1)*(lb+1)
nr_wez=np.arange(0,l_wez).reshape(lb+1,la+1)

k=0

wez_elem = np.zeros((la*lb,4))

for i in range (0,lb):
    for j in range (0,la):
        m_pom=nr_wez[np.ix_([i,i+1],[j,j+1])]
        wez_elem[k,0]=m_pom[1,0]
        wez_elem[k,1]=m_pom[1,1]
        wez_elem[k,2]=m_pom[0,1]
        wez_elem[k,3]=m_pom[0,0]
        k=k+1

#nr_st_sw=[]

for i in range (0,4):
    pom1=np.array(([wez_elem[:,i]])).transpose()
    if i==0:
        nr_st_sw=np.column_stack((3*(pom1)+0,3*(pom1)+1,3*(pom1)+2))
    else:
        nr_st_sw=np.column_stack((nr_st_sw,3*(pom1)+0,3*(pom1)+1,3*(pom1)+2))

Kij = np.zeros((3,3))

def K_ij(E,H,Ni,a,i,j):
    D=(E*H**3)/(12*(1-Ni)**2)
    if (i==0)|(i==3):
        ksi_i=-1
    else:
        ksi_i=1
    if (i==0)|(i==1):
        eta_i=-1
    else:
        eta_i=1
    if (j==0)|(j==3):
        ksi_j=-1
    else:
        ksi_j=1
    if (j==0)|(j==1):
        eta_j=-1
    else:
        eta_j=1
    
    Kij[0,0]=((D*ksi_i*ksi_j*eta_i*eta_j)/((a**2)*700))*(1014+525*(ksi_i*ksi_j+eta_i*eta_j))
    Kij[0,1]=-((D*ksi_i*eta_i*eta_j)/(a*700))*(384+175*(ksi_i*ksi_j+3*eta_i*eta_j)+210*Ni*(1+ksi_i*ksi_j))
    Kij[0,2]=-((D*ksi_i*ksi_j*eta_i)/(a*700))*(384+175*(3*ksi_i*ksi_j+eta_i*eta_j)+210*Ni*(1+eta_i*eta_j))
    Kij[1,1]=((D*eta_i*eta_j)/700)*(349+295*ksi_i*ksi_j+175*(3+ksi_i*ksi_j)*eta_i*eta_j)
    Kij[1,2]=((D*ksi_j*eta_i)/700)*(104+175*(ksi_i*ksi_j+eta_i*eta_j)+7*Ni*((6+5*ksi_i*ksi_j)*(6+5*eta_i*eta_j)-1))
    Kij[2,2]=((D*ksi_i*ksi_j)/700)*(349+295*eta_i*eta_j+175*(3+eta_i*eta_j)*ksi_i*ksi_j)
    Kij[1,0]=ksi_i*ksi_j*Kij[0,1]
    Kij[2,0]=eta_i*eta_j*Kij[0,2]
    Kij[2,1]=ksi_i*ksi_j*eta_i*eta_j*Kij[1,2]

K_e = np.zeros((2,12,12))

for i in range(0,4):
    for j in range(0,4):
        pom4=3*i
        pom5=3*j
        K_ij(E,h1,Ni,a,i,j)
        K_e[0,pom4:(pom4+3),pom5:(pom5+3)]=Kij
        K_ij(E,h2,Ni,a,i,j)
        K_e[1,pom4:(pom4+3),pom5:(pom5+3)]=Kij

Q1_1=np.array([[q1*a**2], [-(1./3.)*q1*a**3*(-1)], [-(1./3.)*q1*a**3*(-1)]])
Q1_2=np.array([[q1*a**2], [-(1./3.)*q1*a**3*1], [-(1./3.)*q1*a**3*(-1)]])
Q1_3=np.array([[q1*a**2], [-(1./3.)*q1*a**3*1], [-(1./3.)*q1*a**3*1]])
Q1_4=np.array([[q1*a**2], [-(1./3.)*q1*a**3*(-1)], [-(1./3.)*q1*a**3*1]])
Q_e1=np.concatenate((Q1_1,Q1_2,Q1_3,Q1_4))
Q2_1=np.array([[q2*a**2], [-(1./3.)*q2*a**3*(-1)], [-(1./3.)*q2*a**3*(-1)]])
Q2_2=np.array([[q2*a**2], [-(1./3.)*q2*a**3*1], [-(1./3.)*q2*a**3*(-1)]])
Q2_3=np.array([[q2*a**2], [-(1./3.)*q2*a**3*1], [-(1./3.)*q2*a**3*1]])
Q2_4=np.array([[q2*a**2], [-(1./3.)*q2*a**3*(-1)], [-(1./3.)*q2*a**3*1]])
Q_e2=np.concatenate((Q2_1,Q2_2,Q2_3,Q2_4))
Q_e=np.concatenate((Q_e1,Q_e2),1) #column_stack

pom2=np.shape(nr_st_sw)
pom3=np.zeros((pom2[0],2))

pom3[48:54,0]=1
pom3[60:66,0]=1
pom3[72:78,0]=1
pom3[84:90,0]=1
pom3[96:102,0]=1
pom3[108:114,0]=1
pom3[120:126,0]=1
pom3[132:138,0]=1

pom3[72:192,1]=1

nr_st_sw=np.column_stack((nr_st_sw,pom3))

K=np.zeros((3*l_wez,3*l_wez))

for element in range (0,le):
    for wiersz in range (0,12):
        for kolumna in range (0,12):
            wierszG=nr_st_sw[element,wiersz]
            kolumnaG=nr_st_sw[element,kolumna]
            grubosc=nr_st_sw[element,12]
            K[wierszG,kolumnaG]=K[wierszG,kolumnaG]+K_e[grubosc,wiersz,kolumna]

Q=np.zeros((3*l_wez,1))

for element in range (0,le):
    for wiersz in range (0,12):
        wierszG=nr_st_sw[element,wiersz]
        obciazenie=nr_st_sw[element,13]
        Q[wierszG]=Q[wierszG]+Q_e[wiersz,obciazenie]

for i in range (0,625,39):
    for j in range(0,3):
        K[i+j,:]=0
        K[:,i+j]=0
        K[i+j,i+j]=1
        Q[i+j]=0

for i in range (36,661,39):
    for j in range (0,3):
        K[i+j,:]=0
        K[:,i+j]=0
        K[i+j,i+j]=1
        Q[i+j]=0
    
for i in range (3,34,3):
    for j in range (0,2):
        K[i+j,:]=0
        K[:,i+j]=0
        K[i+j,i+j]=1
        Q[i+j]=0

for i in range (627,658,3):
    for j in range (0,2):
        K[i+j,:]=0
        K[:,i+j]=0
        K[i+j,i+j]=1
        Q[i+j]=0
        

u = np.linalg.solve(K,Q) # u=K\Q

w=u[::3,:]
v_x=u[1::3,:]
v_y=u[2::3,:]
m_w=w.reshape(lb+1,la+1).transpose()
m_v_x=v_x.reshape(lb+1,la+1).transpose()
m_v_y=v_y.reshape(lb+1,la+1).transpose()

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
from pylab import *

X=np.linspace(0,13,13)
Y=np.linspace(0,17,17)
X, Y = np.meshgrid(X, Y)

m_w=m_w.transpose()
m_v_x=m_v_x.transpose()
m_v_y=m_v_y.transpose()

Wariant I

fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X, Y, -m_w, rstride=1, cstride=1)
plt.show()
savefig('m_w1.png')

fig = plt.figure()
ax = Axes3D(fig)
#ax.plot_surface(X, Y, m_v_x, rstride=1, cstride=1, cmap=cm.jet)
ax.plot_wireframe(X, Y, m_v_x, rstride=1, cstride=1)
plt.show()
savefig('m_v_x1.png')

fig = plt.figure()
ax = Axes3D(fig)
#ax.plot_surface(X, Y, m_v_y, rstride=1, cstride=1, cmap=cm.jet)
ax.plot_wireframe(X, Y, m_v_y, rstride=1, cstride=1)
plt.show()
savefig('m_v_y1.png')

K=np.zeros((3*l_wez,3*l_wez))

for element in range (0,le):
    for wiersz in range (0,12):
        for kolumna in range (0,12):
            wierszG=nr_st_sw[element,wiersz]
            kolumnaG=nr_st_sw[element,kolumna]
            grubosc=nr_st_sw[element,12]
            K[wierszG,kolumnaG]=K[wierszG,kolumnaG]+K_e[grubosc,wiersz,kolumna]

Q=np.zeros((3*l_wez,1))

for element in range (0,le):
    for wiersz in range (0,12):
        wierszG=nr_st_sw[element,wiersz]
        obciazenie=nr_st_sw[element,13]
        Q[wierszG]=Q[wierszG]+Q_e[wiersz,obciazenie]

for i in range (156,469,39):
    for j in range(0,3):
        K[i+j,:]=0
        K[:,i+j]=0
        K[i+j,i+j]=1
        Q[i+j]=0
   
for i in range (0,37,3):
    for j in range (0,2):
        K[i+j,:]=0
        K[:,i+j]=0
        K[i+j,i+j]=1
        Q[i+j]=0

for i in range (642,661,3):
    for j in range (0,2):
        K[i+j,:]=0
        K[:,i+j]=0
        K[i+j,i+j]=1
        Q[i+j]=0
        

u = np.linalg.solve(K,Q) # u=K\Q

w=u[::3,:]
v_x=u[1::3,:]
v_y=u[2::3,:]
m_w=w.reshape(lb+1,la+1).transpose()
m_v_x=v_x.reshape(lb+1,la+1).transpose()
m_v_y=v_y.reshape(lb+1,la+1).transpose()

m_w=m_w.transpose()
m_v_x=m_v_x.transpose()
m_v_y=m_v_y.transpose()

Wariant II

fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X, Y, -m_w, rstride=1, cstride=1)
plt.show()
savefig('m_w2.png')

fig = plt.figure()
ax = Axes3D(fig)
#ax.plot_surface(X, Y, m_v_x, rstride=1, cstride=1, cmap=cm.jet)
ax.plot_wireframe(X, Y, m_v_x, rstride=1, cstride=1)
plt.show()
savefig('m_v_x2.png')

fig = plt.figure()
ax = Axes3D(fig)
#ax.plot_surface(X, Y, m_v_y, rstride=1, cstride=1, cmap=cm.jet)
ax.plot_wireframe(X, Y, m_v_y, rstride=1, cstride=1)
plt.show()
savefig('m_v_y2.png')