>	restart:with(LinearAlgebra):with(plots):

 

>	#*****************************************

 

>	# Vibration of Mindlin rectangular plates

 

>	#----------------------------------------

 

>	# Created: Milan Batista (milan.batista@fpp.uni-lj.si)

 

>	# History: 21.3.2010

 

>

#

 

>	# This worksheat impliments the Ritz method for

 

>	# calculation of frequency factors for

 

>	# rectangular elastic thick plates. The details

 

>	# of method may be found in

 

>

#

 

>	# Liew, Wang et al - Vibration of Mindlin Plates

 

>	#                    Elsevier 1998, pp 93-95

 

>

#

 

>	# Note: in equation 4.26a kij,cd alpha is deleted

 

>

#

 

>	# The implementation conist of two procedures:

 

>	#  calc - calculate eigenvalues and eigenvectors

 

>	#  post - form solution for particular vibration mode

 

>

#

 

>	# User interface begins at DATA

 

>

#

 

>	#*****************************************

 

>	# Procedures

 

>	#*********************************************

 

>	calc:=proc(pmin,pmax,pstp,rtol,opt) local L1, L2, L3, R1, R2,R3, B1, B2, B3, T1, T2, T3, ff, i, j, k, K, Kcc, Kcd, Kce, Kdd, Kde, Kee, lc, Lp,M, Mcc, Mdd, Mee, m, N1,NN, q, t1, U, x, xx; global BC, p,alpha, tau, nu, kappa2,N,L,V,RE,ix,c,d,e,phi01, phi02, phi03,phi1,phi2, phi3;

 

>	L1:=0:L2:=0:L3:=0: # Free

 

>	R1:=0:R2:=0:R3:=0: # Free

 

>	B1:=0:B2:=0:B3:=0: # Free

 

>	T1:=0:T2:=0:T3:=0: # Free

 

>	if substring(BC,1)='S' then L1:=1;L2:=0;L3:=1;fi;

 

>	if substring(BC,1)='C' then L1:=1;L2:=1;L3:=1;fi; #------------ Left

 

>	if substring(BC,2)='S' then R1:=1;R2:=0;R3:=1;fi;

 

>	if substring(BC,2)='C' then R1:=1;R2:=1;R3:=1;fi; #------------ Right

 

>	if substring(BC,3)='S' then B1:=1;B2:=1;B3:=0;fi;

 

>	if substring(BC,3)='C' then B1:=1;B2:=1;B3:=1;fi; #------------- Bottom

 

>	if substring(BC,4)='S' then T1:=1;T2:=1;T3:=0;fi;

 

>	if substring(BC,4)='C' then T1:=1;T2:=1;T3:=1;fi; #------------- Top

 

>	#-------------------------------

 

>	# total number of coefficients

 

>	#------------------------------

 

>	lc:=0;for p from max(1,pmin) to pmax by max(1,pstp) do; lc:=lc+1;N:=(p+1)*(p+2)/2;

 

>	#-------------------------------

 

>	# Allocate space for coefficients

 

>	#----------------------------------

 

>	c:=vector(N,0):

 

>	d:=vector(N,0):

 

>	e:=vector(N,0):

 

>	#------------------------

 

>	# Allocate space for shape functions

 

>	#-----------------------------------

 

>	phi1:=vector(N,0):

 

>	phi2:=vector(N,0):

 

>	phi3:=vector(N,0):

 

>	#--------------------------------

 

>	# Define basic function

 

>	#---------------------------

 

>	phi01:=(xi+1)**L1*(eta+1)**B1*(xi-1)**R1*(eta-1)**T1;

 

>

 

>	phi02:=(xi+1)**L2*(eta+1)**B2*(xi-1)**R2*(eta-1)**T2;;

 

>

 

>	phi03:=(xi+1)**L3*(eta+1)**B3*(xi-1)**R3*(eta-1)**T3;;

 

>	#---------------------------------

 

>	# Define shape functions

 

>	#--------------------------

 

>	for q from 0 to p do;for i from 0 to q do; m:=(q+1)*(q+2)/2-i;phi1[m]:=xi**i*eta**(q-i)*phi01;od;od;

 

>	for q from 0 to p do;for i from 0 to q do; m:=(q+1)*(q+2)/2-i;phi2[m]:=xi**i*eta**(q-i)*phi02;od;od;

 

