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

조회 수: 10 (최근 30일)
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
Alan Stevens 2020년 12월 1일
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
Vincent Pipitone
Vincent Pipitone 2020년 12월 1일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Walter Roberson
Walter Roberson 2020년 12월 1일
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
  댓글 수: 2
Vincent Pipitone
Vincent Pipitone 2020년 12월 8일
Your answer helped me with the next version of the program. Should I make a new thread about it?

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Stress and Strain에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by