Question 3.9.4: Consider the problem of Example 3.6.1 and compute the respon......
Consider the problem of Example 3.6.1 and compute the response numerically. Compare the result to the analytical solution computed in Example 3.6.1 by plotting both the analytical solution and the numerical solution on the same graph.
The equation of motion to be solved in this example is
10 \ddot{x}(t)+10 x(t)=\frac{t}{t_1}-\left(\frac{t-t_1}{t_1}\right) \Phi\left(t-t_1\right) (3.110)
where Φ is the Heaviside step function used to “turn on” the second term in the forcing function at t > t_1. The parameter t1 is used to control how steep the disturbance is. See Figure 3.16 for a plot of the driving force. The analytical solution is given in equation (3.78) to be
x(t)=\frac{F_0}{k}\left(\frac{t}{t_1}-\frac{\sin \omega_n t}{\omega_n t_1}\right)-\frac{F_0}{k}\left(\frac{t-t_1}{t_1}-\frac{\sin \omega_n\left(t-t_1\right)}{\omega_n t_1}\right) \Phi\left(t-t_1\right) (3.111)
Again, the equation of motion can be solved numerically in the form given by equation (3.110) directly in Mathematica or by putting equation (3.110) into state-space form and solving Runge–Kutta in Mathcad or Matlab. The plots of both the numerical solution and analytical solution are given in Figure 3.25.
The Mathcad code for this solution follows:
| x0 := 0 v0 := 0 m := 10 t1 := 1 k := 10 \omega n:=\sqrt{\frac{k}{m}} ωn = 1 F0 := 1 tp :=\frac{1}{\omega n } \cdot \operatorname{atan}\left(\frac{1-\cos (\omega n \cdot t 1)}{\sin (\omega n \cdot t 1)}\right) \quad \omega n \cdot tp =0.5 \quad tp =0.5 x m(t):=\frac{F 0}{k} \cdot\left[1+\frac{1}{\omega n \cdot t 1} \cdot \sqrt{2 \cdot(1-\cos (\omega n \cdot t 1))}\right] x a(t):=\frac{F 0}{k} \cdot\left(\frac{t}{t 1}-\frac{\sin (\omega n \cdot t)}{\omega n \cdot t 1}\right)-\frac{F 0}{k}\left[\frac{t-t 1}{t 1}-\frac{\sin [\omega n \cdot(t-t 1)]}{\omega n \cdot t 1}\right] \Phi(t-t 1) f(t):=\frac{t}{t 1} \cdot \frac{F 0}{m}-\frac{t-t 1}{t 1} \cdot \frac{F 0}{m} \cdot \Phi\left(t-t_1\right) \quad x:=\left[\begin{array}{l} 0 \\ 0 \end{array}\right] D(t, x):=\left[\begin{array}{c} X_1 \\ -\left(\omega n^2 \cdot X_0\right)+f(t) \end{array}\right] \quad x m(0)=0.196 Z := rkfixed (X, 0, 15, 2000, D) t := Z^{<0>} x := \overrightarrow{xa(t)} xn := Z^{<1>} F := \overrightarrow{f(t)} Xmax := \overrightarrow{xm(t)} |
Script File
The Matlab code for solving for the response is
clear all
%% analytical solution
t=0:0.01:15;
m=10; k=10; Fo=1; t1=1;
wn=sqrt(k/m);
Heaviside=stepfun(t, t1);% define Heaviside Step function for 0<t<15
for i=1:max(length(t)),
xt(i)=Fo/k*(t(i)/t1 - sin(wn*t(i))/wn/tl) - Fo/k*((t(i)-t1)/t1 –
sin(wn*(t(i)-t1))/wn*t1)*Heaviside(i);
end
plot(t,xt,'- -'); hold on
%% Numerical Solution
xo=[0; 0];
ts=[0 15];
[t, x]=ode45('f', ts, xo);
plot(t, x(:,1)); hold off
%--------------------------------------------
function v=f(t, x)
m=10; k=10; wn=sqrt(k/m); Fo=1; t1=1;
v=[x(2); x(1).*-wn^2 + t/t1*Fo/m-(t-t1)/tl*Fo/m*stepfun(t, t1)];
The Mathematica code for solving for the response is
| In[1]:= 6 6PlotLegends' In[2]:= x0 = 0; v0 = 0; m = 10; k = 10; \omega n =\sqrt{\frac{ k }{ m }} ; t1 = 1; F0 = 1; tp = 0.5; \operatorname{In}[10]:=\operatorname{xanal}[ t ]=\frac{ F 0}{ k } *\left(\frac{ t }{ t 1}-\frac{\operatorname{Sin}\left[\omega n ^* t \right]}{\omega n ^* t 1}\right)-\frac{ F 0}{ k } *\left(\frac{ t - t 1 }{ t1 }-\frac{\operatorname{Sin}\left[\omega n ^*( t - t 1)\right]}{\omega n ^* t 1}\right) * UnitStep[t - t1]; xnum = NDSolve[{10 * x''[t] + 10 * x[t] == \frac{ t }{ t 1}-\left(\frac{ t - t 1}{ t 1}\right) * UnitStep[t - t1], x[0] == x0, x'[0] == v0}, x[t], {t, 0, 15}]; Plot[{Evaluate[x[t] /. xnum], xanal[t]}, {t, 0, 15}, PlotStyle → {RGBColor[1, 0, 0], RBColor[0, 1, 0]}, PlotLegend → {"Numerical", "Analytical"}, LegendPosition → {1, 0}, LegendSize S {1, 0.5}]; |