
Define PDE EulerBernoulli Beam



Parametrs of piezoelectric and cantilever beam


> 
Ys := 70*10^9: # Young's Modulus structure

> 
Yp := 11.1*10^10: # Young's Modulus pieazo

> 
ha := 0.00125: # Position

> 
hb := 0.001: # Position

> 
hc := 0.0015: # Position

> 
d31 := 180*10^(12): # Piezoelectric constant

> 
b := 0.01: #Width of the beam

> 
epsilon33 := 15.92*10^(9):

> 
hp :=0.00025: # Position

> 
hpc := 0.00125: # Position

> 
YI := b*(Ys*(hb^3 ha^3)+Yp*(hc^3hb^3))/3: # Bending stiffness of the composit cross section

> 
cs := 0.564: # The equivqlent coefficient of strain rate damping

> 
ca := 0: # Viscous air damping coefficient

> 
Ibeam := (b * tb^3 )/12: # The equivalent moment of inertia

> 
m := 0.101: # Mass of the structure

> 
upsilon :=  Yp*d31*b*(hc^2hb^2)/(2*hp): # Coupling term

> 
lb := 0.57:# Length of the structure (Cantilever Beam)

> 
lp := 0.05:# Length of the Piezoelectric

> 
R:= 10000: # Shunted resistor


Electrical circuit equation


> 
PDE1:=(epsilon33 * b*lp / hp) * diff(v(t), t) + (v(t)/R)+ int(d31*Yp*hpc*b* diff(w(x, t),$(x, 2))*diff(w(x, t), t),x = 0..lp)=0;


(1.1.1.1) 
> 




PDE Equation


> 
fn := 3.8:# Direct Excitation frequency;

> 
wb(x,t) := 0.01*sin(fn*2*Pi*t):#Direct Excitation;

> 
plot(wb(x,t),t = 0 .. 0.25*Pi,labels = [t,wb], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 15], linestyle = [longdash], axesfont = ["HELVETICA", "ROMAN", 10], legendstyle = [font = ["HELVETICA", 10], location = right],color = black);

> 
FunctionAdvisor(definition, Dirac(n,x));


(1.2.1) 
> 
PDE2 := YI*diff(w(x, t),$(x, 4))+ cs*Ibeam*diff(w(x, t),$(x, 4))*diff(w(x, t), t)+ ca* diff(w(x, t), t) + m * diff(w(x, t),$(t, 2))+ upsilon*v(t)*(Dirac(1,x) Dirac(1,xlp) ) =m*diff(wb(x, t),$(t, 2))ca*diff(wb(x, t), t);#PDE


(1.2.2) 
> 
N := 20:#NUMBER OF NODE POINT

> 
bc1 := dw(xmin, t) = 0:

> 
bc2 := dw(xmax, t) = 0:




Maple's pdsolve command


> 
bcs := { w(x,0)=0 , D[2](w)(x,0)=0 , w(0, t) = rhs(bc1), D[1](w)(0, t)= rhs(bc1), D[1,1](w)(lb,t) = rhs(bc2), D[1,1,1](w)(lb,t) = rhs(bc2)}; # Boundary conditions for PDE2.


(2.1) 
> 
PDES := pdsolve(PDE2, bcs, numeric, time = t, range = 0 .. xmax, indepvars = [x, t], spacestep = (1/1000)*xmax, timestep = (1/1000)*tmax);


(2.2) 
> 
PDES:plot3d(t = 0 .. tmax, x = 0 .. xmax, axes = boxed, orientation = [120, 40], shading = zhue, transparency = 0.3);