>	for q from 0 to p do;for i from 0 to q do; m:=(q+1)*(q+2)/2-i;phi3[m]:=xi**i*eta**(q-i)*phi03;od;od;

 

>	#-------------------------------

 

>	# Allocate space for matrices

 

>	#-------------------------------

 

>	Kcc:=matrix(N,N,0):

 

>	Kcd:=matrix(N,N,0):

 

>	Kce:=matrix(N,N,0):

 

>	#...........

 

>	Kdd:=matrix(N,N,0):

 

>	Kde:=matrix(N,N,0):

 

>	#.................

 

>	Kee:=matrix(N,N,0):

 

>	#.........................

 

>	Mcc:=matrix(N,N,0):

 

>	Mdd:=matrix(N,N,0):

 

>	Mee:=matrix(N,N,0):

 

>	#-------------------------------------------

 

>	# Define matrices

 

>	#-----------------------

 

>	t1:=6*(1-nu)*kappa2/4/tau**2;

 

>	for i from 1 to N do; for j from i to N do; Kcc[i,j]:=t1*int(int(alpha*diff(phi1[i],xi)*diff(phi1[j],xi)+1/alpha*diff(phi1[i],eta)*diff(phi1[j],eta),xi=-1..1),eta=-1..1);od;od;

 

>	for i from 1 to N do; for j from 1 to N do; Kcd[i,j]:=t1*int(int(diff(phi1[i],xi)*phi2[j],xi=-1..1),eta=-1..1);od;od;

 

>	for i from 1 to N do; for j from 1 to N do; Kce[i,j]:=t1/alpha*int(int(diff(phi1[i],eta)*phi3[j],xi=-1..1),eta=-1..1);od;od;

 

>	#............................

 

>	for i from 1 to N do; for j from i to N do; Kdd[i,j]:=int(int(alpha*diff(phi2[i],xi)*diff(phi2[j],xi)+(1-nu)/2/alpha*diff(phi2[i],eta)*diff(phi2[j],eta)+t1/alpha*phi2[i]*phi2[j],xi=-1..1),eta=-1..1);od;od;

 

>	for i from 1 to N do; for j from 1 to N do; Kde[i,j]:=int(int(nu*diff(phi2[i],xi)*diff(phi3[j],eta)+(1-nu)/2*diff(phi2[i],eta)*diff(phi3[j],xi),xi=-1..1),eta=-1..1);od;od;

 

>	#..................

 

>	for i from 1 to N do; for j from i to N do; Kee[i,j]:=int(int(1/alpha*diff(phi3[i],eta)*diff(phi3[j],eta)+(1-nu)/2*alpha*diff(phi3[i],xi)*diff(phi3[j],xi)+t1/alpha*phi3[i]*phi3[j],xi=-1..1),eta=-1..1);od;od;

 

>	#.........................

 

>	for i from 1 to N do; for j from i to N do; Mcc[i,j]:=1/16/alpha*int(int(phi1[i]*phi1[j],xi=-1..1),eta=-1..1);od;od;

 

>	for i from 1 to N do; for j from i to N do; Mdd[i,j]:=tau**2/48/alpha*int(int(phi2[i]*phi2[j],xi=-1..1),eta=-1..1);od;od;

 

>	for i from 1 to N do; for j from i to N do; Mee[i,j]:=tau**2/48/alpha*int(int(phi3[i]*phi3[j],xi=-1..1),eta=-1..1);od;od;

 

>	#---------------------------------

 

>	# Assemble matrices

 

>	#--------------------------------

 

>

NN:=3*N;

 

>	K:=Matrix(NN,NN,0,datatype=float):

 

>	M:=Matrix(NN,NN,0,datatype=float):

 

>	L:=Vector(NN,0):    # Eigenvalues

 

>	V:=matrix(NN,NN,0): # Eigenvectors

 

>	#--------------

 

>

# Data

 

>	#------------------------

 

>	#alpha:=1;tau:=0.1;nu:=0.3;kappa2:=5/6; #0.86667;#evalf(Pi**2/12); #<================================

 

>	#.....................

 

>	for i from 1 to N do; for j from i to N do; K[i,j]:=evalf(Kcc[i,j]);od;od;

 

