I'm building a parametric model of a mass damper system. The majority of this work is simple enough in simulink however I have a complex damper which requires six nonlinear equations to be solved.
I have managed to use 'fsolve' and 'vpasolve' to do this however, when I try and apply bounds to the system to recreate the limits of the physical system the vpa system breaks down. This is problem number one.
One of the things that's really confusing me is that the my equation states that the sum of the flows is equal to the velocity multiplied by the area of by supporting rod. However when I sum these flows together after the solver has completed the value is a factor of 10 out from the value given in the equation.
I've moved to vpasolve because I could constrain/apply bound values in fsolve (not that i'm having much luck at doing this with vpasolve).
Code attached below, any help, mild abuse for lack of programming prowess or general feedback would be greatly appreciated. Cheers Jon
A_rod = 0.010;
A_b = 0.020;
A_0 = 0.022;
A_v = 0.020;
a = 0.003;
b = 0.005;
Cd = 0.60;
Cf = 0.30;
dens = 1000;
D_v = 0.10;
D_p = 0.110;
L = 0.010;
k = 5000;
Mu = 0.50;
F_sp = 5;
x_dot = 5;
Pc = 300000;
syms Q_b Q_v Q_lp y Pr Pv
[Q_b_sol, Q_v_sol, Q_lp_sol, y_sol, Pr_sol, Pv_sol] = vpasolve(...
A_rod*x_dot == Q_b - Q_v - Q_lp ...
Q_b == (A_b * Cd * sqrt( (2*(Pc-Pr))/dens) ) , ...
Q_v == (A_0 * Cd * sqrt( (2*(Pc-Pv))/dens) ) , ...
Q_v == ( (a*pi*D_v*y)*Cd )*sqrt( (2*(Pv-Pr))/dens ) , ...
k*y == ( (Pv-Pr)*A_v ) - ( dens * Cf * ((Q_v/A_0)^2) ) + F_sp , ...
Q_lp ==( ( ( ((Pc-Pr)*b^3)/(12*Mu*L) ) + ( x_dot*(b/2) ) )*pi*D_p ), ...
[Q_b Q_v Q_lp y Pr Pv])
flow = Q_b_sol+Q_v_sol+Q_lp_sol;
disp(' all done ');