How to solve a symbolic angle for an equation without complex numbers

10 views (last 30 days)
Vincent Pipitone
Vincent Pipitone on 1 Dec 2020
Commented: Walter Roberson on 8 Dec 2020
This program takes a diameter of a nitinol coil and calculates the distance the spring would be pulled apart at a force (F_A). The problem is that using the syms function with these complicated formulas gives me complex numbers when I do vpa(solve(function==value, angle_to_be_solved_for). Also, I can't use just a solve() to get a value I can put into a tand() function to get the other values afterwards. To summarize, I need a way to solve for the angle in degrees without getting a complex number. The code in question is bolded and everything else feeds in values.
clear;
clc;
Input Values
%Shear_max=560*(10^6); %yield stress from data sheet
g_max_A=172.369*(10^6); %max cyclic shear stress austentite MPa (25ksi)
g_max_M=68.9476*(10^6); %max cyclic shear stress martensite MPa (10ksi)
G_M=10.8*(10^9); %Shear Modulus of the cold Material
G_A=28.8*(10^9); %Shear Modulus of the hot Material
v=0.33; %Poisson's ratio
%Wire and Coil Diameter
D=input('The Mean diameter of the Spring in mm; ');
D=D/1000;
d_net=0;
length_f=0.121;
%Iterating final length to find a diameter for force
while length_f >= 0.01 & d_net <= 0.07857
length_f=length_f-0.001;
dl=105;
%dl=input('The Displacement of the Spring in mm; ');
dl=dl/1000;
angle_f=90;
g=zeros(40);
while (g<=0.08) & (angle_f>=0)
angle_f=angle_f-0.1; %iterated final angle
pitch_f=tand(angle_f)*pi*D; %Coil pitch
%Length of unsprung spring set to forearm
length_i=length_f-dl; %Final length = initial + displacement
n=length_f/pitch_f; %Number of loops in coil
angle_i=atand(length_i/(pi*n*D));
%Shear Strain and Detwinned Martensite Fraction
d=150*(10^-6); %The diameter of the Wire (microns converted to m)
%Calculate martensite percent from range of shear strain
g=d./D.*cosd(angle_i).^2.*(sind(angle_f-(2*angle_i))+...
sind(angle_i))./(cosd(angle_f-(2*angle_i)).^2.*...
(cosd(angle_f-(2*angle_i)).^2.+...
(sind(angle_f-(2*angle_i)).^2)./(1+v)));
end
angle_range=angle_i:0.1:angle_f;
g_range=d./D.*cosd(angle_i).^2.*(sind(angle_range-(2*angle_i))+...
sind(angle_i))./(cosd(angle_range-(2*angle_i)).^2.*...
(cosd(angle_range-(2*angle_i)).^2.+...
(sind(angle_range-(2*angle_i)).^2)./(1+v)));
%plot(angle_range,g_range)
%Defining critical shear stresses and variable shear stress
gcr_s=0.01;
gcr_f=0.12;
%Detwinned martensite fraction
ESg=0.5.*cos(pi.*(g_range-gcr_f)./(gcr_s-gcr_f))+0.5;
eigth=round(length(ESg)*0.01/g);
ESg(1:eigth)=0;
%plot(g_range,ESg)
%Final angles of austentite and martensite phase finish (spring range)
syms angle_A_f
%Displacement of the two metal phases for 100g hystersis force
F_A=G_A*d^4/(8*D^3*n)*cosd(angle_i)^3/(cosd(angle_A_f)^2*(cosd(angle_A_f)...
^2+sind(angle_A_f)^2/(1+v)))*n*pi*D*tand(angle_A_f);
F_M=G_M*d^4/(8*D^3*n)*(cosd(angle_i)^3)/(cosd(angle_f)^2*(cosd(angle_f)...
^2+sind(angle_f)^2/(1+v)))*dl-pi*d^3/(8*D)*G_M*ESg*0.06; %residual g
angle_A_f=(solve(F_A==.981, angle_A_f)); %MARKED
d_A=n*pi*D*tand(angle_A_f); %MARKED
d_net=dl-d_A;
end
disp(d_net)

Sign in to answer this question.

Accepted Answer

Alan Stevens
Alan Stevens on 1 Dec 2020
You could find the angle with fzero:
%Input Values
%Shear_max=560*(10^6); %yield stress from data sheet
g_max_A=172.369*(10^6); %max cyclic shear stress austentite MPa (25ksi)
g_max_M=68.9476*(10^6); %max cyclic shear stress martensite MPa (10ksi)
G_M=10.8*(10^9); %Shear Modulus of the cold Material
G_A=28.8*(10^9); %Shear Modulus of the hot Material
v=0.33; %Poisson's ratio
%Wire and Coil Diameter
%D=input('The Mean diameter of the Spring in mm; ');
D = 5;
D=D/1000;
d_net=0;
length_f=0.121;
%Iterating final length to find a diameter for force
while (length_f >= 0.01) & (d_net <= 0.07857)
length_f=length_f-0.001;
dl=105;
%dl=input('The Displacement of the Spring in mm; ');
dl=dl/1000;
angle_f=90;
g=zeros(40);
while (g<=0.08) & (angle_f>=0)
angle_f=angle_f-0.1; %iterated final angle
pitch_f=tand(angle_f)*pi*D; %Coil pitch
%Length of unsprung spring set to forearm
length_i=length_f-dl; %Final length = initial + displacement
n=length_f/pitch_f; %Number of loops in coil
angle_i=atand(length_i/(pi*n*D));
%Shear Strain and Detwinned Martensite Fraction
d=150*(10^-6); %The diameter of the Wire (microns converted to m)
%Calculate martensite percent from range of shear strain
g=d./D.*cosd(angle_i).^2.*(sind(angle_f-(2*angle_i))+...
sind(angle_i))./(cosd(angle_f-(2*angle_i)).^2.*...
(cosd(angle_f-(2*angle_i)).^2.+...
(sind(angle_f-(2*angle_i)).^2)./(1+v)));
end
angle_range=angle_i:0.1:angle_f;
g_range=d./D.*cosd(angle_i).^2.*(sind(angle_range-(2*angle_i))+...
sind(angle_i))./(cosd(angle_range-(2*angle_i)).^2.*...
(cosd(angle_range-(2*angle_i)).^2.+...
(sind(angle_range-(2*angle_i)).^2)./(1+v)));
%plot(angle_range,g_range)
%Defining critical shear stresses and variable shear stress
gcr_s=0.01;
gcr_f=0.12;
%Detwinned martensite fraction
ESg=0.5.*cos(pi.*(g_range-gcr_f)./(gcr_s-gcr_f))+0.5;
eigth=round(length(ESg)*0.01/g);
ESg(1:eigth)=0;
%plot(g_range,ESg)
%Final angles of austentite and martensite phase finish (spring range)
angle_A_f0 = 30;
angle_A_f = fzero(@(angle_A_f) findangle(angle_A_f,G_A,d,D,angle_i,v,n),angle_A_f0);
%Displacement of the two metal phases for 100g hystersis force
%F_A=G_A*d^4/(8*D^3*n)*cosd(angle_i)^3/(cosd(angle_A_f)^2*(cosd(angle_A_f)...
% ^2+sind(angle_A_f)^2/(1+v)))*n*pi*D*tand(angle_A_f);
%F_M=G_M*d^4/(8*D^3*n)*(cosd(angle_i)^3)/(cosd(angle_f)^2*(cosd(angle_f)...
% ^2+sind(angle_f)^2/(1+v)))*dl-pi*d^3/(8*D)*G_M*ESg*0.06; %residual g
%angle_A_f=(solve(F_A==.981, angle_A_f)); %MARKED
d_A=n*pi*D*tand(angle_A_f); %MARKED
d_net=dl-d_A;
end
disp(angle_A_f)
disp(d_net)
function Z = findangle(angle_A_f,G_A,d,D,angle_i,v,n)
Z = G_A*d^4/(8*D^3*n)*cosd(angle_i)^3/(cosd(angle_A_f)^2*(cosd(angle_A_f)...
^2+sind(angle_A_f)^2/(1+v)))*n*pi*D*tand(angle_A_f) - 0.981;
end
  1 Comment
Vincent Pipitone
Vincent Pipitone on 1 Dec 2020
It looks like the angle can be found that way, by setting the function equal to zero and solving for it. I was able to find the austentite angle by this method and it seems to make sense, thank you.

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 1 Dec 2020
Ran in:
You can use vpasolve() with range [-inf inf]
Caution: if you copy this code back, fix the entry for D back to use input()... the testing facility I used does not accept input() so I used rand()
%Shear_max=560*(10^6); %yield stress from data sheet
g_max_A=172.369*(10^6); %max cyclic shear stress austentite MPa (25ksi)
g_max_M=68.9476*(10^6); %max cyclic shear stress martensite MPa (10ksi)
G_M=10.8*(10^9); %Shear Modulus of the cold Material
G_A=28.8*(10^9); %Shear Modulus of the hot Material
v=0.33; %Poisson's ratio
%Wire and Coil Diameter
%D=input('The Mean diameter of the Spring in mm; ');
D = rand() * 10; %<==== FOR ONLINE TESTING
D=D/1000;
d_net=0;
length_f=0.121;
%Iterating final length to find a diameter for force
while length_f >= 0.01 & d_net <= 0.07857
length_f=length_f-0.001;
dl=105;
%dl=input('The Displacement of the Spring in mm; ');
dl=dl/1000;
angle_f=90;
g=zeros(40);
while (g<=0.08) & (angle_f>=0)
angle_f=angle_f-0.1; %iterated final angle
pitch_f=tand(angle_f)*pi*D; %Coil pitch
%Length of unsprung spring set to forearm
length_i=length_f-dl; %Final length = initial + displacement
n=length_f/pitch_f; %Number of loops in coil
angle_i=atand(length_i/(pi*n*D));
%Shear Strain and Detwinned Martensite Fraction
d=150*(10^-6); %The diameter of the Wire (microns converted to m)
%Calculate martensite percent from range of shear strain
g=d./D.*cosd(angle_i).^2.*(sind(angle_f-(2*angle_i))+...
sind(angle_i))./(cosd(angle_f-(2*angle_i)).^2.*...
(cosd(angle_f-(2*angle_i)).^2.+...
(sind(angle_f-(2*angle_i)).^2)./(1+v)));
end
angle_range=angle_i:0.1:angle_f;
g_range=d./D.*cosd(angle_i).^2.*(sind(angle_range-(2*angle_i))+...
sind(angle_i))./(cosd(angle_range-(2*angle_i)).^2.*...
(cosd(angle_range-(2*angle_i)).^2.+...
(sind(angle_range-(2*angle_i)).^2)./(1+v)));
%plot(angle_range,g_range)
%Defining critical shear stresses and variable shear stress
gcr_s=0.01;
gcr_f=0.12;
%Detwinned martensite fraction
ESg=0.5.*cos(pi.*(g_range-gcr_f)./(gcr_s-gcr_f))+0.5;
eigth=round(length(ESg)*0.01/g);
ESg(1:eigth)=0;
%plot(g_range,ESg)
%Final angles of austentite and martensite phase finish (spring range)
syms angle_A_f
%Displacement of the two metal phases for 100g hystersis force
F_A=G_A*d^4/(8*D^3*n)*cosd(angle_i)^3/(cosd(angle_A_f)^2*(cosd(angle_A_f)...
^2+sind(angle_A_f)^2/(1+v)))*n*pi*D*tand(angle_A_f);
F_M=G_M*d^4/(8*D^3*n)*(cosd(angle_i)^3)/(cosd(angle_f)^2*(cosd(angle_f)...
^2+sind(angle_f)^2/(1+v)))*dl-pi*d^3/(8*D)*G_M*ESg*0.06; %residual g
angle_A_f=(vpasolve(F_A==.981, angle_A_f, [-inf inf])); % <====== HERE
d_A=n*pi*D*tand(angle_A_f);
d_net=dl-d_A;
end
disp(vpa(d_net))
91.943581578108842623752027152327

Categories

Find more on Stress and Strain in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!