>	for i from 1 to N do; for j from 1 to N do; K[i,j+N]:=evalf(Kcd[i,j]);od;od;

 

>	for i from 1 to N do; for j from 1 to N do; K[i,j+2*N]:=evalf(Kce[i,j]);od;od;

 

>	#.....................

 

>	for i from 1 to N do; for j from i to N do; K[i+N,j+N]:=evalf(Kdd[i,j]);od;od;

 

>	for i from 1 to N do; for j from 1 to N do; K[i+N,j+2*N]:=evalf(Kde[i,j]);od;od;

 

>	#...................

 

>	for i from 1 to N do; for j from i to N do; K[i+2*N,j+2*N]:=evalf(Kee[i,j]);od;od;

 

>	#...........................

 

>	for i from 1 to N do; for j from i to N do; M[i,j]:=evalf(Mcc[i,j]);od;od;

 

>	for i from 1 to N do; for j from i to N do; M[i+N,j+N]:=evalf(Mdd[i,j]);od;od;

 

>	for i from 1 to N do; for j from i to N do; M[i+2*N,j+2*N]:=evalf(Mee[i,j]);od;od;

 

>	#-------------------------------

 

>	# Fill lower part by symetry

 

>	#----------------------

 

>	for i from 1 to NN-1 do; for j from i+1 to NN do; K[j,i]:=K[i,j];M[j,i]:=M[i,j];od;od;

 

>	#print(K);

 

>	#print(M);

 

>	#-----------------------------

 

>

# Solve

 

>	#-------------------------

 

>	U:=Eigenvectors(K,M):

 

>	L:=Re(U[1]): # eigenvalues

 

>	V:=Re(U[2]): # eigenvectors

 

>	ff:=1;if opt=1 then ff:=Pi**2 fi;for i from 1 to NN do; L[i]:=evalf(Re(sqrt(L[i]))/ff);od:

 

>	#--------------------------------------

 

>	# Sort eigenvalues

 

>	#----------------------------

 

>	ix:=vector(NN,0): # indices for eigenvectors

 

>	for i from 1 to NN do;ix[i]:=i;od:

 

>

for i from 1 to NN-1 do;k:=i;x:=L[i];xx:=i;for j from i+1 to NN do;if L[j]<x then k:=j;x:=L[j];xx:=j; fi;od;L[k]:=L[i];L[i]:=x;ix[k]:=i;ix[i]:=xx;od:unassign('x'):print(p,L[1],L[2],L[3],L[4],L[5],L[6]); if lc > 1 then for i from 1 to N1 do; RE[i]:=abs(L[i]-Lp[i])/abs(L[i]); od; if max(RE[1..6])<rtol then print("rerr",RE[1],RE[2],RE[3],RE[4],RE[5],RE[6]); return fi;fi;if p < pmax then unassign('Lp','RE');N1:=N;Lp:=vector(N1,L[1..N1]);RE:=Vector(N1,0) fi;

 

>	#---------------

 

>	# Free space

 

>	#------------------

 

>	unassign('Kcc','Kcd','Kce','Kdd','Kde','Kee','Mcc','Mdd','Mee','K','M','U'); od;p:=p-pstp;

 

>	print("*** Error: rtol test faild:",RE[1],RE[2],RE[3],RE[4],RE[5],RE[6]);

 

>

end proc:

 

>	#==== END PROC CALC ======

 

>	post:=proc(J) local q,i, m, IX;global ix,V,N,L,c,d,e,w,px,py,phi1,phi2,phi3,Mx,My,Mxy,Qx,Qy, alpha, nu;

 

>	IX:=ix[J];print(L[J]):# Mode index and eigenvalue

 

>	#--------- Parameters

 

>	c:=V[1..N,IX]: # coefficients for w

 

>	d:=V[N+1..2*N,IX]: # coefficients for phi1

 

>	e:=V[2*N+1..3*N,IX]: # coefficients for phi2

 

>	#---------- Functions

 

>	w:=0:for q from 0 to p do: for i from 0 to q do: m:=(q+1)*(q+2)/2-i:w:=w+c[m]*phi1[m]:od;od:

 

>	px:=0:for q from 0 to p do: for i from 0 to q do: m:=(q+1)*(q+2)/2-i:px:=px+d[m]*phi2[m]:od;od:

 

