need help to calculate a dependent variable outside of function using ode45

3 views (last 30 days)
md mayeen uddin
md mayeen uddin on 3 Sep 2021
Commented: Walter Roberson on 4 Sep 2021
Below is my function and the code to call that function using ode45 for a first order non-homogenius equation. THe code works fine , and the answer is also correct. But i just want to calculate kf outside of the program, to decrease the runtime
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = FunctionOfNacheilzone(N,q) % Function to calculate stress in lagging zone
Radius_R = 21; % Radius of idle Roll
n = 1;
%___________________________________________________________________________
Wire_dH = 8; % Diameter of hot rolled copper
Wire_d0 = 2; % Reduced diameter for roll drwaing operation
meu = 0.1; % Friction acting in between wire and rolls
C = 522; % Coefficient of Ludwik equation
M = 0.368; % Coefficiant of Ludwik equation
%___________________________________________________________________________
h_1 = [1.7900 ; 1.7000; 1.6000; 1.5000; 1.4000]; % Final height after reduction
h_0 = [2; 2; 2; 2; 2 ]; % Initi height of copper wire
b_1 = [2.0400; 2.0900; 2.1400; 2.2100; 2.2700]; % Final width after roll drawing
b_0 = [2; 2; 2; 2; 2 ]; % Initial width of copper wire
Del_h = abs(h_1(n,1) - h_0(n,1)); % Relative height reduction
l_d = sqrt(Radius_R.*Del_h); % Bite angle
b = ((b_1(n,1)-b_0(n,1))./l_d)*(l_d-N) +b_0(n,1); % Width as a function of position shift in XY plane
h = h_1(n,1)+(N^2/Radius_R); % Height as a function of position shift in XY plane
pi_b = log(b./b_0(n,1)); % Degree of deformation along with width
pi_h = log (h./h_0(n,1)); % Degree of deformation along with height
pi_l = (-pi_b - pi_h); % Degree of deformation along with length
pi_eq = sqrt((2/3).*(((pi_l.^2)+(pi_b.^2)+(pi_h.^2)))); % Equivalent degree of deformation
pi_0 = abs(log(Wire_d0^2/Wire_dH^2)); % Initial degree of deformation
pi_1 = pi_eq + pi_0; % Increase of degree of deformation with increase of span
kf = C.*(pi_1).^(M); % Kf values at every point by Ludwik equation
%_________________________________________________________________________
% the ordinary differential equation of first order to determine stress in
% the lagging zone
dqdN = -(((2./(h_1+N^2/Radius_R)).*((tan(N/Radius_R)-tan((N/Radius_R)-atan(meu))))).*q)...
-((2./(h_1+N^2/Radius_R)).*kf.*tan((N/Radius_R)-atan(meu)));
f = dqdN;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 1;
Radius_R = 21; % Radius of idle roll
h_1 = h1(1:end,:); % Height after reduction
h_0 = h0(1:end,:); % Initial height of copper wire
Delta_h = abs(h_1(n,1) - h_0(n,1)); % Reltive height reduction
In_alf = 0; % Final condition of bite angle for lagging zone
Fi_alf = sqrt(Radius_R.*Delta_h); % Initial condition of bite angle
Nspan = [Fi_alf:-0.0001:In_alf]; % Span of bite angle
%____________________________________________________
q0 = [0;0;0;0;0]; % Intital condition of stress in the lagging zone for ODE
[Nout, q] = ode45(@FunctionOfNacheilzone,Nspan,q0); % Calling the ode function using MATLAB built-in function

Answers (1)

Walter Roberson
Walter Roberson on 3 Sep 2021
kf = 522*(((2*(log(51/50 - N/105) + log(N^2/42 + 179/200))^2)/3 + (2*log(N^2/42 + 179/200)^2)/3 + (2*log(51/50 - N/105)^2)/3)^(1/2) + 6243314768165359/2251799813685248)^(46/125)
That needs N to calculate, so you cannot pre-calculate it.
A few parts of it can be pre-calculated, but not so much.
  2 Comments
Walter Roberson
Walter Roberson on 4 Sep 2021
No, that does not help.
syms N
syms q [5 1]
Radius_R = 21; % Radius of idle Roll
n = 1;
%___________________________________________________________________________
Wire_dH = 8; % Diameter of hot rolled copper
Wire_d0 = 2; % Reduced diameter for roll drwaing operation
meu = 0.1; % Friction acting in between wire and rolls
C = 522; % Coefficient of Ludwik equation
M = 0.368; % Coefficiant of Ludwik equation
%___________________________________________________________________________
h_1 = [1.7900 ; 1.7000; 1.6000; 1.5000; 1.4000]; % Final height after reduction
h_0 = [2; 2; 2; 2; 2 ]; % Initi height of copper wire
b_1 = [2.0400; 2.0900; 2.1400; 2.2100; 2.2700]; % Final width after roll drawing
b_0 = [2; 2; 2; 2; 2 ]; % Initial width of copper wire
Del_h = abs(h_1(n,1) - h_0(n,1)); % Relative height reduction
l_d = sqrt(Radius_R.*Del_h); % Bite angle
b = ((b_1(n,1)-b_0(n,1))./l_d)*(l_d-N) +b_0(n,1); % Width as a function of position shift in XY plane
h = h_1(n,1)+(N^2/Radius_R); % Height as a function of position shift in XY plane
pi_b = log(b./b_0(n,1)); % Degree of deformation along with width
pi_h = log (h./h_0(n,1)); % Degree of deformation along with height
pi_l = (-pi_b - pi_h); % Degree of deformation along with length
pi_eq = sqrt((2/3).*(((pi_l.^2)+(pi_b.^2)+(pi_h.^2)))); % Equivalent degree of deformation
pi_0 = abs(log(Wire_d0^2/Wire_dH^2)); % Initial degree of deformation
pi_1 = pi_eq + pi_0; % Increase of degree of deformation with increase of span
kf = C.*(pi_1).^(M); % Kf values at every point by Ludwik equation
simplify(kf)
ans = 
%_________________________________________________________________________
% the ordinary differential equation of first order to determine stress in
% the lagging zone
dqdN = -(((2./(h_1+N^2/Radius_R)).*((tan(N/Radius_R)-tan((N/Radius_R)-atan(meu))))).*q)...
-((2./(h_1+N^2/Radius_R)).*kf.*tan((N/Radius_R)-atan(meu)));
f = dqdN;
f
f = 
Your kf depends upon N, which is the current "time" (or at least the independent variable.) Pre-computing would only work if kf was independent of the inputs N and q, or if it depended on N and q in "simple" ways so that you could break out part of the computation . But kf does not depend upon N and q in simple ways: that entire sqrt() expression needs to be evaluated according to the current value of the independent variable.
Please also notice the in the above output of f, which has . If 51/50 - N/105 is negative then the log() will be complex. Can we prove that N/105 will be less than 51/50 ?
If we make guesses that h1 is the same as h_1 and h0 is the same as h_0 then
Radius_R = 21; % Radius of idle roll
Delta_h = abs(h_1(n,1) - h_0(n,1)); % Reltive height reduction
In_alf = 0; % Final condition of bite angle for lagging zone
Fi_alf = sqrt(Radius_R.*Delta_h); % Initial condition of bite angle
Nspan = [Fi_alf:-0.0001:In_alf]; % Span of bite angle
max(Nspan)/105
ans = 0.0200
In this case the answer is Yes, N/105 is comfortably less than 51/50 ... provided that our guess was right about h1 vs h_1 and h0 vs h_0

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by