I'm suppose to plot the response of the system but i'm getting wrong values
1 次查看(过去 30 天)
显示 更早的评论
i was trying to solve the system , but i'm getting wrong va;ues , what am i missing?here is my code
clc;
clear all;
close all;
% Define system parameters
params = struct();
% System parameters
params.M1 = 6e-3; % tip Mass of beam I in kg
params.M2 = 5e-3; % tip Mass of beam II in kg
params.L1 = 7e-2; % Length of beam I in meters
params.L2 = 4e-2; % Length of beam II in meters
params.w1 = 2e-2; % Width of beam I in meters
params.w2 = 2e-2; % Width of beam II in meters
params.h1 = 1e-3; % Thickness of beam I in meters
params.h2 = 1e-3; % Thickness of beam II in meters
params.E1 = 69e9; % Young's modulus of beam I in Pascals
params.E2 = 2.344e9; % Young's modulus of beam II in Pascals
params.c1 = 0.1; % Damping of beam I in Ns/m
params.c2 = 0.1; % Damping of beam II in Ns/m
params.c12 = 0; % Couplind damping in Ns/m
params.c21 = 0; % Couplind damping in Ns/m
params.theta = 0.4e-3; % Piezoelectric constant in m/V
params.R = 10e6; % Resistance in ohms
params.Cp = 72e-9; % Piezoelectric capacity in farads
params.F0 = 5; % Amplitude of base excitation force in Newton
params.omega_a = 10; % Angular frequency of base excitation in rad/s
params.epsilon_0 = 8.854e-12; % Vacuum permittivity in F/m
params.tau = 1e-3; % Time constant in seconds
params.epsilon_r = 1; % Relative permittivity
params.do_2 = 1e-3; % Constant
params.sigma = 1e-9; % Conductivity in S/m
params.rho1 = 2700; % Density for Beam I in kg/m^3
params.rho2 = 1220; % Density for Beam II in kg/m^3
params.lt = 1e-2; % the length of the Tribo electric material in m
params.wt = 1e-2; % the width of the Tribo electric material in m
params.S = params.lt * params.wt; % contact area of the tribo material
params.gi = 6e-1; % the gap in m
params.I1 = (1/12) * params.w1 * (params.h1^3); % moment of inertia for first beam
params.I2 = (1/12) * params.w2 * (params.h2^3); % moment of inertia for second beam
params.alpha = 2*params.epsilon_0*params.epsilon_r*params.S;
params.alpha2 = params.epsilon_0*params.R*params.S;
params.pt = 2e-2 ; %thickness of pizo
params.sigma = 1e-6; %surfce charge density for PDMS
params.gap = 6e-3; % the gap For PDMS
params.alpha3 = ((params.pt/params.epsilon_0)+params.do_2-params.gap)+((params.sigma/(params.R*params.epsilon_0))*(params.do_2-params.gap));
% Calculate equivalent masses
params.m1 = params.rho1 * params.L1 * params.w1 * params.h1;
params.m2 = params.rho2 * params.L2 * params.w2 * params.h2;
params.Meq1 = params.M1 + (33 * params.m1 * params.L1) / 140;
params.Meq2 = params.M2 + (33 * params.m2 * params.L2) / 140;
% Calculate ko
ko = 6 * params.E1 * params.I1 / ((4 * params.E1 * params.I1 * params.L1^3 * params.L2^3) + (3 * params.E2 * params.I2 * params.L1^4 * params.L2^2));
% Calculate K matrix
K = ko * [2 * ((2 * params.E1 * params.I1 * params.L2^3) + 6 * params.E2 * params.I2 * (params.L1 * params.L2^2 - params.L1^2 * params.L2) + (2 * params.E2 * params.I2 * params.L1^3)), ...
(3 * params.E2 * params.I2 * params.L1^2 * params.L2) - (2 * params.E2 * params.I2 * params.L1^3);
(3 * params.E2 * params.I2 * params.L1^2 * params.L2) - (2 * params.E2 * params.I2 * params.L1^3), ...
2 * params.E2 * params.I2 * params.L1^3];
% Extract individual elements
params.k11 = K(1, 1);
params.k12 = K(1, 2);
params.k21 = K(2, 1);
params.k22 = K(2, 2);
% Define state-space matrices
A = [0, 1, 0, 0, 0, 0;
-(params.k11/params.Meq1), -(params.c1/params.Meq1), -(params.k12/params.Meq1), 0, 0, -(params.theta/params.Meq1);
0, 0, 0, 1, 0, 0;
-(params.k21/params.Meq2), 0, -(params.k22/params.Meq2), -(params.c21/params.Meq2), 0, 0;
0, 0, 0, 0, 0, 1/params.R;
0, 0, 0, 0, 0, -(params.theta/params.Cp)];
% Define time span
tspan = [0 20]; % Time span from 0 to 10 seconds
% Define initial conditions
x0 = [1;1;1;1;1;1];
% Define input function (if any)
a = @(t) params.F0*sin(params.omega_a*t); % No input function for now
% Define the function for the system of first-order differential equations
dxdt = @(t, x) [x(2);
-(params.c1/params.Meq1)*x(2) - (params.c12/params.Meq1)*x(4) - (params.k11/params.Meq1)*x(1) - (params.k12/params.Meq1)*x(3) - (params.theta/params.Meq1)*x(6) - a(t);
x(4);
-(params.c2/params.Meq2)*x(4) - (params.c21/params.Meq2)*x(2) - (params.k22/params.Meq2)*x(3) - (params.k21/params.Meq2)*x(1) - ((x(5))^2)/(params.alpha)- a(t);
((x(5))/params.alpha2) * (params.alpha3 + params.do_2 - params.gap) + params.alpha3 * (params.do_2 - params.gap);
x(6)/params.R + params.Cp * x(6)];
% Solve the system using ode45
[t, x] = ode45(dxdt, tspan, x0);
% Extracting the results
y1 = x(:,1);
y1dot = x(:,2);
y2 = x(:,3);
y2dot = x(:,4);
q = x(:,5);
v = x(:,6);
% Plot the results if needed
% For example, plot y1 and y2 over time
plot(t, q, 'b', t, v, 'r');
xlabel('Time');
ylabel('Displacement');
legend('y1', 'y2');
title('Displacement Response');
This the equations :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1620683/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1620688/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1620693/image.png)
2 个评论
Torsten
2024-2-18
编辑:Torsten
2024-2-18
Why does params.gap equal x(3) in equation (5) ?
And note that in equation (5), you divide by params.alpha2 which is in the order of 1e-12. This will blow up the time derivative tremendously.
And you will have to define your dxdt in a function instead of a function handle because the differential equation for y2 depends on a condition.
回答(1 个)
Torsten
2024-2-19
编辑:Torsten
2024-2-19
alpha is a constant , and this is the value of it ,
I doubt that the rate of change at t=0 for x(5) should be 2.55e17 as you get when you call dxdt for the vector of initial values (see below). The main contribution is a factor of 1e12 caused by params.alpha2.
params.gap = 6e-6 not x3
The equation for q2 in your mathematical description is written with y2, not with params.gap.
and for dxdt , what do you mean by a function , you mean seperate one ? if so , can you please guide me through the process of definning it and call to the ode45
By defining dxdt in a function file as below, you can easily use an if-condition to define the equation for y2 state-dependent:
% Define time span
tspan = [0 20]; % Time span from 0 to 10 seconds
% Define initial conditions
x0 = [1;1;1;1;1;1];
% Solve the system using ode45
fun(0,x0)
[t, x] = ode45(@fun, tspan, x0);
% Extracting the results
y1 = x(:,1);
y1dot = x(:,2);
y2 = x(:,3);
y2dot = x(:,4);
q = x(:,5);
v = x(:,6);
% Plot the results if needed
% For example, plot y1 and y2 over time
plot(t, q, 'b', t, v, 'r');
xlabel('Time');
ylabel('Displacement');
legend('y1', 'y2');
title('Displacement Response');
function dxdt = fun(t,x)
% Define system parameters
params = struct();
% System parameters
params.M1 = 6e-3; % tip Mass of beam I in kg
params.M2 = 5e-3; % tip Mass of beam II in kg
params.L1 = 7e-2; % Length of beam I in meters
params.L2 = 4e-2; % Length of beam II in meters
params.w1 = 2e-2; % Width of beam I in meters
params.w2 = 2e-2; % Width of beam II in meters
params.h1 = 1e-3; % Thickness of beam I in meters
params.h2 = 1e-3; % Thickness of beam II in meters
params.E1 = 69e9; % Young's modulus of beam I in Pascals
params.E2 = 2.344e9; % Young's modulus of beam II in Pascals
params.c1 = 0.1; % Damping of beam I in Ns/m
params.c2 = 0.1; % Damping of beam II in Ns/m
params.c12 = 0; % Couplind damping in Ns/m
params.c21 = 0; % Couplind damping in Ns/m
params.theta = 0.4e-3; % Piezoelectric constant in m/V
params.R = 10e6; % Resistance in ohms
params.Cp = 72e-9; % Piezoelectric capacity in farads
params.F0 = 5; % Amplitude of base excitation force in Newton
params.omega_a = 10; % Angular frequency of base excitation in rad/s
params.epsilon_0 = 8.854e-12; % Vacuum permittivity in F/m
params.tau = 1e-3; % Time constant in seconds
params.epsilon_r = 1; % Relative permittivity
params.do_2 = 1e-3; % Constant
params.sigma = 1e-9; % Conductivity in S/m
params.rho1 = 2700; % Density for Beam I in kg/m^3
params.rho2 = 1220; % Density for Beam II in kg/m^3
params.lt = 1e-2; % the length of the Tribo electric material in m
params.wt = 1e-2; % the width of the Tribo electric material in m
params.S = params.lt * params.wt; % contact area of the tribo material
params.gi = 6e-1; % the gap in m
params.I1 = (1/12) * params.w1 * (params.h1^3); % moment of inertia for first beam
params.I2 = (1/12) * params.w2 * (params.h2^3); % moment of inertia for second beam
params.alpha = 2*params.epsilon_0*params.epsilon_r*params.S;
params.alpha2 = params.epsilon_0*params.R*params.S;
params.pt = 2e-2 ; %thickness of pizo
params.sigma = 1e-6; %surfce charge density for PDMS
params.gap = 6e-3; % the gap For PDMS
params.alpha3 = ((params.pt/params.epsilon_0)+params.do_2-params.gap)+((params.sigma/(params.R*params.epsilon_0))*(params.do_2-params.gap));
% Calculate equivalent masses
params.m1 = params.rho1 * params.L1 * params.w1 * params.h1;
params.m2 = params.rho2 * params.L2 * params.w2 * params.h2;
params.Meq1 = params.M1 + (33 * params.m1 * params.L1) / 140;
params.Meq2 = params.M2 + (33 * params.m2 * params.L2) / 140;
% Calculate ko
ko = 6 * params.E1 * params.I1 / ((4 * params.E1 * params.I1 * params.L1^3 * params.L2^3) + (3 * params.E2 * params.I2 * params.L1^4 * params.L2^2));
% Calculate K matrix
K = ko * [2 * ((2 * params.E1 * params.I1 * params.L2^3) + 6 * params.E2 * params.I2 * (params.L1 * params.L2^2 - params.L1^2 * params.L2) + (2 * params.E2 * params.I2 * params.L1^3)), ...
(3 * params.E2 * params.I2 * params.L1^2 * params.L2) - (2 * params.E2 * params.I2 * params.L1^3);
(3 * params.E2 * params.I2 * params.L1^2 * params.L2) - (2 * params.E2 * params.I2 * params.L1^3), ...
2 * params.E2 * params.I2 * params.L1^3];
% Extract individual elements
params.k11 = K(1, 1);
params.k12 = K(1, 2);
params.k21 = K(2, 1);
params.k22 = K(2, 2);
% Define state-space matrices
A = [0, 1, 0, 0, 0, 0;
-(params.k11/params.Meq1), -(params.c1/params.Meq1), -(params.k12/params.Meq1), 0, 0, -(params.theta/params.Meq1);
0, 0, 0, 1, 0, 0;
-(params.k21/params.Meq2), 0, -(params.k22/params.Meq2), -(params.c21/params.Meq2), 0, 0;
0, 0, 0, 0, 0, 1/params.R;
0, 0, 0, 0, 0, -(params.theta/params.Cp)];
% Define input function (if any)
a = params.F0*sin(params.omega_a*t); % No input function for now
% Define the function for the system of first-order differential equations
dxdt = [x(2);
-(params.c1/params.Meq1)*x(2) - (params.c12/params.Meq1)*x(4) - (params.k11/params.Meq1)*x(1) - (params.k12/params.Meq1)*x(3) - (params.theta/params.Meq1)*x(6) - a;
x(4);
-(params.c2/params.Meq2)*x(4) - (params.c21/params.Meq2)*x(2) - (params.k22/params.Meq2)*x(3) - (params.k21/params.Meq2)*x(1) - ((x(5))^2)/(params.alpha)- a;
((x(5))/params.alpha2) * (params.alpha3 + params.do_2 - params.gap) + params.alpha3 * (params.do_2 - params.gap);
x(6)/params.R + params.Cp * x(6)];
end
4 个评论
Sam Chak
2024-2-19
@Mohammad Adeeb, I noticed that the state matrix
in your code is unused. Could you please clarify this? Additionally, have you verified if the differential equations are correctly written in the code? If possible, it would be helpful if you could provide an image of the expected graph for
and
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1621198/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1621203/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1621208/image.png)
Sam Chak
2024-2-19
I have some concerns about the reliability of the differential equations presented in your image. Typically, reputable journals or publishers adhere to standard mathematical notations. It appears that the equations might have been typed using Microsoft Equation Editor or LaTeX without being thoroughly proofread by professionals.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1621278/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1621283/image.png)
If the differential equations in your image are indeed correct, they can be rearranged as follows:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1621288/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1621293/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1621298/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1621303/image.png)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!