>	py:=0:for q from 0 to p do: for i from 0 to q do: m:=(q+1)*(q+2)/2-i:py:=py+e[m]*phi3[m]:od;od:

 

>	#------------------------------------------------------

 

>	w:=subs(xi=x*alpha,eta=y,w): # physical dimensions

 

>	px:=subs(xi=x*alpha,eta=y,px): # physical dimensions

 

>	py:=subs(xi=x*alpha,eta=y,py): # physical dimensions

 

>	#================= Stress resultants

 

>	Mx:=0:Mx:=diff(px,x)+nu*diff(py,y):

 

>	My:=0:My:=diff(py,y)+nu*diff(px,x):

 

>	Mxy:=0:Mxy:=diff(px,y)+diff(py,x):

 

>	Qx:=0:Qx:=px+diff(w,x):

 

>	Qy:=0:Qy:=py+diff(w,y):

 

>	#------------------------------------

 

>

#

 

>

#

 

>

end proc:

 

>	#==== END PROC POST ======

 

>	#*************************************

 

>

#  DATA

 

>	#**********************

 

>	alpha:=0.5;tau:=0.15;nu:=0.3;kappa2:=0.86667;  # b/a ratio, h/b ratio, Poisson ratio, Mindlin factor squared <======== INPUT

 

(1)

 

(1)

 

(1)

 

(1)

 

>	BC:=SCFS; # Boundary conditions: Left Right Bottom Top / S-simple F-free C-Clamped <================================== INPUT

 

(2)

 

>	p1:=2;p2:= 10; p3:=1;   # degree of polynomial space: from p1 to p2 by p3          <================================== INPUT

 

(3)

 

(3)

 

(3)

 

>	rtol:=1e-3; # relative tolerance to stop                                           <================================== INPUT

 

(4)

 

>	opt:=0;  # frequency factor devided by Pi^2: 0=No, 1=yes                           <================================== INPUT

 

(5)

 

>	#******************

 

>

# Solve

 

>	#*****************

 

>	calc(p1,p2,p3,rtol,opt); # out: p, freqency factor 1..6

 

(6)

 

(6)

 

(6)

 

(6)

 

(6)

 

(6)

 

(6)

 

(6)

 

(6)

 

(6)

 

>	#*****************************

 

>	# Graphics

 

>	#*****************************

 

>	J:=6; # Input Mode number <========================================= INPUT

 

(7)

 

>	#----------------------------

 

>	post(J): # out:  frequency factor

 

(8)

 

>	# Contour plot of deflection

 

>	contourplot(evalc(w),x=-1/alpha..1/alpha,y=-1..1,grid=[50,50],filled=true,axes='none',scaling=constrained,coloring=[yellow,red]);

 

 

>	# Animation of plate vibration

 

>	animate( plot3d, [cos(t)*evalc(w),x=-1/alpha...1/alpha,y=-1..1], t=0..2*Pi, style=patch );

 

 

 

>	# Distribution of deflection

 

>	plot3d(w,x=-1/alpha...1/alpha,y=-1..1, grid=[50,50],style=patch,title="w",axes=boxed );

 

 

>	# Distribution of bending moment Mx

 

>	plot3d(Mx,x=-1/alpha...1/alpha,y=-1..1, grid=[50,50],style=patch,title="Mx",axes=boxed );

 

 

>	# Distribution of bending moment My

 

>	plot3d(My,x=-1/alpha...1/alpha,y=-1..1, grid=[50,50],style=patch,title="My",axes=boxed );

 

 

>	# Distribution of twisting moment Mxy

 

>	plot3d(Mxy,x=-1/alpha...1/alpha,y=-1..1, grid=[50,50],style=patch,title="Mxy",axes=boxed );

 

 

>	# Distribution of shear force Qx

 

>	plot3d(Qx,x=-1/alpha...1/alpha,y=-1..1, grid=[50,50],style=patch,title="Qx",axes=boxed );

 

 

>	# Distribution of shear force Qy

 

>	plot3d(Qy,x=-1/alpha...1/alpha,y=-1..1, grid=[50,50],style=patch,title="Qy",axes=boxed );

 

 

>	#### END #####

 

>