Trouble finding solution of unknown variable in difficult integration problem

7 views (last 30 days)
I am trying to solve for 'h' using the first equation above. My script produces a result, however it is unexpected. The solution of 'h' should increase with 's', yet this is not the case. I have attached the script below, any help would be much appreciated.
clc; clear all; close all;
%%%%%%%%%% Physical Parameters %%%%%%%%%%
r=0.1;
b=0.15;
W=60;
Kc=0.99*1000;
Kphi=1528.4*1000;
n=1.1;
c=1.04*1000;
phi=28;
A1=(Kc/b+Kphi);
k=0.6;
Beta=0.0872665;
kx=0.043*Beta+0.036;
%%%%%%%%%% Defining equations %%%%%%%%%%
s=1
syms h
syms Pheta
PhetaF=acos(1-h/r);
PhetaR=acos(1-k*h/r);
jx=r*(PhetaF-Pheta-(1-s)*(sin(PhetaF)-sin(Pheta)));
Sig=A1*r^n*(cos(Pheta)-cos(PhetaF));
Taux=(c+Sig*tand(phi))*(1-exp(-jx/kx));
%%%%%%%%%% Solving for h %%%%%%%%%%
fun=Taux*sin(Pheta)+Sig*cos(Pheta);
eqn=W==r*b*vpaintegral(fun,Pheta,[PhetaR,PhetaF]);
sol=vpasolve(eqn,h,.1)

Accepted Answer

Alan Stevens
Alan Stevens on 30 Jan 2021
I used a different approach - see below. I could only see an increase in h, with increasing s, as a step change around s = 4.5.
s = 1:0.1:10;
h0 = 0.1; % Use h0 = 0.15 to get a different set of results.
for i = 1:numel(s)
h(i) = fzero(@(h) Zfun(h,s(i)), h0);
end
plot(s,h,'o'),grid
xlabel('s'),ylabel('h')
legend(['initial guess h0 = ' num2str(h0)])
disp(h)
function Z = Zfun(h,s)
r=0.1;
n=1.1;
rn = r^n;
b=0.15;
W=60;
Kc=0.99*1000;
Kphi=1528.4*1000;
k = 0.6;
c=1.04*1000;
tanphi=tand(28);
Beta=0.0872665;
kx=0.043*Beta+0.036;
sigA = (Kc/b+Kphi)*rn;
thetar = acos(1-k*h/r);
thetaf = acos(1-h/r);
sig = @(theta) sigA*(cos(theta) - cos(thetaf));
jx = @(theta) r*(thetaf - theta - (1-s)*(sin(thetaf) - sin(theta)));
taux = @(theta) (c + sig(theta).*tanphi).*(1 - exp(-jx(theta)/kx));
kern = @(theta) taux(theta).*sin(theta) + sig(theta).*cos(theta);
Z = r*b*integral(kern,thetar,thetaf) - W;
end
This results in
  5 Comments
Alan Stevens
Alan Stevens on 31 Jan 2021
Ok. Assuming you want the original values of h, then you can do the following:
s = 1:0.1:10;
h0 = 0.1; % Use h0 = 0.15 to get a different set of results.
for i = 1:numel(s)
h(i) = fzero(@(h) Zfun(h,s(i)), h0);
Fx(i) = integralfunction2(h(i),s(i));
end
plot(s,h,'o'),grid
xlabel('s'),ylabel('h')
legend(['initial guess h0 = ' num2str(h0)])
%disp(h)
figure
subplot(2,1,1)
plot(h,Fx,'o'),grid
xlabel('h'),ylabel('Fx')
subplot(2,1,2)
plot(s,Fx,'o'),grid
xlabel('s'),ylabel('Fx')
function Z = Zfun(h,s)
W=60;
Z = integralfunction(h,s) - W;
end
function Irb = integralfunction(h,s)
r=0.1;
n=1.1;
rn = r^n;
b=0.15;
Kc=0.99*1000;
Kphi=1528.4*1000;
k = 0.6;
c=1.04*1000;
tanphi=tand(28);
Beta=0.0872665;
kx=0.043*Beta+0.036;
sigA = (Kc/b+Kphi)*rn;
thetar = acos(1-k*h/r);
thetaf = acos(1-h/r);
sig = @(theta) sigA*(cos(theta) - cos(thetaf));
jx = @(theta) r*(thetaf - theta - (1-s)*(sin(thetaf) - sin(theta)));
taux = @(theta) (c + sig(theta).*tanphi).*(1 - exp(-jx(theta)/kx));
kern = @(theta) taux(theta).*sin(theta) + sig(theta).*cos(theta);
Irb = r*b*integral(kern,thetar,thetaf);
end
function Irb2 = integralfunction2(h,s)
r=0.1;
n=1.1;
rn = r^n;
b=0.15;
Kc=0.99*1000;
Kphi=1528.4*1000;
k = 0.6;
c=1.04*1000;
tanphi=tand(28);
Beta=0.0872665;
kx=0.043*Beta+0.036;
sigA = (Kc/b+Kphi)*rn;
thetar = acos(1-k*h/r);
thetaf = acos(1-h/r);
sig1 = @(theta) sigA*(cos(theta) - cos(thetaf));
jx = @(theta) r*(thetaf - theta - (1-s)*(sin(thetaf) - sin(theta)));
taux1 = @(theta) (c + sig1(theta).*tanphi).*(1 - exp(-jx(theta)/kx));
kern1 = @(theta) taux1(theta).*cos(theta)-sig1(theta).*cos(theta);
Irb2 = r*b*integral(kern1,thetar,thetaf);
end

Sign in to comment.

More Answers (1)

Alex Sha
Alex Sha on 30 Jan 2021
Hi, there are two solutions as below:
1: h=0.073660616949704
2: h=0.165217787421255
  2 Comments
Alex Sha
Alex Sha on 31 Jan 2021
If give s form 1 to 4, there will be two set of solutions, however, both of them, 'h' will decrease with 's':
solution 1:
s h
1 0.073660616949704
1.5 0.0711236790825779
2 0.0690063620895268
2.5 0.0672422418190603
3 0.0657674528467493
3.5 0.0645275378178806
4 0.0634782731873953
Solution 2:
s h
1 0.165217787441994
1.5 0.156485557183266
2 0.14921049827769
2.5 0.143711115798623
3 0.139599112575786
3.5 0.136470016301514
4 0.134031387697928

Sign in to comment.

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!