syms beta alpha L L2 x gamma real

g = 0.85;
Kt = 2.65;
E = 19e9; % Pa
b = 1/4*0.0254;

% CHANGE THESE PARAMETERS %%%%%%%%
l = 8e-2; %rigid part, should be less than 85% of l2
l2 = 15e-2; % compliant part
h = 0.5e-2; %5mm
I = b*h^3/12;
K = 4*g*Kt*E*I/l2;% 2x stiffness of one beam to account for both sides

%find angle of flexure relative to x (vertical)
eq1 = L/sin(beta) == gamma*L2/sin(alpha+beta); %use law of sines to solve for beta
eq2 = L/sin(beta) == ((gamma*L2 - L) + x)/sin(alpha);
sol = solve([eq1, eq2], [beta, alpha]);

a = subs(abs(sol.alpha), [L, L2, gamma], [l, l2, g]);
b1 = subs((sol.beta(1)), [L, L2, gamma], [l, l2, g]); % angle for PRBM
b2 = subs((sol.beta(2)), [L, L2, gamma], [l, l2, g]);
a_func = matlabFunction(a)
b1_func = matlabFunction(b1) %convert symbolic solution to function
b2_func = matlabFunction(b2)

dx = 0.001;
xspan = 0:dx:2*l-dx;

%calculate angle of rigid component relative to x (vertical)
th =  b1 + a ; %only outputs -pi/2 to pi/2 and we need 0 to pi so need to split
th_func = matlabFunction(th);
hold off;
plot(xspan, th_func(xspan))
hold on
% plot(xspan, b1_func(xspan));legend('th', 'beta')
[~,bound] = max(th_func(xspan)); %% find spliting point
xspan1 = 0:dx:xspan(bound);
xspan2 = xspan(bound)+dx:dx:2*l-dx;

Theta = [max(th_func(xspan1)), max(pi-th_func(xspan2))]; % make a lookup table for theta values then interpolate
Theta_func = @(xq) interp1(xspan,Theta, xq) ; %interpolate

U = @(x) K.*b1_func(x).^2/2 %energy
F = @(x) K.*b1_func(x) % force exerted by flexure onto rigid component at linkage
Fx = @(x) F(x).*cos(Theta_func(x)) % force exerted by rigid component
Kx = @(x) -diff(Fx(xspan)) %numerically differentiate force to get stiffness

hold off;
figure;
plot(xspan, U(xspan))
xlabel('x')
ylabel('Energy (J)')

figure;
plot(xspan, Fx(xspan))
xlabel('x')
ylabel('Force (N)')

figure;
plot(xspan, abs(F(xspan).*l2.*(h/2)./I./2).*10^(-9)) %My/I = FL(t/2)/I then divide by 2 since the force is for both beams
xlabel('x')
ylabel('bending stress (GPa)')

figure;
plot(xspan(1:end-1), Kx(xspan))
xlabel('x')
ylabel('Stiffness (N/m)')

Published with MATLAB® R2023b