Energy Harvesting of Cantilever Beam with MSMA alloy

Hi, I want to write a code for the paper I attached. This papar is about harvest energy from a cantilver beam coupled with MSMA. I wrote a code but I do not know if it is correct or not and it gives some errors. Could you please someone help me in this regard? Thank you. The code:

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
integral_3 = k1_B * integral((phiu_x * phiw_xx), 0, L_B);

Operator '*' is not supported for operands of type 'function_handle'.

integral_4 = k1_M * integral((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
integral_7 = k0_B * integral((phiu_x * phiw_x2), 0, L_B);
integral_8 = k0_M * integral((phiu_x * phiw_x2), 0, L_M);

% For K_www
integral_9 = k1_B * integral((phiw_x2 * phiw_xx), 0, L_B);
integral_10 = k1_M * integral((phiw_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
integral_15 = I1_B * integral((phiu * phiw_x), 0, L_B);
integral_16 = I1_M * integral((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, x, 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))); 
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

[Hamed]

Also, if someone can write the correct code based on the paper I attached, please tell me and we can talk about it. This is important. Thanks

Answer 1

[Abhishek]

Hi,

It looks like you're running into a common issue when working with function handles in MATLAB. The error you're encountering stems from attempting to use the ‘*’ operator with function handles, which isn't supported by MATLAB. Instead, you will want to first create a new function handle by taking the element-wise product of the two function handle.

>> FE = @(x) phiu_x(x) .* phiw_xx(x);

The 'FE' will be the new function handle. Pass 'FE' to the 'integral()' function, where it will be evaluated from x = 0 to x = L_B;

>> integral_3 = k1_B * integral(FE,0,L_B);

This is a good approach since the above is using ‘.*’, which is the element-wise multiplication operator. This operator is necessary when you want to multiply arrays element-by-element, and it's especially crucial when you're working with vectors or matrices.

Please refer the documentation for more information: https://www.mathworks.com/help/releases/R2022b/matlab/matlab_prog/array-vs-matrix-operations.html

To resolve this issue, please follow these steps:

1. For getting the product of two function handles, first evaluate the function handles on data.
2. Perform element-wise multiplication using the ‘.*’ operator.

I have made changes, where the old code performed multiplication of two function handles. Please replace your sections of your old code with this:

>> 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);
>> integral_5 = k2_B * integral(phiw_xx2, 0, L_B);
>> integral_6 = k2_M * integral(phiw_xx2, 0, L_M);
>> 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);
>> 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);
>> integral_11 = k0_B * integral(phiw_x4, 0, L_B);
>> integral_12 = k0_M * integral(phiw_x4, 0, L_M);
>> integral_13 = I0_B * integral(phiu2, 0, L_B);
>> integral_14 = I0_M * integral(phiu2, 0, L_M);
>> 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);
>> integral_17 = I0_B * integral(phiw2, 0, L_B);
>> integral_18 = I0_M * integral(phiw2, 0, L_B);
>> integral_19 = integral(phiu_x, 0, L_M);
>> integral_20 = integral(phiw_xx, 0, L_M); 

By making the above changes, you should be able to resolve the error and ensure your code runs correctly. If you have any further questions or need additional help, feel free to ask. Hope these helps.

[Hamed]

Thank you for your reply dear Abhishek.

I applied all your revisions and it works. Thank you so much. There is also another 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

[Abhishek]

Hi @Hamed, can you please create a new quesiton on MATLAB Answers for that particular error?

[Hamed]

I create a post with this name: Energy Harvesting of Cantilever Beam with MSMA alloy1

Thank you