Discussion Closed This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.

error in output

Please login with a confirmed email address before reporting spam

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');


0 Replies Last Post 4 feb 2018, 14:00 GMT-5
COMSOL Moderator

Hello Mangai Padma

Your Discussion has gone 30 days without a reply. If you still need help with COMSOL and have an on-subscription license, please visit our Support Center for help.

If you do not hold an on-subscription license, you may find an answer in another Discussion or in the Knowledge Base.

Note that while COMSOL employees may participate in the discussion forum, COMSOL® software users who are on-subscription should submit their questions via the Support Center for a more comprehensive response from the Technical Support team.