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

31 views (last 30 days)
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)

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 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.

Walter Roberson on 1 Dec 2020
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

Vincent Pipitone on 8 Dec 2020
Walter Roberson on 8 Dec 2020
Sure