In theory of vibrations, Duhamel's integral is a way of calculating the response of linear systems and structures to arbitrary time-varying external perturbation. It is named after Jean-Marie Duhamel. It is the convolution of a linear system's impulse response function with the input function.
Introduction
Background The response of a linear, viscously-damped, single degree of freedom (SDOF) system to a time-varying mechanical excitation p(t) is given by the following second-order ordinary differential equation
m d 2 x ( t ) d t 2 + c d x ( t ) d t + k x ( t ) = p ( t ) {\displaystyle m{\frac {d^{2}x(t)}{dt^{2}}}+c{\frac {dx(t)}{dt}}+kx(t)=p(t)}
where m is the (equivalent) mass, x stands for the amplitude of vibration, t for time, c for the viscous damping coefficient, and k for the stiffness of the system or structure. If a system initially rests at its equilibrium position, from where it is acted upon by a unit-impulse at the instance t=0, i.e., p(t) in the equation above is a Dirac delta function δ(t), x ( 0 ) = d x d t | t = 0 = 0 {\textstyle x(0)=\left.{\frac {dx}{dt}}\right|_{t=0}=0} , then by solving the differential equation one can get a fundamental solution (known as a unit-impulse response function)
h ( t ) = { 1 m ω d e − ς ω n t sin ω d t , t > 0 0 , t < 0 {\displaystyle h(t)={\begin{cases}{\frac {1}{m\omega _{d}}}e^{-\varsigma \omega _{n}t}\sin \omega _{d}t,&t>0\\0,&t<0\end{cases}}}
where ς = c 2 k m {\textstyle \varsigma ={\frac {c}{2{\sqrt {km}}}}} is called the damping ratio of the system, ω n = k m {\textstyle \omega _{n}={\sqrt {\frac {k}{m}}}} is the natural angular frequency of the undamped system (when c=0) and ω d = ω n 1 − ς 2 {\textstyle \omega _{d}=\omega _{n}{\sqrt {1-\varsigma ^{2}}}} is the angular frequency when damping effect is taken into account (when c ≠ 0 {\displaystyle c\neq 0} ). If the impulse happens at t=τ instead of t=0, i.e. p ( t ) = δ ( t − τ ) {\displaystyle p(t)=\delta (t-\tau )} , the impulse response is
h ( t − τ ) = 1 m ω d e − ς ω n ( t − τ ) sin [ ω d ( t − τ ) ] {\displaystyle h(t-\tau )={\frac {1}{m\omega _{d}}}e^{-\varsigma \omega _{n}(t-\tau )}\sin[\omega _{d}(t-\tau )]} ， t ≥ τ {\displaystyle t\geq \tau }
Conclusion Regarding the arbitrarily varying excitation p(t) as a superposition of a series of impulses:
p ( t ) ≈ ∑ τ < t p ( τ ) ⋅ Δ τ ⋅ δ ( t − τ ) {\displaystyle p(t)\approx \sum _{\tau <t}{p(\tau )\cdot \Delta \tau \cdot \delta }(t-\tau )}
then it is known from the linearity of system that the overall response can also be broken down into the superposition of a series of impulse-responses:
x ( t ) ≈ ∑ τ < t p ( τ ) ⋅ Δ τ ⋅ h ( t − τ ) {\displaystyle x(t)\approx \sum _{\tau <t}{p(\tau )\cdot \Delta \tau \cdot h}(t-\tau )}
Letting Δ τ → 0 {\displaystyle \Delta \tau \to 0} , and replacing the summation by integration, the above equation is strictly valid
x ( t ) = ∫ 0 t p ( τ ) h ( t − τ ) d τ {\displaystyle x(t)=\int _{0}^{t}{p(\tau )h(t-\tau )d\tau }}
Substituting the expression of h(t-τ) into the above equation leads to the general expression of Duhamel's integral
x ( t ) = 1 m ω d ∫ 0 t p ( τ ) e − ς ω n ( t − τ ) sin [ ω d ( t − τ ) ] d τ {\displaystyle x(t)={\frac {1}{m\omega _{d}}}\int _{0}^{t}{p(\tau )e^{-\varsigma \omega _{n}(t-\tau )}\sin[\omega _{d}(t-\tau )]d\tau }}
Mathematical proof The above SDOF dynamic equilibrium equation in the case p(t)=0 is the homogeneous equation:
d 2 x ( t ) d t 2 + c ¯ d x ( t ) d t + k ¯ x ( t ) = 0 {\displaystyle {\frac {d^{2}x(t)}{dt^{2}}}+{\bar {c}}{\frac {dx(t)}{dt}}+{\bar {k}}x(t)=0} , where c ¯ = c m , k ¯ = k m {\displaystyle {\bar {c}}={\frac {c}{m}},{\bar {k}}={\frac {k}{m}}}
The solution of this equation is:
x h ( t ) = C 1 ⋅ e − 1 2 ( c ¯ + c ¯ 2 − 4 ⋅ k ¯ ) t + C 2 ⋅ e 1 2 ( − c ¯ + c ¯ 2 − 4 ⋅ k ¯ ) t {\displaystyle x_{h}(t)=C_{1}\cdot e^{-{\frac {1}{2}}\left({\bar {c}}+{\sqrt {{\bar {c}}^{2}-4\cdot {\bar {k}}}}\right)t}+C_{2}\cdot e^{{\frac {1}{2}}\left(-{\bar {c}}+{\sqrt {{\bar {c}}^{2}-4\cdot {\bar {k}}}}\right)t}}
The substitution: A = 1 2 ( c ¯ − c ¯ 2 − 4 k ¯ ) , B = 1 2 ( c ¯ + c ¯ 2 − 4 k ¯ ) , P = c ¯ 2 − 4 k ¯ , P = B − A {\displaystyle A={\frac {1}{2}}\left({\bar {c}}-{\sqrt {{\bar {c}}^{2}-4{\bar {k}}}}\right),\;B={\frac {1}{2}}\left({\bar {c}}+{\sqrt {{\bar {c}}^{2}-4{\bar {k}}}}\right),\;P={\sqrt {{\bar {c}}^{2}-4{\bar {k}}}},\;P=B-A} leads to:
x h ( t ) = C 1 e − B ⋅ t + C 2 e − A ⋅ t {\displaystyle x_{h}(t)=C_{1}e^{-B\cdot t}\;+\;C_{2}e^{-A\cdot t}}
One partial solution of the non-homogeneous equation: d 2 x ( t ) d t 2 + c ¯ d x ( t ) d t + k ¯ x ( t ) = p ¯ ( t ) {\textstyle {\frac {d^{2}x(t)}{dt^{2}}}+{\bar {c}}{\frac {dx(t)}{dt}}+{\bar {k}}x(t)={\bar {p}}(t)} , where p ¯ ( t ) = p ( t ) m {\textstyle {\bar {p}}(t)={\frac {p(t)}{m}}} , could be obtained by the Lagrangian method for deriving partial solution of non-homogeneous ordinary differential equations. This solution has the form:
x p ( t ) = ∫ p ( t ) ¯ ⋅ e A t d t ⋅ e − A t − ∫ p ( t ) ¯ ⋅ e B t d t ⋅ e − B t P {\displaystyle x_{p}(t)={\frac {\int {{\bar {p(t)}}\cdot e^{At}dt}\cdot e^{-At}-\int {{\bar {p(t)}}\cdot e^{Bt}dt}\cdot e^{-Bt}}{P}}}
Now substituting: ∫ p ( t ) ¯ ⋅ e A t d t | t = z = Q z , ∫ p ( t ) ¯ ⋅ e B t d t | t = z = R z {\textstyle \left.\int {{\bar {p(t)}}\cdot e^{At}dt}\right|_{t=z}=Q_{z},\left.\int {{\bar {p(t)}}\cdot e^{Bt}dt}\right|_{t=z}=R_{z}} ,where ∫ x ( t ) d t | t = z {\textstyle \left.\int x(t)dt\right|_{t=z}} is the primitive of x(t) computed at t=z, in the case z=t this integral is the primitive itself, yields:
x p ( t ) = Q t ⋅ e − A t − R t ⋅ e − B t P {\displaystyle x_{p}(t)={\frac {Q_{t}\cdot e^{-At}-R_{t}\cdot e^{-Bt}}{P}}}
Finally the general solution of the above non-homogeneous equation is represented as:
x ( t ) = x h ( t ) + x p ( t ) = C 1 ⋅ e − B t + C 2 ⋅ e − A t + Q t ⋅ e − A t − R t ⋅ e − B t P {\displaystyle x(t)=x_{h}(t)+x_{p}(t)=C_{1}\cdot e^{-Bt}+C_{2}\cdot e^{-At}+{\frac {Q_{t}\cdot e^{-At}-R_{t}\cdot e^{-Bt}}{P}}}
with time derivative:
d x d t = − A e − A t ⋅ C 2 − B e − B t ⋅ C 1 + 1 P [ Q t ˙ ⋅ e − A t − A Q t ⋅ e − A t − R t ˙ ⋅ e − B t + B R t ⋅ e − B t ] {\displaystyle {\frac {dx}{dt}}=-Ae^{-At}\cdot C_{2}-Be^{-Bt}\cdot C_{1}+{\frac {1}{P}}\left[{\dot {Q_{t}}}\cdot e^{-At}-AQ_{t}\cdot e^{-At}-{\dot {R_{t}}}\cdot e^{-Bt}+BR_{t}\cdot e^{-Bt}\right]} , where Q t ˙ = p ( t ) ⋅ e A t , R t ˙ = p ( t ) ⋅ e B t {\displaystyle {\dot {Q_{t}}}=p(t)\cdot e^{At},{\dot {R_{t}}}=p(t)\cdot e^{Bt}}
In order to find the unknown constants C 1 , C 2 {\displaystyle C_{1},C_{2}} , zero initial conditions will be applied:
x ( t ) | t = 0 = 0 : C 1 + C 2 + Q 0 ⋅ 1 − R 0 ⋅ 1 P = 0 {\displaystyle x(t)|_{t=0}=0:C_{1}+C_{2}+{\frac {Q_{0}\cdot 1-R_{0}\cdot 1}{P}}=0} ⇒ C 1 + C 2 = R 0 − Q 0 P {\displaystyle C_{1}+C_{2}={\frac {R_{0}-Q_{0}}{P}}}
d x d t | t = 0 = 0 : − A ⋅ C 2 − B ⋅ C 1 + 1 P ⋅ [ − A ⋅ Q 0 + B ⋅ R 0 ] = 0 {\displaystyle \left.{\frac {dx}{dt}}\right|_{t=0}=0:-A\cdot C_{2}-B\cdot C_{1}+{\frac {1}{P}}\cdot [-A\cdot Q_{0}+B\cdot R_{0}]=0} ⇒ A ⋅ C 2 + B ⋅ C 1 = 1 P ⋅ [ B ⋅ R 0 − A ⋅ Q 0 ] {\displaystyle A\cdot C_{2}+B\cdot C_{1}={\frac {1}{P}}\cdot [B\cdot R_{0}-A\cdot Q_{0}]}
Now combining both initial conditions together, the next system of equations is observed:
C 1 + C 2 = R 0 − Q 0 P B ⋅ C 1 + A ⋅ C 2 = 1 P ⋅ [ B ⋅ R 0 − A ⋅ Q 0 ] | C 1 = R 0 P C 2 = − Q 0 P {\displaystyle \left.{\begin{alignedat}{5}C_{1}&&\;+&&\;C_{2}&&\;=&&\;{\frac {R_{0}-Q_{0}}{P}}&\\B\cdot C_{1}&&\;+&&\;A\cdot C_{2}&&\;=&&\;{\frac {1}{P}}\cdot [B\cdot R_{0}-A\cdot Q_{0}]\end{alignedat}}\right|{\begin{alignedat}{5}C_{1}&&\;=&&\;{\frac {R_{0}}{P}}&\\C_{2}&&\;=&&\;-{\frac {Q_{0}}{P}}\end{alignedat}}}
The back substitution of the constants C 1 {\displaystyle C_{1}} and C 2 {\displaystyle C_{2}} into the above expression for x(t) yields:
x ( t ) = Q t − Q 0 P ⋅ e − A ⋅ t − R t − R 0 P ⋅ e − B ⋅ t {\displaystyle x(t)={\frac {Q_{t}-Q_{0}}{P}}\cdot e^{-A\cdot t}-{\frac {R_{t}-R_{0}}{P}}\cdot e^{-B\cdot t}}
Replacing Q t − Q 0 {\displaystyle Q_{t}-Q_{0}} and R t − R 0 {\displaystyle R_{t}-R_{0}} (the difference between the primitives at t=t and t=0) with definite integrals (by another variable τ) will reveal the general solution with zero initial conditions, namely:
x ( t ) = 1 P ⋅ [ ∫ 0 t p ¯ ( τ ) ⋅ e A τ d τ ⋅ e − A t − ∫ 0 t p ¯ ( τ ) ⋅ e B τ d τ ⋅ e − B t ] {\displaystyle x(t)={\frac {1}{P}}\cdot \left[\int _{0}^{t}{{\bar {p}}(\tau )\cdot e^{A\tau }d\tau }\cdot e^{-At}-\int _{0}^{t}{{\bar {p}}(\tau )\cdot e^{B\tau }d\tau }\cdot e^{-Bt}\right]}
Finally substituting c = 2 ξ ω m , k = ω 2 m {\displaystyle c=2\xi \omega m,\;k=\omega ^{2}m} , accordingly c ¯ = 2 ξ ω , k ¯ = ω 2 {\displaystyle {\bar {c}}=2\xi \omega ,{\bar {k}}=\omega ^{2}} , where ξ<1 yields:
P = 2 ω D i , A = ξ ω − ω D i , B = ξ ω + ω D i {\displaystyle P=2\omega _{D}i,\;A=\xi \omega -\omega _{D}i,\;B=\xi \omega +\omega _{D}i} , where