Discussion Forum

Hello everyone,

Iam designing piezoelectric energy harvester by analytical modelling technique, for determining voltage in matlab coding. While trying to execute it, iam getting voltage output, but the oscillations are not proper, and the time is increasing randomly upto 8000 seconds. I need to set the voltage [y axis] within -8 millivolts to +8 millivolts and time[x axis] within 0 to 10 milliseconds. I couldn't set the voltage and time in coding. I have attached my coding here. Please, tell me where i need to change in the cosding to get correct output.

clear all; close all;

%Beam properties

rho=5150; %beam density Ebeam=70e9; %young's modulus of beam Lbeam=0.003; %beam length width=0.001; %beam and piezo width Lf=0.200; %length to forcing f(t) thickbeam=0.0002; %beam thick Area=thickbeamwidth; %area I=widththickbeam^3/12; %Inertia

%Piezo peoperties g31=-9.5e-3; %piezoelectric coefficient Epiezo=50e9; %young's modulus of piezo thickpiezo=0.00005e-3; %piezo thick Lpiezo=0.003; %piezo length d31=-320e-12; %m/volt dielectric constant Rsource = 330000; %ohm source resistance c=(EbeamI/(rhoArea))^.5;

%Natural Frequencies

B = 0:0.01:50; cheq = cos(LbeamB).cosh(LbeamB); plot(B,cheq) axis([0 30 -8 8]) ezplot('-1',[0 30 -5 5]) axis([0 30 -5 5]) xlabel('Beta values'); title('Graph of characteristic equation') beta = [4.2;7.2;12.2;18.2;23.2]; wn=beta.^2c; fn=wn/(2*pi);

%ModeShapes

syms x for i=1:5, ModeShape(i) = cosh(beta(i)x)-cos(beta(i)x)-(sinh(beta(i)Lbeam)... -sin(beta(i)Lbeam))/(cosh(beta(i)Lbeam)+cos(beta(i)Lbeam))... (sinh(beta(i)x)-sin(beta(i)*x)); end

%Harmonic Forcing Function

w=wn; zeta=0.03; wd=wn(1-zeta^2)^.5; syms t tt Force=(0.35)sin(w(1)tt)1/(rhoArea)subs(ModeShape,x,Lf); for i=1:5 q(i)=1/wd(i)exp(-zetawn(i)t)... *int(Force(i)exp(zetawn(i)tt)sin(wd(i)(t-tt)),tt,0,t); end

%Develop voltage equations

for i=1:5 before_w3(i)=q(i)*ModeShape(i); end

w3=sum(before_w3); K_x_t=vpa(diff(w3,x,2),5); Kavg=vpa((1/Lpiezo)int(K_x_t,x,0,Lpiezo),5); tau = 2:.01:80; Moment=EbeamI.Kavg; psi=(Ebeamthickbeam)/(Epiezothickpiezo); T=thickbeam/thickpiezo; Veb_sym=(6g31psi(1+T))/(widththickpiezo(1+psi^2T^2+2psi(2+3T+2T^2))).Moment; Veb=subs(Veb_sym,'t',tau);

plot(Veb,'g'); xlabel('time'); ylabel('voltage');