Energy Harvesting of Cantilever Beam with MSMA alloy1

2 views (last 30 days)
Hamed
Hamed on 19 Feb 2025
Edited: Torsten on 19 Feb 2025
Ran in:
There is an error for line 189 (M_ww) and I think after solve this error, the same error will be for line 193 (Fbar). The error is:
Error using subs
Expected input number 1, S, to be one of these types:
sym
Instead its type was function_handle.
Error in sym/subs (line 69)
validateattributes(F, {'sym'}, {}, 'subs', 'S', 1);
Could you please help me? Thank you
clc; clear;
% MSMA and Material Parameters
epsilon_r_max = 6.1 / 100; % Maximum strain from martensitic rearrangement (6.1%)
epsilon_0 = 3 / 100; % Pre-strain (3%)
M_sat = 574360;%0.72; % Saturation magnetization (T)
ro_K1 = 1.9e5; % Magnetic crystal anisotropy energy (J/m^3)
ro_M = 7900; % Density of the MSMA (kg/m^3)
mio_0 = 4 * pi * 1e-7; % Permeability of free space (H/m)
D_xx = 0.04; % Demagnetization factor in x-direction
D_zz = 0.48; % Demagnetization factor in z-direction
L_M = 0.02; % Length of the MSMA (m)
E_M = 9e9; % Elastic modulus of the MSMA (Pa)
w_M = 0.003; % Width of the MSMA (m)
h_M = 0.003; % Thickness of the MSMA (m)
A_M = w_M * h_M; % MSMA Area (m^2)
H_app = 636619;%0.65; % Applied external magnetic field (A/m)
sigma_s12 = -1.9e6; % the stress value at the start points of variant reorientation from 1 to 2 (Pa)
sigma_f12 = -1.4e6; % the stress value at the finish points of variant reorientation from 1 to 2 (Pa)
sigma_s21 = -3e6; % the stress value at the start points of variant reorientation from 2 to 1 (Pa)
sigma_f21 = -4e6; % the stress value at the finish points of variant reorientation from 2 to 1 (Pa)
S = 1 / 5.291e-11; % effective elastic compliance (Pa)
% Beam Properties
ro_B = 2700; % Density of the beam (kg/m^3)
E_B = 70e9; % Elastic modulus of the beam (Pa)
L_B = 0.1; % Length of the beam (m)
w_B = 0.01; % Width of the beam (m)
h_B = 0.0012; % Thickness of the beam (m)
A_B = w_B * h_B; % Beam Area (m^2)
d = 0.0055; % Step distance for MSMA connection (m)
% Proof Mass Properties
M_p = 8e-3; % Proof mass at beam tip (kg)
I_p = 4e-4; %(1/12) * M_p * (w_B^2 + h_B^2); % Moment of inertia of proof mass (kg·m²)
% Electrical Parameters for Coil
L = 8.24e-3; % Inductance of the coil (H)
R = 500; % Resistance (Ohms)
C = 6.1e-3; % Capacitance (F)
N_coil = 2500; % Number of coil turns
A_coil = 0.0000025; % Area of the coil (m^2)
% Angles (Converted to Radians)
teta_3_0 = deg2rad(54.54); % Rotation angle for zeta = 0
teta_3_1 = deg2rad(43.27); % Rotation angle for zeta = 1
teta_4_0 = deg2rad(0.807); % Rotation angle for zeta = 0
teta_4_1 = deg2rad(0); % Rotation angle for zeta = 1
% Magnetic Fields
H_z_0 = 412040;%0.52; % Magnetic field in z-direction when zeta = 0 (T)
H_z_1 = 360920;%0.45; % Magnetic field in z-direction when zeta = 1 (T)
H_x_0 = -13326;%-0.017; % Magnetic field in x-direction when zeta = 0 (T)
H_x_1 = 0; % Magnetic field in x-direction when zeta = 1 (T)
A_hyst = 14840; % Hysteresis constant A (Pa)
C_hyst = 45340; % Hysteresis constant C (Pa)
B1_hyst = 6913; % Hysteresis constant B1 (Pa)
B2_hyst = 7625; % Hysteresis constant B2 (Pa)
Y_hyst = 56430; % Threshold stress for phase transformation (Pa)
% Define time vector
t = linspace(0, 0.6, 100); % Time from 0 to 1 sec, with 100 points
% Force parameters
F0 = 10; % Magnitude of impulse force (N)
t1 = 0.1; % Start time of impulse
t2 = 0.5; % End time of impulse
F = F0 * (t >= t1 & t <= t2);
% Plot the force
figure;
plot(t, F, 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Force (N)');
title('Rectangular Impulse Force');
grid on;
% Zeta calculations
zeta_s12 = (1/A_hyst) * ((sigma_s12 * epsilon_r_max) - (mio_0 * M_sat * H_x_0 * (cos(teta_3_0) + sin(teta_4_0))) + (mio_0 * M_sat * H_z_0 * (cos(teta_4_0) - sin(teta_3_0))) + (ro_K1 * (sin(teta_3_0)^2 - sin(teta_4_0)^2)) - B1_hyst - B2_hyst - Y_hyst);
zeta_f12 = (1/A_hyst) * ((sigma_f12 * epsilon_r_max) - (mio_0 * M_sat * H_x_1 * (cos(teta_3_1) + sin(teta_4_1))) + (mio_0 * M_sat * H_z_1 * (cos(teta_4_1) - sin(teta_3_1))) + (ro_K1 * (sin(teta_3_1)^2 - sin(teta_4_1)^2)) - B1_hyst - B2_hyst - Y_hyst);
zeta_s21 = (1/C_hyst) * ((sigma_s21 * epsilon_r_max) - (mio_0 * M_sat * H_x_1 * (cos(teta_3_1) - sin(teta_4_1))) + (mio_0 * M_sat * H_z_1 * (cos(teta_4_1) - sin(teta_3_1))) + (ro_K1 * (sin(teta_3_1)^2 - sin(teta_4_1)^2)) - B1_hyst + B2_hyst + Y_hyst);
zeta_f21 = (1/C_hyst) * ((sigma_f21 * epsilon_r_max) - (mio_0 * M_sat * H_x_0 * (cos(teta_3_0) - sin(teta_4_0))) + (mio_0 * M_sat * H_z_0 * (cos(teta_4_0) - sin(teta_3_0))) + (ro_K1 * (sin(teta_3_0)^2 - sin(teta_4_0)^2)) - B1_hyst + B2_hyst + Y_hyst);
% Parameters
I0_B = ro_B * A_B;
I1_B = (ro_B * A_B ^2) / 2;
I0_M = ro_M * A_M;
I1_M = (ro_M * A_M ^2) / 2;
k0_B = E_B * A_B;
k1_B = (E_B * A_B ^2) / 2;
k2_B = (E_B * A_B ^3) / 3;
k0_M = E_M * A_M;
k1_M = (E_M * A_M ^2) / 2;
k2_M = (E_M * A_M ^3) / 3;
syms x beta_u beta_w
eq1 = ((k0_B)*(M_p)*(beta_u)*(sin(beta_u))*L_B)-((I0_B)*(cos(beta_u))*L_B) == 0;
eq2 = 1 + (cos(beta_w)*L_B*cosh(beta_w)*L_B) + (((beta_w)*(M_p)/(I0_B))*((cos(beta_w)*L_B*sinh(beta_w)*L_B) - (sin(beta_w)*L_B*cosh(beta_w)*L_B))) - ...
((((beta_w)^3)*(I_p)/(I0_B))*((cosh(beta_w)*L_B*sin(beta_w)*L_B) + (sinh(beta_w)*L_B*cos(beta_w)*L_B))) + ...
((((beta_w)^4)*(M_p)*(I_p)/((I0_B)^2))*(1 - ((cos(beta_w)*L_B*L_B*cosh(beta_w))))) == 0;
[beta_u, beta_w] = vpasolve([eq1, eq2], [beta_u, beta_w]);
% Display results
fprintf('Solution for beta_u: %f\n', beta_u);
Solution for beta_u: -226.194671
fprintf('Solution for beta_w: %f\n', beta_w);
Solution for beta_w: 98.999533
% Define x array
N_points = 100;
%x = linspace(0, L_B, N_points);
% Define the functions
phiu = matlabFunction(sin(beta_u * x));
phiw = matlabFunction((cos(beta_w * x)) - (cosh(beta_w) * x) + (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B))))) * ((sin(beta_w * x)) - (sinh(beta_w * x))));
phiu2 = matlabFunction((sin(beta_u * x))^2);
phiw2 = matlabFunction(((cos(beta_w * x)) - (cosh(beta_w) * x) + (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B))))) * ((sin(beta_w * x)) - (sinh(beta_w * x))))^2);
% Compute derivatives
phiu_x = matlabFunction(beta_u * cos(beta_u * x));
phiu_xx = matlabFunction((-beta_u^2) * sin(beta_u * x));
phiu_x2 = matlabFunction((beta_u^2 * cos(beta_u * x)^2));
phiw_x = matlabFunction((-beta_w * sin(beta_w * x)) + ((beta_w * cos(beta_w * x) - beta_w * cosh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - (beta_w * sinh(beta_w * x)));
phiw_xx = matlabFunction(((-beta_w ^2) * cos(beta_w * x)) + (((-beta_w ^2) * sin(beta_w * x) - (beta_w ^2) * sinh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - ((beta_w ^2) * cosh(beta_w * x)));
phiw_x2 = matlabFunction(((-beta_w * sin(beta_w * x)) + ((beta_w * cos(beta_w * x) - beta_w * cosh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - (beta_w * sinh(beta_w * x)))^2);
phiw_x4 = matlabFunction(((-beta_w * sin(beta_w * x)) + ((beta_w * cos(beta_w * x) - beta_w * cosh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - (beta_w * sinh(beta_w * x)))^4);
phiw_xx2 = matlabFunction((((-beta_w ^2) * cos(beta_w * x)) + (((-beta_w ^2) * sin(beta_w * x) - (beta_w ^2) * sinh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - ((beta_w ^2) * cosh(beta_w * x)))^2);
% Define the integrals
% For K_uu
integral_1 = k0_B * integral(phiu_x2, 0, L_B);
integral_2 = k0_M * integral(phiu_x2, 0, L_M);
% For K_uw
product_phiu_x_phiw_xx = @(x) phiu_x(x) .* phiw_xx(x);
integral_3 = k1_B * integral(product_phiu_x_phiw_xx,0,L_B);
integral_4 = k1_M * integral(product_phiu_x_phiw_xx,0,L_M);
% For K_ww
integral_5 = k2_B * integral(phiw_xx2, 0, L_B);
integral_6 = k2_M * integral(phiw_xx2, 0, L_M);
% For K_uww
product_phiu_x_phiw_x2 = @(x) phiu_x(x) .* phiw_x2(x);
integral_7 = k0_B * integral(product_phiu_x_phiw_x2,0,L_B);
integral_8 = k0_M * integral(product_phiu_x_phiw_x2,0,L_M);
% For K_www
product_phiu_x2_phiw_xx = @(x) phiu_x2(x) .* phiw_xx(x);
integral_9 = k1_B * integral(product_phiu_x2_phiw_xx,0,L_B);
integral_10 = k1_M * integral(product_phiu_x2_phiw_xx,0,L_M);
% For K_wwww
integral_11 = k0_B * integral(phiw_x4, 0, L_B);
integral_12 = k0_M * integral(phiw_x4, 0, L_M);
% For M_uu
integral_13 = I0_B * integral(phiu2, 0, L_B);
integral_14 = I0_M * integral(phiu2, 0, L_M);
% For M_uw
product_phiu_phiw_x = @(x) phiu(x) .* phiw_x(x);
integral_15 = I1_B * integral(product_phiu_phiw_x,0,L_B);
integral_16 = I1_M * integral(product_phiu_phiw_x,0,L_M);
% For M_ww
integral_17 = I0_B * integral(phiw2, 0, L_B);
integral_18 = I0_M * integral(phiw2, 0, L_B);
% For Kbar_u
integral_19 = integral(phiu_x, 0, L_M);
% For Kbar_w
integral_20 = integral(phiw_xx, 0, L_M);
% For Kbar_ww
integral_21 = integral(phiw_x, 0, L_M);
% Stiffness and Mass Matrix Calculation
K_uu = integral_1 + integral_2;
K_uw = -2 .* (integral_3 + integral_4);
K_ww = integral_5 + integral_6;
K_uww = integral_7 + integral_8;
K_www = -integral_9 - integral_10;
K_wwww = 1/4 .* (integral_11 + integral_12);
M_uu = integral_13 + integral_14 + (M_p * L_B^2 * sin(beta_u)^2);
M_uw = -2 .* (integral_15 + integral_16);
M_ww = integral_17 + integral_18 + (M_p * (subs(phiw2, x, L_B))) + (I_p * (subs(phiw_x2, x, L_B)));
Error using subs
Expected input number 1, S, to be one of these types:

sym

Instead its type was function_handle.

Error in sym/subs (line 69)
validateattributes(F, {'sym'}, {}, 'subs', 'S', 1);
Kbar_u =k0_M * (epsilon_0 - ((1 - zeta_f12) * epsilon_r_max)) * (integral_19);
Kbar_w =-k1_M * (epsilon_0 - ((1 - zeta_f12) * epsilon_r_max)) * (integral_20);
Kbar_ww =k0_M * (epsilon_0 - ((1 - zeta_f12) * epsilon_r_max)) * (integral_21);
Fbar = F * subs(phiw, x, L_B);
% Time span for simulation
tspan = linspace(0, 0.01, 50); % Adjust as needed
% Initial conditions [q_u(0), dq_u/dt(0), q_w(0), dq_w/dt(0)]
y0 = [0; 0; 0; 0]; % Replace with actual values
opts = odeset('RelTol',1e-6,'AbsTol',1e-8);
% Solve ODEs using ode45
[t, Y] = ode45(@(t, y) system_eqs(t, y, M_uu, M_ww, M_uw, ...
K_uu, K_ww, K_uw, K_uww, K_wwww, ...
Kbar_u, Kbar_w, Kbar_ww, Fbar), tspan, y0);
% Extract solutions
q_u = Y(:, 1);
q_w = Y(:, 3);
% Plot results
figure;
plot(t, q_u, 'b-', 'LineWidth', 2); hold on;
plot(t, q_w, 'r--', 'LineWidth', 2);
xlabel('Time (t)');
ylabel('q_u, q_w');
legend('q_u(t)', 'q_w(t)');
title('Solution of q_u and q_w over time');
grid on;
% Compute Displacements
[X, T] = meshgrid(x, t);
U = sin(beta_u * X) .* interp1(t, q_u, T, 'linear', 'extrap');
W = (cos(beta_w * X) - cosh(beta_w * X)) .* interp1(t, q_w, T, 'linear', 'extrap');
% Plot Results
figure;
plot(t, q_u, 'b-', 'LineWidth', 2); hold on;
plot(t, q_w, 'r--', 'LineWidth', 2);
xlabel('Time (t)'); ylabel('q_u, q_w');
legend('q_u(t)', 'q_w(t)');
title('Solution of q_u and q_w over time'); grid on;
figure;
surf(T, X, U, 'EdgeColor', 'none');
xlabel('Time (t)'); ylabel('Position (x)'); zlabel('u(x,t)');
title('Longitudinal Displacement u(x,t)');
colorbar; view(3);
figure;
surf(T, X, W, 'EdgeColor', 'none');
xlabel('Time (t)'); ylabel('Position (x)'); zlabel('w(x,t)');
title('Transverse Displacement w(x,t)');
colorbar; view(3);
% Compute voltage using numerical differentiation
V = -gradient((mio_0 * M_sat * (1 - zeta) * cos(teta_3_0) - zeta * sin(teta_4_0)), tspan) * N_coil * A_coil;
% Plot voltage over time
figure;
plot(t, V, 'r', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Voltage (V)');
title('Induced Voltage Over Time');
grid on;
P = abs((V^2) / R);
subplot(4, 1, 3);
plot(t, P);
title('Power Output');
xlabel('Time (s)'); ylabel('Power (W)');
% Function defining the system of equations
function dydt = system_eqs(~, y, M_uu, M_ww, M_uw, ...
K_uu, K_ww, K_uw, K_uww, K_wwww, ...
Kbar_u, Kbar_w, Kbar_ww, Fbar)
% Extract variables
q_u = y(1); dq_u = y(2);
q_w = y(3); dq_w = y(4);
% Define mass matrix
M = [M_uu, 0.5*M_uw;
0.5*M_uw, M_ww];
% Define stiffness and nonlinear terms
K1 = [K_uu*q_u + 0.5*K_uw*q_w + K_uww*q_w^2 - Kbar_u;
K_ww*q_w + 0.5*K_uw*q_u + K_uww*q_u*q_w + 2*K_wwww*q_w^3 - (Kbar_w + Kbar_ww*q_w + Fbar)];
% Solve for accelerations
accel = M \ (-K1);
% Define dydt
dydt = [dq_u; accel(1); dq_w; accel(2)];
end

Sign in to answer this question.

Accepted Answer

Torsten
Torsten on 19 Feb 2025
Edited: Torsten on 19 Feb 2025
Ran in:
I inserted "warning('off')" at the start of your code because your matrix M in your function "system_eqs" is singular and MATLAB complains about it.
I made a lot of modifications to your code - so check whether the lines marked with a "%changed" make sense to you.
As long as the warnings appear and "ode45" and "ode15s" give different results, I wouldn't trust in them.
I also lost motivation to care about the final error in the graphical representation of the results - you should be able to manage this part on your own.
warning('off') %changed
clc; clear;
% MSMA and Material Parameters
epsilon_r_max = 6.1 / 100; % Maximum strain from martensitic rearrangement (6.1%)
epsilon_0 = 3 / 100; % Pre-strain (3%)
M_sat = 574360;%0.72; % Saturation magnetization (T)
ro_K1 = 1.9e5; % Magnetic crystal anisotropy energy (J/m^3)
ro_M = 7900; % Density of the MSMA (kg/m^3)
mio_0 = 4 * pi * 1e-7; % Permeability of free space (H/m)
D_xx = 0.04; % Demagnetization factor in x-direction
D_zz = 0.48; % Demagnetization factor in z-direction
L_M = 0.02; % Length of the MSMA (m)
E_M = 9e9; % Elastic modulus of the MSMA (Pa)
w_M = 0.003; % Width of the MSMA (m)
h_M = 0.003; % Thickness of the MSMA (m)
A_M = w_M * h_M; % MSMA Area (m^2)
H_app = 636619;%0.65; % Applied external magnetic field (A/m)
sigma_s12 = -1.9e6; % the stress value at the start points of variant reorientation from 1 to 2 (Pa)
sigma_f12 = -1.4e6; % the stress value at the finish points of variant reorientation from 1 to 2 (Pa)
sigma_s21 = -3e6; % the stress value at the start points of variant reorientation from 2 to 1 (Pa)
sigma_f21 = -4e6; % the stress value at the finish points of variant reorientation from 2 to 1 (Pa)
S = 1 / 5.291e-11; % effective elastic compliance (Pa)
% Beam Properties
ro_B = 2700; % Density of the beam (kg/m^3)
E_B = 70e9; % Elastic modulus of the beam (Pa)
L_B = 0.1; % Length of the beam (m)
w_B = 0.01; % Width of the beam (m)
h_B = 0.0012; % Thickness of the beam (m)
A_B = w_B * h_B; % Beam Area (m^2)
d = 0.0055; % Step distance for MSMA connection (m)
% Proof Mass Properties
M_p = 8e-3; % Proof mass at beam tip (kg)
I_p = 4e-4; %(1/12) * M_p * (w_B^2 + h_B^2); % Moment of inertia of proof mass (kg·m²)
% Electrical Parameters for Coil
L = 8.24e-3; % Inductance of the coil (H)
R = 500; % Resistance (Ohms)
C = 6.1e-3; % Capacitance (F)
N_coil = 2500; % Number of coil turns
A_coil = 0.0000025; % Area of the coil (m^2)
% Angles (Converted to Radians)
teta_3_0 = deg2rad(54.54); % Rotation angle for zeta = 0
teta_3_1 = deg2rad(43.27); % Rotation angle for zeta = 1
teta_4_0 = deg2rad(0.807); % Rotation angle for zeta = 0
teta_4_1 = deg2rad(0); % Rotation angle for zeta = 1
% Magnetic Fields
H_z_0 = 412040;%0.52; % Magnetic field in z-direction when zeta = 0 (T)
H_z_1 = 360920;%0.45; % Magnetic field in z-direction when zeta = 1 (T)
H_x_0 = -13326;%-0.017; % Magnetic field in x-direction when zeta = 0 (T)
H_x_1 = 0; % Magnetic field in x-direction when zeta = 1 (T)
A_hyst = 14840; % Hysteresis constant A (Pa)
C_hyst = 45340; % Hysteresis constant C (Pa)
B1_hyst = 6913; % Hysteresis constant B1 (Pa)
B2_hyst = 7625; % Hysteresis constant B2 (Pa)
Y_hyst = 56430; % Threshold stress for phase transformation (Pa)
% Define time vector
T = linspace(0, 0.6, 100); % Time from 0 to 1 sec, with 100 points % changed
% Force parameters
F0 = 10; % Magnitude of impulse force (N)
t1 = 0.1; % Start time of impulse
t2 = 0.5; % End time of impulse
F = @(t)F0 * (t >= t1 & t <= t2); %changed
% Plot the force
figure;
plot(T, F(T), 'LineWidth', 2); %changed
xlabel('Time (s)');
ylabel('Force (N)');
title('Rectangular Impulse Force');
grid on;
% Zeta calculations
zeta_s12 = (1/A_hyst) * ((sigma_s12 * epsilon_r_max) - (mio_0 * M_sat * H_x_0 * (cos(teta_3_0) + sin(teta_4_0))) + (mio_0 * M_sat * H_z_0 * (cos(teta_4_0) - sin(teta_3_0))) + (ro_K1 * (sin(teta_3_0)^2 - sin(teta_4_0)^2)) - B1_hyst - B2_hyst - Y_hyst);
zeta_f12 = (1/A_hyst) * ((sigma_f12 * epsilon_r_max) - (mio_0 * M_sat * H_x_1 * (cos(teta_3_1) + sin(teta_4_1))) + (mio_0 * M_sat * H_z_1 * (cos(teta_4_1) - sin(teta_3_1))) + (ro_K1 * (sin(teta_3_1)^2 - sin(teta_4_1)^2)) - B1_hyst - B2_hyst - Y_hyst);
zeta_s21 = (1/C_hyst) * ((sigma_s21 * epsilon_r_max) - (mio_0 * M_sat * H_x_1 * (cos(teta_3_1) - sin(teta_4_1))) + (mio_0 * M_sat * H_z_1 * (cos(teta_4_1) - sin(teta_3_1))) + (ro_K1 * (sin(teta_3_1)^2 - sin(teta_4_1)^2)) - B1_hyst + B2_hyst + Y_hyst);
zeta_f21 = (1/C_hyst) * ((sigma_f21 * epsilon_r_max) - (mio_0 * M_sat * H_x_0 * (cos(teta_3_0) - sin(teta_4_0))) + (mio_0 * M_sat * H_z_0 * (cos(teta_4_0) - sin(teta_3_0))) + (ro_K1 * (sin(teta_3_0)^2 - sin(teta_4_0)^2)) - B1_hyst + B2_hyst + Y_hyst);
% Parameters
I0_B = ro_B * A_B;
I1_B = (ro_B * A_B ^2) / 2;
I0_M = ro_M * A_M;
I1_M = (ro_M * A_M ^2) / 2;
k0_B = E_B * A_B;
k1_B = (E_B * A_B ^2) / 2;
k2_B = (E_B * A_B ^3) / 3;
k0_M = E_M * A_M;
k1_M = (E_M * A_M ^2) / 2;
k2_M = (E_M * A_M ^3) / 3;
syms x beta_u beta_w
eq1 = ((k0_B)*(M_p)*(beta_u)*(sin(beta_u))*L_B)-((I0_B)*(cos(beta_u))*L_B) == 0;
eq2 = 1 + (cos(beta_w)*L_B*cosh(beta_w)*L_B) + (((beta_w)*(M_p)/(I0_B))*((cos(beta_w)*L_B*sinh(beta_w)*L_B) - (sin(beta_w)*L_B*cosh(beta_w)*L_B))) - ...
((((beta_w)^3)*(I_p)/(I0_B))*((cosh(beta_w)*L_B*sin(beta_w)*L_B) + (sinh(beta_w)*L_B*cos(beta_w)*L_B))) + ...
((((beta_w)^4)*(M_p)*(I_p)/((I0_B)^2))*(1 - ((cos(beta_w)*L_B*L_B*cosh(beta_w))))) == 0;
[beta_u, beta_w] = vpasolve([eq1, eq2], [beta_u, beta_w]);
beta_u = double(beta_u); %changed
beta_w = double(beta_w); %changed
% Display results
fprintf('Solution for beta_u: %f\n', beta_u);
Solution for beta_u: -226.194671
fprintf('Solution for beta_w: %f\n', beta_w);
Solution for beta_w: 98.999533
% Define x array
N_points = 100;
%x = linspace(0, L_B, N_points);
% Define the functions
phiu = matlabFunction(sin(beta_u * x));
phiw = matlabFunction((cos(beta_w * x)) - (cosh(beta_w) * x) + (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B))))) * ((sin(beta_w * x)) - (sinh(beta_w * x))));
phiu2 = matlabFunction((sin(beta_u * x))^2);
phiw2 = matlabFunction(((cos(beta_w * x)) - (cosh(beta_w) * x) + (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B))))) * ((sin(beta_w * x)) - (sinh(beta_w * x))))^2);
% Compute derivatives
phiu_x = matlabFunction(beta_u * cos(beta_u * x));
phiu_xx = matlabFunction((-beta_u^2) * sin(beta_u * x));
phiu_x2 = matlabFunction((beta_u^2 * cos(beta_u * x)^2));
phiw_x = matlabFunction((-beta_w * sin(beta_w * x)) + ((beta_w * cos(beta_w * x) - beta_w * cosh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - (beta_w * sinh(beta_w * x)));
phiw_xx = matlabFunction(((-beta_w ^2) * cos(beta_w * x)) + (((-beta_w ^2) * sin(beta_w * x) - (beta_w ^2) * sinh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - ((beta_w ^2) * cosh(beta_w * x)));
phiw_x2 = matlabFunction(((-beta_w * sin(beta_w * x)) + ((beta_w * cos(beta_w * x) - beta_w * cosh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - (beta_w * sinh(beta_w * x)))^2);
phiw_x4 = matlabFunction(((-beta_w * sin(beta_w * x)) + ((beta_w * cos(beta_w * x) - beta_w * cosh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - (beta_w * sinh(beta_w * x)))^4);
phiw_xx2 = matlabFunction((((-beta_w ^2) * cos(beta_w * x)) + (((-beta_w ^2) * sin(beta_w * x) - (beta_w ^2) * sinh(beta_w * x)) * (((sin(beta_w) * L_B) - (sinh(beta_w) * L_B) + (((beta_w * M_p)/(I0_B)) * ((cos(beta_w) * L_B) - (cosh(beta_w) * L_B)))) / ((cos(beta_w) * L_B) + (cosh(beta_w) * L_B) - (((beta_w * M_p)/(I0_B)) * ((sin(beta_w) * L_B) - (sinh(beta_w) * L_B)))))) - ((beta_w ^2) * cosh(beta_w * x)))^2);
% Define the integrals
% For K_uu
integral_1 = k0_B * integral(phiu_x2, 0, L_B);
integral_2 = k0_M * integral(phiu_x2, 0, L_M);
% For K_uw
product_phiu_x_phiw_xx = @(x) phiu_x(x) .* phiw_xx(x);
integral_3 = k1_B * integral(product_phiu_x_phiw_xx,0,L_B);
integral_4 = k1_M * integral(product_phiu_x_phiw_xx,0,L_M);
% For K_ww
integral_5 = k2_B * integral(phiw_xx2, 0, L_B);
integral_6 = k2_M * integral(phiw_xx2, 0, L_M);
% For K_uww
product_phiu_x_phiw_x2 = @(x) phiu_x(x) .* phiw_x2(x);
integral_7 = k0_B * integral(product_phiu_x_phiw_x2,0,L_B);
integral_8 = k0_M * integral(product_phiu_x_phiw_x2,0,L_M);
% For K_www
product_phiu_x2_phiw_xx = @(x) phiu_x2(x) .* phiw_xx(x);
integral_9 = k1_B * integral(product_phiu_x2_phiw_xx,0,L_B);
integral_10 = k1_M * integral(product_phiu_x2_phiw_xx,0,L_M);
% For K_wwww
integral_11 = k0_B * integral(phiw_x4, 0, L_B);
integral_12 = k0_M * integral(phiw_x4, 0, L_M);
% For M_uu
integral_13 = I0_B * integral(phiu2, 0, L_B);
integral_14 = I0_M * integral(phiu2, 0, L_M);
% For M_uw
product_phiu_phiw_x = @(x) phiu(x) .* phiw_x(x);
integral_15 = I1_B * integral(product_phiu_phiw_x,0,L_B);
integral_16 = I1_M * integral(product_phiu_phiw_x,0,L_M);
% For M_ww
integral_17 = I0_B * integral(phiw2, 0, L_B);
integral_18 = I0_M * integral(phiw2, 0, L_B);
% For Kbar_u
integral_19 = integral(phiu_x, 0, L_M);
% For Kbar_w
integral_20 = integral(phiw_xx, 0, L_M);
% For Kbar_ww
integral_21 = integral(phiw_x, 0, L_M);
% Stiffness and Mass Matrix Calculation
K_uu = integral_1 + integral_2;
K_uw = -2 .* (integral_3 + integral_4);
K_ww = integral_5 + integral_6;
K_uww = integral_7 + integral_8;
K_www = -integral_9 - integral_10;
K_wwww = 1/4 .* (integral_11 + integral_12);
M_uu = integral_13 + integral_14 + (M_p * L_B^2 * sin(beta_u)^2);
M_uw = -2 .* (integral_15 + integral_16);
M_ww = integral_17 + integral_18 + M_p * phiw2(L_B) + I_p * phiw_x2(L_B); %changed
Kbar_u =k0_M * (epsilon_0 - ((1 - zeta_f12) * epsilon_r_max)) * (integral_19);
Kbar_w =-k1_M * (epsilon_0 - ((1 - zeta_f12) * epsilon_r_max)) * (integral_20);
Kbar_ww =k0_M * (epsilon_0 - ((1 - zeta_f12) * epsilon_r_max)) * (integral_21);
Fbar = @(t)F(t) * phiw(L_B); %changed
% Time span for simulation
tspan = linspace(0, 0.01, 100); % Adjust as needed %changed
% Initial conditions [q_u(0), dq_u/dt(0), q_w(0), dq_w/dt(0)]
y0 = [0; 0; 0; 0]; % Replace with actual values
opts = odeset('RelTol',1e-6,'AbsTol',1e-8);
% Solve ODEs using ode45
[t, Y] = ode45(@(t, y) system_eqs(t, y, M_uu, M_ww, M_uw, ...
K_uu, K_ww, K_uw, K_uww, K_wwww, ...
Kbar_u, Kbar_w, Kbar_ww, Fbar), tspan, y0);
% Extract solutions
q_u = Y(:, 1);
q_w = Y(:, 3);
% Plot results
figure;
plot(t, q_u, 'b-', 'LineWidth', 2); hold on;
plot(t, q_w, 'r--', 'LineWidth', 2);
xlabel('Time (t)');
ylabel('q_u, q_w');
legend('q_u(t)', 'q_w(t)');
title('Solution of q_u and q_w over time');
grid on;
% Compute Displacements
[X, T] = meshgrid(x, t);
U = sin(beta_u * X) .* interp1(t, q_u, T, 'linear', 'extrap');
W = (cos(beta_w * X) - cosh(beta_w * X)) .* interp1(t, q_w, T, 'linear', 'extrap');
% Plot Results
figure;
plot(t, q_u, 'b-', 'LineWidth', 2); hold on;
plot(t, q_w, 'r--', 'LineWidth', 2);
xlabel('Time (t)'); ylabel('q_u, q_w');
legend('q_u(t)', 'q_w(t)');
title('Solution of q_u and q_w over time'); grid on;
figure;
surf(T, X, U, 'EdgeColor', 'none');
Error using surf (line 71)
Z must be a matrix, not a scalar or vector.
xlabel('Time (t)'); ylabel('Position (x)'); zlabel('u(x,t)');
title('Longitudinal Displacement u(x,t)');
colorbar; view(3);
figure;
surf(T, X, W, 'EdgeColor', 'none');
xlabel('Time (t)'); ylabel('Position (x)'); zlabel('w(x,t)');
title('Transverse Displacement w(x,t)');
colorbar; view(3);
% Compute voltage using numerical differentiation
V = -gradient((mio_0 * M_sat * (1 - zeta) * cos(teta_3_0) - zeta * sin(teta_4_0)), tspan) * N_coil * A_coil;
% Plot voltage over time
figure;
plot(t, V, 'r', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Voltage (V)');
title('Induced Voltage Over Time');
grid on;
P = abs((V^2) / R);
subplot(4, 1, 3);
plot(t, P);
title('Power Output');
xlabel('Time (s)'); ylabel('Power (W)');
% Function defining the system of equations
function dydt = system_eqs(t, y, M_uu, M_ww, M_uw, ...
K_uu, K_ww, K_uw, K_uww, K_wwww, ...
Kbar_u, Kbar_w, Kbar_ww, Fbar) %changed
% Extract variables
q_u = y(1); dq_u = y(2);
q_w = y(3); dq_w = y(4);
% Define mass matrix
M = [M_uu, 0.5*M_uw;
0.5*M_uw, M_ww];
% Define stiffness and nonlinear terms
K1 = [K_uu*q_u + 0.5*K_uw*q_w + K_uww*q_w^2 - Kbar_u;
K_ww*q_w + 0.5*K_uw*q_u + K_uww*q_u*q_w + 2*K_wwww*q_w^3 - (Kbar_w + Kbar_ww*q_w + Fbar(t))]; %changed
% Solve for accelerations
accel = M \ (-K1);
% Define dydt
dydt = [dq_u; accel(1); dq_w; accel(2)];
end

More Answers (0)

Categories

Find more on Correlation and Convolution in Help Center and File Exchange

Tags

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!