restart; # Bernoulli equation eq := (1/2)*rho[p]*(v-u)^2+Yp=(1/2)*rho[t]*u^2+Rt: pen := solve(eq, u); # We take the second solution # because of u < v and simplify it u_v[1] := simplify(pen[2]): simplify(series(u_v[1],rho[p]=rho[t]),symbolic): u_v[2] := subs({rho[t]=rho,rho[p]=rho},op(1,%)); # Definition of initial value problem deq := diff(l(v), v)=l(v)*rho[p]*(v-u(v))/Yp: incond := l(v0) = L: len[1] := simplify (rhs (dsolve({deq, incond}, l(v)))); subs(u=unapply(u_v[2],v),{rho[p]=rho,rho[t]=rho},value(len[1])): len[2] := simplify (expand (%)); u_v[3] := simplify (subs(Rt=Yp, u_v[1]), symbolic): subs(u=unapply(u_v[3],v),len[1]): l := simplify(%, symbolic): pe := simplify(int(l*u_v[3], v=0..v0)); simplify(limit(pe,rho[t]=rho[p]),symbolic); vm := sqrt(2*(m-1)*Yp/rho): simplify (subs(Rt=m*Yp,len[2]*u_v[2])): pe := (rho/Yp)*int(%, v=vm..v0): pe1:=factor(simplify(expand(limit(pe,m=1)))); pe3 := simplify(expand(subs(m=3, pe))); restart: # Bernoulli equation for same densities rho=rho[t]=rho[p] eqe := (1/2)*rho*(v(t)-u)^2+Yp=(1/2)*rho*u^2+Rt: u := solve(eqe,u): deqe1 := rho*l(t)*diff(v(t), t)=-Yp: deqe2 := diff(l(t), t)=-(v(t)-u): incond := l(0)=L, v(0)=v0: sol := dsolve({deqe1, deqe2, incond}, {l(t), v(t)}, series): # Polynomial approximation of v(t) and l(t) v := subs(sol,v(t)); l := subs(sol,l(t)): # Polynomial approximation of u(t) uappr := convert(series(v/2+(Yp-Rt)/(rho*v),t=0,6),polynom): # First few coefficients of uappr series(uappr,t=0,2); rho := 7810: # Target and penetrator density L := 0.25: # Length of rod # Target resistance Rt = 900*10^6 # Strength of projectile Yp = 900*10^6 u := subs(Rt=900*10^6,Yp=900*10^6,uappr): plot3d(u,t=0..0.0004,v0=1000..2500, axes=boxed,orientation=[-30,60]); restart; assume(mp>0); de := -k(mu)*DP(P)*diff(DP(P), P)/(a + c(mu)*DP(P)^2) =1/(mp - mu*P); Dsu:=[diff(DP(P),P)=diff(P(t),t,t)/diff(P(t),t),DP(P)=diff(P(t),t),P=P(t)]; map(u->simplify(u*denom(lhs(de))*denom(rhs(de))),de): deq:=subs(Dsu,%); ic := DP(0) = vp: pen := {dsolve({de,ic},DP(P))}; sol:=map(u->`if`(op(1,rhs(u))=-1,NULL,rhs(u)),pen)[]: su_1:=c(mu)/mu/k(mu)=epsilon(mu)/2; bs_1:=epsilon(mu)=solve(su_1,epsilon(mu)): t_1 := simplify(subs(su_1,mp=Mp,sol)): su_2:=(Mp-mu*P)^(epsilon(mu))=A^epsilon(mu)*Mp^(epsilon(mu)); bs_2:=A=expand(solve(su_2,A)): t_2:=expand((subs(su_2,numer(t_1)^2/a/c(mu)))): dP:=subs(bs_2,sqrt(convert(t_2,horner,A^epsilon(mu))*a/c(mu))); Pf:=simplify(solve(dP=0,P)); # Projectile and target data and the final penetration depth rho[t] := 7810: Mp := 1.491: Yt := 9.70*10^8: rho[p] := rho[t]: dp := 0.0242: Yp := 15.81*10^8: Ap := Pi*dp^2/4: v0 := 1457.9: pf := 0.245: # 1st step - Determination of mup and graph penetration depth # vs erosion rate mu, the definition of constants ap, a, b # and initial penetration velocity ap := 1/(2*rho[p]*Ap): a := 3*Yt*Ap: b := rho[t]*Ap/2: vp := v0/2+(Yt-Yp)/(rho[t]*v0): # The definition of functions k, c, eps and Pf of mu k := unapply(1+ap*mu, mu): c := unapply(b-ap*mu^2,mu): epsilon:=unapply(rhs(bs_1),mu): # - Graph of influence of mu on the penetration depth plot(Pf,mu=0..6); # Determination of limit Pfzero := limit(Pf,mu=0): # Solution of the equation Pf(mu)=pf mup := fsolve(Pf=pf,mu): # 2nd step - Graph of the penetration velocity # along the penetration path plot(subs(mu=mup,dP),P=0..pf); # 3rd step - Graph of the time dependence of the penetration tf := 0.0006: # The evaluation of the time dependence # of the penetration and its velocity deq:=subs(mp=Mp,mu=mup,deq): ini := P(0) = 0, D(P)(0) = vp: F := dsolve({deq, ini}, P(t), numeric): plots[odeplot](F, [t, P(t)], 0..0.0006); # 4th step - The evaluation of the influence # of target strength and impact velocity vp:=vi + (y - Yp)/(rho[t]*vi): plot3d(subs(mu=mup,Pf),vi=1400..2200,y=900*10^6..2000*10^6,orientation=[-125,60], axes=boxed,color=white);