############################################################### # # This file contains the example and all necessary routines # and informations for the paper # # `Symbolic Computation of Longitudinal Wave Propagation' # # by Bin Hu and Peter Eberhard # # submitted on November 2, 1999 to the Journal # `Computer Methods in Applied Mechanics and Engineering' # #-------------------------------------------------------------- # # Computing and visualizing wave propagations for a thin rod # at one end fixed and one end struck by a moving rigid body # #-------------------------------------------------------------- # displacement #-------------------------------------------------------------- L2 := 2*L : u := (x,t)-> f(c*t-x)-f(c*t+x-L2): #-------------------------------------------------------------- # velocity #-------------------------------------------------------------- v := (x,t) -> c*df(c*t-x) -c*df(c*t+x-L2) : #-------------------------------------------------------------- # strain #-------------------------------------------------------------- epsilon := (x,t)-> -df(c*t-x)-df(c*t+x-L2) : #-------------------------------------------------------------- # stress #-------------------------------------------------------------- sigma := (x,t)->Emod*epsilon(x,t) : #-------------------------------------------------------------- # The functions df(y) and f(y) are two wave form functions # which are piecewise defined in two intervals: # during the impact 0 <= y <= ctc, # and after the impact y > ctc #-------------------------------------------------------------- ctc := c*tc : df := proc(y) if y <= ctc then df_contact(y) else df_free(y) fi end: f := proc(y) if y <= ctc then f_contact(y) else f_free(y) fi end: #-------------------------------------------------------------- # During the impact, the function f'(y) is described # by df_contact(y) #-------------------------------------------------------------- df_contact := proc(y) local yi, n, z: yi := y/L2: n := floor(yi): z := yi - n: v2/c*s(n,z,beta); end: #-------------------------------------------------------------- # and the function f(y) is described by f_contact(y) #-------------------------------------------------------------- u0 := L2/beta*v2/c: f_contact := proc(y) local yi, n, z: yi := y/L2: n := floor(yi): z := yi - n: u0*ss(n,z,beta) : end: #-------------------------------------------------------------- # The function s(n,z,beta) is used to constitute # function df_contact(y) #-------------------------------------------------------------- s := (n,z,beta) -> exp(-beta*z)*Ps(n,z,beta) : #-------------------------------------------------------------- # Polynomial Ps(n,z,beta) is used to constitute the s(n,z,beta) #-------------------------------------------------------------- Ps := proc(n,z,beta) local P,xx; option remember; if n = -1 then 0 elif n = 0 then 1 else xx:='xx': simplify(Ps(n-1,z,beta) + b*(Ps(n-1,1,beta) - Ps(n-2,1,beta)) -2*beta*int(Ps(n-1,xx,beta),xx=0..z)) ; fi end: #-------------------------------------------------------------- # The function ss(n,z,beta) is used to constitute the function f2(y) #-------------------------------------------------------------- ss := proc(n,z,beta) local xx; option remember; if n < 0 then 0 else xx:='xx': ss(n-1,1,beta) + simplify(beta*int(s(n,xx,beta),xx=0..z)) : simplify(subs(exp(-beta)=b,%)) fi end : #-------------------------------------------------------------- # The function ss(n,z,beta) can be written as # exp(-beta*z)*Pss(n,z,beta)+Css(n) #-------------------------------------------------------------- Pss := proc(n,z,beta) local d : diff(subs(exp(-beta*z)=d, ss(n,z,beta)),d) end: Css := proc(n) subs(exp(-beta*z)=0, ss(n,z,beta)) end: #-------------------------------------------------------------- # searching the instant tc and computing the maximum stress #-------------------------------------------------------------- contactduration := proc(beta) local searching, i, p0, p1, P, z, nc, zc: global e, tc, stress_max : searching := true: stress_max := 0 : for i from 1 while (searching) do p0 := evalf(Ps(i,0,beta)) : stress_max := max(stress_max, p0) : p0 := evalf(p0 + Ps(i-1,0,beta)); p1 := evalf(Ps(i,1,beta)+Ps(i-1,1,beta)); if p1 < 0 then z := 'z': P:= evalf(Ps(i,z,beta)+Ps(i-1,z,beta)): nc := i; zc := fsolve(P=0,z,0..1,maxsols=1) : searching := false : fi : od: #-------------------------------------------------------------- # the dimensionless maximum stress #-------------------------------------------------------------- stress_max := 2*stress_max : #-------------------------------------------------------------- # the coefficient of restitution #-------------------------------------------------------------- e := evalf(exp(-beta*zc)*(-Ps(nc,zc,beta)+Ps(nc-1,zc,beta))) : #-------------------------------------------------------------- # the dimensionless duration of contact #-------------------------------------------------------------- tc := (nc+zc): end : #-------------------------------------------------------------- # After the impact, the function f'(y) and f(y) are described # by df_free(y) and f_free(y), respectively. #-------------------------------------------------------------- df_free := proc(y) local iy, yi: yi := y - ctc : iy := ceil(yi/L2): yi := yi-iy*L2+ctc: (-1)^(iy)*df_contact(yi) end: f_free := proc(y) local iy, yi: yi := y - ctc : iy := ceil(yi/L2): yi := yi-iy*L2+ctc: (-1)^(iy)*f_contact(yi)+(iy mod 2)*constant1 end : #-------------------------------------------------------------- # Analytical expressions for Ps(n,z) and Pss(n,z) # # To avoid very long computation times, they are computed # here only for n < 10 and alpha >= 0.05 #-------------------------------------------------------------- # #output of function Ps(i,z), Pss(i,z) and coefficients Css(i) # for i from 0 to 9 do printf(`%3s%d%4s\n`, `Ps(`,i, `,z) = `) ; Ps(i,z,beta); printf(`%4s%d%4s\n`, `Pss(`,i, `,z) = `) ; Pss(i,z,beta); printf(`%4s%d%4s\n`, `Css(`,i, `) = `) ; Css(i); od; # #output of recursive expressions for Ps(i,z) and Pss(i,z) ) #with b = exp(-b # for i from 5 to 9 do printf(`%3s%d%8s%d%8s\n`, `Ps(`,i, `,z)-bPs(`,i-1,`,z+1) = `) ; simplify(Ps(i,z,beta)-b*Ps(i-1,z+1,beta)); printf(`%4s%d%9s%d%8s\n`, `Pss(`,i, `,z)-bPss(`,i-1,`,z+1) = `) ; simplify(Pss(i,z,beta)-b*Pss(i-1,z+1,beta)); od; #------------------------------------------------------------- # Computing the limiting mass ratio alpha #------------------------------------------------------------- printf(`\nLimiting mass ratio alpha \n`); b:= exp(-beta): broot := infinity : printf(`%12s%12s%16s\n`, `tc (T)`, `alpha`,` 1/alpha`) ; for i from 1 to 5 do broot := fsolve(eval(Ps(i,1,beta)+Ps(i-1,1,beta))=0,beta,0..broot) : printf(`%9d%16.6f%16.6f\n`,i+1,broot/2, 2/broot) : od: #--------------------------------------------------------------- # Numerically Computing tc, e and sigma_max as functions of the # mass ratio 1/alpha, for example 1/alpha from 0.5 to 50 # by step 0.5 #-------------------------------------------------------------- # Number of digits carried in floats is set to 8 #-------------------------------------------------------------- Digits := 8 : printf(`\nDuration of conatct, coefficient of restitution and maximum stress \n`); #-------------------------------------------------------------- #The results cabe be saved in the file results.dat #-------------------------------------------------------------- # writeto(`results.dat`) ; #-------------------------------------------------------------- jmax := 20 : # the maximum inverse of the mass ratio considered, jamx = max(1/alpha) jpoints := 50 : # computation points from 0 to jmax j := `j` : printf(`\n alpha 1/alpha tc(T) e sigma_max\n`); for j from 1 by 1 to jpoints do alpha := evalf(jpoints/j/jmax) : beta := 2*alpha : b := exp(-beta) : contactduration(beta) : #------------------------------------------------------------- # printing numerical results #------------------------------------------------------------- printf(`%12.6f%12.6f%12.6f%12.6f%12.6f\n`, alpha,1/alpha, tc,e, stress_max); #------------------------------------------------------------- # printing numerical results to an array AP for plotting #------------------------------------------------------------- AP[j,1] := alpha : AP[j,2] := 1/alpha : AP[j,3] := tc : AP[j,4] := e : AP[j,5] := stress_max #------------------------------------------------------------ od : #------------------------------------------------------------ # writeto(terminal); #------------------------------------------------------------- # Plotting these numerical results #------------------------------------------------------------- with(plots) : #-------------------------------------------------------- # Plotting the duration of contact tc vs. 1/alpha #------------------------------------------------------- plotfig := plot([seq([AP[j,2],AP[j,3]],j=1..jpoints)], thickness = 3, resolution=600): plottex1 := textplot([jmax, 0, `1/a`], font=[SYMBOL,18], align={ABOVE}): plottex2 := textplot([[0.5, 8,'t'], [1,7.9,'c'],[1.5,7.95,`/T`]], font=[TIMES,ROMAN, 18],align={ABOVE,RIGHT}): display({plotfig,plottex1,plottex2}); #-------------------------------------------------------- # Plotting the coefficient of restitution e vs. 1/alpha #------------------------------------------------------- plotfig :=plot([seq([AP[j,2],AP[j,4]],j=1..jpoints)], thickness = 3, resolution=600): plottex1 := textplot([jmax, 0, `1/a`], font=[SYMBOL,18], align={ABOVE}): e := 'e' : plottex2 := textplot([0.5, 1,`e`], font=[TIMES,ROMAN,18],align={ABOVE,RIGHT}): display({plotfig,plottex1,plottex2}); #-------------------------------------------------------- # Plotting the maximum stress sigma_max vs. 1/alpha #------------------------------------------------------- plotfig := plot([seq([AP[j,2],AP[j,5]],j=1..jpoints)], thickness = 3, resolution=600) : plottex1 := textplot([[1.5,7,`s /s`], [jmax, 0, `1/a`]], font=[SYMBOL,18], align={ABOVE}): plottex2 := textplot([2, 6.9,`max 0`], font=[TIMES,ROMAN, 12],align={ABOVE}): display({plotfig,plottex1,plottex2}); #------------------------------------------------------------ # Computing physical quantities for a given system #------------------------------------------------------------ printf(`Physical parameters: \n`); #------------------------------------------------------------- # Material parameter #------------------------------------------------------------- pho := 1 ; Emod := 1 ; #------------------------------------------------------------- # Wave propagation velocity #------------------------------------------------------------- c:= sqrt(Emod/pho) ; #------------------------------------------------------------- # Geometrical data of the rod #------------------------------------------------------------- L := 1 ; A := 1 ; #------------------------------------------------------------- # Mass of the rod #------------------------------------------------------------- m1 := pho*A*L; #------------------------------------------------------------- # Mass of the striking body #------------------------------------------------------------- m2 := 1 ; #------------------------------------------------------------- # Initial velocity of the striking body #------------------------------------------------------------- v2 := 1 ; #------------------------------------------------------------- # Mass ratio #------------------------------------------------------------- alpha := m1/m2 ; beta := 2*alpha : b := exp(-beta) : #-------------------------------------------------------------- # Numerical results about wave propagations in the rod #-------------------------------------------------------------- sigma0 := v2/c*Emod : T := L2/c : contactduration(beta) : ctc := c*tc*T : constant1 := f_contact(ctc)+f_contact(ctc-L2) : #-------------------------------------------------------------- # the maximum stress #-------------------------------------------------------------- stress_max := sigma0*stress_max : lprint(`The maximum stress =`, stress_max ) ; #-------------------------------------------------------------- # The coefficient of restitution #-------------------------------------------------------------- lprint(`The coefficient of restitution = `, e) ; #-------------------------------------------------------------- # The duration of contact #-------------------------------------------------------------- lprint(`The duration of contact =`, tc*T) ; #-------------------------------------------------------------- # Visualizing numerical results #-------------------------------------------------------------- plot(evaln(u(0,t)),t=0..10*T, title=`Displacement at the contact end u(0,t)`, resolution = 600,thickness=3); plot(evaln(v(0,t)),t=0..10*T, title=`Velocity at the contact end v(0,t)`, resolution = 1600,thickness=3); plot(evaln(sigma(L,t)),t=0..10*T, title=`Stress at the fixed end sigma(L,t)`, resolution = 1600,thickness=3); plot3d(evaln(u(x,t)),x=0..L,t=0..10*T, numpoints=1000, grid=[20,20], title = `Displacement of the rod u(x,t)`, style = PATCH, axes = FRAME, orientation=[-145,55] ); plot3d(evaln(sigma(x,t)),x=0..L,t=0..10*T, numpoints=4000, grid=[20,20], title = `Stress of the rod sigma(x,t)`, style = PATCH, axes = FRAME, orientation=[-125,55]); plot3d(evaln(sigma(x,t)),x=0..0.2*L,t=0..2.5*T, numpoints=4000, grid=[20,20], title = `Stress of the rod sigma(x,t)`, style = PATCH, axes = FRAME, orientation=[-155,70] ); #animation of the rod printf(`Animation of the rod \n`); animate(evaln(u(x,t)),x=0..L,t=0..5*T,frames=25, title = `Displacement of the rod u(x,t)`, thickness=6, color=green);