Summer Research Fellowship Programme of India's Science Academies

Wave Motion in an Elastically Connected Double Beam using Super-convergent Finite Element Formulation

Himanshu Giria

Indian Institute of Engineering Science and Technology, Shibpur, India

Guided by:

S Gopalakrishnan

Indian Institute of Science, Bangalore, India


Composites form an integral part of the materials used in aircrafts; gradually replacing the conventional use of metals. Adhesively bonded metal/composite joints are common nowadays, and to check if they can sustain a certain number of flight hours without failing becomes important as they differ from the conventional modes of failure, like buckling, yielding, etc. Wave propagation is used here to study bonded joints. The mathematical model consists of a double beam element with continuous elasticity (stiffness, say, K) in between the two Timoshenko beam-rod elements, capable of deforming in axial (rod-like), transverse (Euler Bernoulli beam-like) and shear deformations (Timoshenko beam speciality). This yields six coupled differential equations, solved to get the deformations in terms of 12 arbitrary constants. These expressions, which cotaine both algebraic and exponential terms, are used to obtain the mass and stiffness matrices, and the expressions being exact give more accurate stiffness and mass matrices as compared to the conventional Finite Element Method. These matrices would be used to do the analysis of of high frequency (ultrasonic) wave propagation. Time marching using the Newmark method is the next thing in the study. It is to be noted that different values of K can be used to model different levels of bond strengths.

Keywords: Wave Propagation, Timoshenko beam-bar, Newmark-Beta scheme


FEMFinite Element Method
FRPFiber Reinforced Polymer
ODEOrdinary Differential Equation
SFEMSpectral Finite Element Method
FEAFinite Element Analysis
uaxial displacement in beam
`wtransverse displacement
φrotation of plane normal u (axial) and w(transvers)displacements of the beam



The study intends to find out the behaviour of adhesively bonded metallic/composite joints by modeling them with a set of two beams, with the bond represented by continuous distribution of springs. Composites are emerging as important materials not only for aerospace applications but also for other fields, due to their distinct advantages like corrosion resistance, higher strength to weight ratio, toughness, etc. This property is useful in aerospace applications, as being light-weight is a great advantage. The reason for this property (hgh strength to weight and stiffness to weight ratios) of FRP composites might be due to the fact the fibers are oriented along a particular direction, due to which, the strength in the fiber direction is maximum, and it is lesser in other directions. In other words, along the fiber axis, the strength is maximum but in any other direction, the strength and stiffness have very less values.However, in many cases the application may demand stiffness in only some specified direction, where we can use composites effectively.  

Composites are used sometimes in combination with metals. In designing a structure, two different entities, like two beams, may need to be joined. In conventional metals, riveting was an option. Today, with composites, it is not much suited as is evident by its lessened use. The common joints today have an adhesive bond in them. So, the necessity to study the physics of such joints is indispensable, especially when safety is needed.

The mode of failure is not the same, as yielding-in-tension and such modes. The wave propagation study is a tool to study these joints, and this is the main aim of the study. The waves used would be having a low wavelength, because loads like sudden impact provide a loading of very high frequency. So conventional FEM would need a very small mesh size, which is not a good idea to go on with.

Finite Element Method

Finite Element Method is simply a numerical method, like Finite DIfference Method, which converts differential equations into algebraic equations. It is used to solve problems of the mathematical model of systems, which may be a mechanical system comprising beams and frames, a thermodynamic system dealing with heat flows, a chemical system comprising chemical reactions, a system in the field of biotechnology, or any other such field. ​​

The need for FEA comes from the limitation of theoretical/analytical methods to include, non-uniform / other complicated details of a complex system, such as, fuselage of an aircraft, wings, etc. The models for using theoretical methods are solved with simplifications with a great deal of assumptions. For example, in analyses of stress concentrations at the windows the fuselage is assumed cylindrical, now the closeness of its shape to a cylinder is indeed questionable. This causes FEA to evolve as it does not use so many simplifications. It uses the method of dividing the structure into small parts. However, loosely speaking, just like any general conservation principle, net effort seems to be conserved as with FEA, the complexity of the solution proceedure is increased and it is because of the modern day fast, efficient and cost-effective computing that this field of numerical computations, which is used by FEM, got a chance to boom.  Today, FEM seems to have become an integral part of analysing systems and structures in a wide diversity of manufacturing and research industries.

The Finite Difference Method, which is based on the Taylor series, is somewhat inaccurate as the error of truncating the terms of Taylor series continues to accumuate and deviates from the correct solution as the complexity of the domain becomes high. This is why FEM came into limelight and is used nowadays so much to places where it can be used. However, it is not used for some cases, for example, say on the Navier Stokes' Equation.

The Euler-Bernoulli beam, say has a differential equation of the form:

d2dx2(EI(x)d2wdx2)=q=theloading/unitlength \frac{d^2}{dx^2}(EI(x)\frac{d^2w}{dx^2}) =q = the\; loading/unit \;length

This in general can be solved if EI(x) is not a non-linear function in x; but if say it is non-linear, then it becomes difficult to solve it analytically, and so what comes out is the idea of using some numerical technique, like FEM, which converts it to an algebraic equation which can then be solved easily.

The actual working of FEM relies on a simple technique of solving for a few finite number of points, instead of the infinite space, and then using simple polynomial functions to interpolate between the points (nodes). Depending on the order of interpolation, the obtained solution may be made closer and nearer to the exact solution, which may be at the cost of increased complexity in the solution.

The properties of stiffness matrix are:

  • Symmetric - This arises due to the modelling assumption that the loads are linearly related to the displacements, which can be proven by Betti- Maxwell theorem.
  • Sparsity - A "sparse" matrix contains a large number of zeros, in the off-diagonal terms. This is reasonable here, because if the node numbering is done in a particular way, the chances that a force applied at node on one side may cause an effect (displacement) on a far-off node so elements is very less (corresponding to those in the corners of the op-right and bottom-left regions of the matrix).
  • Non-negative diagonal: The diagonal terms signify the effect of force and displacement on the same node, and it seems quite unphysical that a leftward force could cuse a right-ward displacement at that point.

Wave Propagation

The study of wave propagation takes a different route because of the impracticability of conventional FEM, which otherwise, is a powerful tool for other structural mechanics problems (and also heat transfer problems and many more as discussed above), to solve for this wave propagation behaviour, especially if the frequency of the wave which is being dealt with is very high (in the study, ultrasonic range). This is because higher frequency implies smaller wavelengths and to capture the effects of small wavelengths, the mesh size needs to be smaller than even a millimetre by many orders of magnitude, which is not an intelligent thing to do. So, we need to shift our analysis in the frequency domain (spectral domain). ​​

The wave input may be generated by a machine; other forms of high frequency wave generation in a material structure includes impact.

Beam Theory

Timoshenko Beam :

The assumptions in this model are lessened with respect to those for Euler- Bernoulli beam, i.e., the plane sections are still assumed to remain plane but not neessarily normal to the longitudinal axis of the beam. This increases the number of equations and unknowns i.e., we now take into account shear deformations of the beam along with flexural deformation, the latter being the only one considered in an Euler- Bernoulli model of the beam. The assumption of constant shear stress across the cross-section​ ​leads to the necessity of a shear correction factor to give results with better accuracy. In our case, we have chosen a rectangular cross-sectioned beam which has a shear correction factor γ close to 5/6.



The main motivation for the work is obtained from ​the textbook by Professor ​​Srinivasan Gopalakrishnan, 2016​. This shows the importance and necessity of wave propagation analysis, starting from ground level, and taking it to a height which very vew are able to grasp, in a smooth, interesting, and yet simple way.

The next thing in study is about the materials used nowadays, namely, composites, and in the book by ​Robert M. Jones 1999​​​, the theory of composites, starting from the answers to the basic questions about what and why of composites, have been presented in a systematic and easy to understand way.

Moving on to the method of study, we have FEM in the base, to which the textbook by ​​​J. N. Reddy​​​ gives a good introduction and idea about the Finite Element Method, and later also proceeds, to some practical applications, not only in solid mechanics, but also in other problems, like heat flows. The book by ​​Klaus-Jürgen Bathe, 1996​​ sets up a way for studying Finite Element Method, depending upon the intention of the seeker, and hence is quite useful in orienting oneself in the right direction. It also provides the algorithm for the Newmark-Beta scheme. A brief detailing of Finite Element Method could be found in the book by ​​Robert D. Cook, et al​​, which, is a good choice to go on with, if mathematical side of FEM is the concern. For studying about the models of beams, rods and also for knowing something more about spectral FEM, the book by ​​James F. Doyle​​ may be used.

The book by ​​​​V. Rajaraman, 1997​​ is helpful in understanding the syntax of FORTRAN 77 with some splendid examples to deepen the understanding of programming in general and FORTRAN77 in particular.


The mathematical representation of various physical phenomena involves differential equations in most cases. In this study, we need to solve a set of coupled differential equations. A coupled system is one in which the dependent variables are present together in a differential equation. For example:

aẍ+by = 0

c ÿ + kx = 0

while an example of uncoupled set of ODEs would be

aẍ+bx = 0

c ÿ + ky = 0

where ẍ = d2xdt2\frac{d^{2^{}}x}{dt^2}  and ÿ=.​  d2ydt2    \frac{d^{2}y}{dt^2}    

If a,b,c,k be constants, then the above set is said to be a linear ODE, and if they are functions of x and/or y they become non-linear.

It is clearly visible how difficult it is to solve a set of coupled differential equations. An approach of eliminating variables leads to an increase in the order of derivatives, sometimes causing some extra constants, which need to be lessened to the d.of. available, by some other means. So, finally in our case, we get the final solution as shown next, without the detailing of the unscuccessfully tried methods of solution.

Steps to Solution

The governing differential equations for a set of two Timoshenko beam-rod elements​​, connected by continuous elastic springs of stiffness K=KF, are shown below:

        A11t  d2wtdx2    B11t  d2φtdx2  =  0\displaystyle \;\;\;\;A_{11}^t\;\frac{d^2w^t}{dx^2}\;-\;B_{11}^t\;\frac{d^2\varphi^t}{dx^2}\;=\;0
    B11td2utdx2  +  D11td2φtdx2  +  γA55t  (dwtdx    φ)  =  0\displaystyle \;\;-B_{11}^t\frac{d^2u^t}{dx^2}\;+\;D_{11}^t\frac{d^2\varphi^t}{dx^2}\;+\;\gamma A_{55}^t\;(\frac{d^{}w^t}{dx^{}}\;-\;\varphi)\;=\;0
  γAt55  (  d2wtdx2    dφtdx)  +  KF(wb    wt)  =  0\displaystyle \;\gamma{A^t}_{55}\;(\;\frac{d^2w^t}{dx^2}\;-\;\frac{d\varphi^t}{dx})\;+\;KF(w^b\;-\;w^t)\;=\;0

The corresponding equations for lower beam may be obtained by substituting b in place of the t's.

The meaning of various symbols may be found in Appendix A

Upon eliminating the terms containing the double derivatives of axial displacement 'u' we get left out with 4 equations among w(transverse displacement) and φ (shear displacement) for each beam, i.e., 4 equations and 4 unknowns, down from 6 equations with 6 unknowns..

This set of 4 equations shown below, is now same as that of Timoshenko beam equations:

Equations for Timoshenko Beam:

G1A1γ1(d2w1dx2dϕ1dx)+K(w1w2)=0,EId2ϕ1dx2+G1A1γ1(dw1dxϕ1)=0,subscript1forupperbeamG2A2γ2(d2w2dx2dϕ2dx)+K(w2w1)=0,EId2ϕ2dx2+G2A2γ2(dw2dxϕ2)=0,subscript2forlowerbeam \begin{array}{l}\begin{array}{l}\begin{array}{l}G_1A_1\gamma_1\;(\;\frac{d^2w_1}{dx^2}\;-\;\frac{d\phi_1}{dx}\;)\;+\;K\;(w_{1\;}-\;w_{2\;})=0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;,\;\;\\EI\frac{d^2\phi_1}{dx^2}\;+\;G_1A_1\gamma_1\;(\frac{\;d^{}w_1}{dx^{}}\;-\phi_1\;)\;=\;0,\;subscript\;_1\;for\;upper\;beam\;\;\end{array}\\G_2A_2\gamma_2\;(\;\frac{d^2w_2}{dx^2}\;-\;\frac{d\phi_2}{dx}\;)\;+\;K\;(\;w_{2\;}-\;w_{1\;})=0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;,\;\;\\EI\frac{d^2\phi_2}{dx^2}\;+\;G_2A_2\gamma_2\;(\frac{\;d^{}w_2}{dx^{}}\;-\phi_2\;)\;=\;0,\;subscript\;_2\;for\;lower\;beam\;\;\end{array}\\\end{array}

Equations in the new notations:

(Dt11-Bt112At11)d2φtdx2+γAt55(dwtdx-φt) =  0
γAt55(d2wtdx2dφtdx  +  KF(wbwt)  =    0\displaystyle \begin{array}{l}\gamma{A^t}_{55}(\frac{d^2w^t}{dx^2}-\frac{d\varphi^t}{dx}\;+\;KF(w^b-w^t)\;=\;\;0\\\end{array}
(Db11 - Bb112Ab11)d2φbdx2+γAb55(dwbdx-φb) =  0
γAb55(d2wbdx2dφbdx  +  KF(wtwb)  =    0\displaystyle \gamma{A^b}_{55}(\frac{d^2w^b}{dx^2}-\frac{d\varphi^b}{dx}\;+\;KF(w^t-w^b)\;=\;\;0

These equations as shown above, may be solved in the same way as the Timoshenko beam

Differentiating equations (5) and (7) and re-arranging them, and then substituting them in equations(6) and (8) we get:

pt  d3φtdx3    pb  d3φbdx3=0-p^t\;\frac{d^3\varphi^t}{dx^3}\;-\;p^b\;\frac{d^3\varphi^b}{dx^3}=0

which on integrating 3 times with respect to x yields:​

          φt=  pbptφb+a2x22+a1x+a0\displaystyle \;\;\;\;\;\varphi^t=\;\frac{-p^b}{p^t}\varphi^b+a_2\frac{x^2}2+a_1x+a_0

where a's denote the 3 constants.

From (4) and (6)

  dwtdx  =    φt      ptγAt55  d2φtdx2    =    φt    std2φtdx2  \displaystyle \;\frac{dw^t}{dx}\;=\;\;\varphi^t\;-\;\frac{\;p^t}{\gamma{A^t}_{55}}\;\frac{\displaystyle d^2\varphi^t}{\displaystyle dx^{2\;}}\;=\;\;\varphi^t\;-\;s^t\frac{\displaystyle d^2\varphi^t}{\displaystyle dx^{2\;}}
      dwbdx  =    φb      pbγAb55  d2φbdx2    =    φb    sbd2φbdx2  \displaystyle \;\;\frac{\;dw^b}{dx}\;=\;\;\varphi^b\;-\;\frac{\;p^b}{\gamma{A^b}_{55}}\;\frac{\displaystyle d^2\varphi^b}{\displaystyle dx^{2\;}}\;=\;\;\varphi^b\;-\;s^b\frac{\displaystyle d^2\varphi^b}{\displaystyle dx^{2\;}}

Differentiating (7) & using the above 3 equations, yields, on simplification :

d4φbdx4    +  (QtQb)d2φbdx2  +  (Ot+Ob)φb  =  Oba2x22+Oba1x+(Oba0    Obsta2)\displaystyle \frac{d^4\varphi^b}{dx^{4\;}}\;+\;(-Q^t-Q^b)\frac{\displaystyle d^2\varphi^b}{\displaystyle dx^2}\;+\;(O^t+O^b)\varphi^b\;=\;O^ba_2\frac{x^2}2+O^ba_1x+(O^ba_0\;-\;O^bs^ta_2)

Now, this may easily yield φb\varphi^bas below, from which φt\varphi^tfound out .

φb=c1em1x+  c2em2x+  c3em1x  +  c4em2x    +β3β2a1x+β32β2a2x2\displaystyle \varphi^b=c_1e^{m_1x}+\;c_2e^{m_2x}+\;c_3e^{-m_1x}\;+\;c_4e^{-m_2x}\;\;+\frac{\beta_3}{\beta_2}a_1x+\frac{\beta_3}{2\beta_2}a_2x^2

By substitution in equations above, 'w', the component of transverse displacements, normal to the longitudinal axis of the beam may be found to be :

wb    =[1m1    sbm1]c1em1x  +  [1m2    sbm2]c2em2x+[1m1  +  sbm1]c1em1x+  [1m2  +  sbm2]c2em2x+[β3β2a0+  β4β2a2  β1β3β22a2+sbβ3β2a2]  x  +    β32β2a1x2+  β36β2a2x3+c5\displaystyle w^{b\;}\;=\lbrack\frac1{m_1}\;-\;s^bm_1\rbrack c_1e^{m_1x}\;+\;\lbrack\frac1{m_2}\;-\;s^bm_2\rbrack c_2e^{m_2x}+\lbrack-\frac1{m_1}\;+\;s^bm_1\rbrack c_1e^{-m_1x}+\;\lbrack-\frac1{m_2}\;+\;s^bm_2\rbrack c_2e^{-m_2x}+\lbrack\frac{\beta_3}{\beta_2}a_0+\;\frac{\beta_4}{\beta_2}a_2-\;\frac{\beta_1\beta_3}{\beta_2^2}a_2+s^b\frac{\beta_3}{\beta_2}a_2\rbrack\;x\;+\;\;\frac{\beta_3}{2\beta_2}a_1x^2+\;\frac{\beta_3}{6\beta_2}a_2x^3+c_5

This can be substituted back to the original governing equations, to get the axial components of the displacements u for each beam to be :

ub=  Bb11Ab11[c1em1x+  c2em2x+  c3em1x  +  c4em2x    +β32β2a2x2]  +  a3x+a4\displaystyle u^b=\;\frac{{B^b}_{11}}{{A^b}_{11}}\lbrack c_1e^{m_1x}+\;c_2e^{m_2x}+\;c_3e^{-m_1x}\;+\;c_4e^{-m_2x}\;\;+\frac{\beta_3}{2\beta_2}a_2x^2\rbrack\;+\;a_3x+a_4
  ut=  Bt11At11[c1em1x+c2em2x+c3em1x+c4em2x  ](pbpt)  +Bt11At11[(pbpt)β32β2+12]+a5x+a6    \displaystyle \;u^t=\;\frac{{B^t}_{11}}{{A^t}_{11}}\lbrack c_1e^{m_1x}+c_2e^{m_2x}+c_3e^{-m_1x}+c_4e^{-m_2x}\;\rbrack(-\frac{p^b}{p^t})\;+\frac{{B^t}_{11}}{{A^t}_{11}}\lbrack(-\frac{p^b}{p^t})\frac{\beta_3}{2\beta_2}+\frac12\rbrack+a_5x+a_6\;\;

Discussion (Nature of equations)

The resulting solution of the unknowns contained both algebraic polynomial and exponential parts, since they represented the exact solutions for the model instead of only polyniomials, which is generally the case with the conventional Finite Element Method. This exactness of solution makes this method more accurate than the FEM in terms of stiffness & mass matrices. This is the novelty of this method, and despite being simple, it gives solutions ​that are close to exact. The inaccuracy may still be present due to the following:

  • Error in domain modelling
  • Numerical errors

The third error that is not involvd in this method is

3. Approximation of solution.

Hence, the name Super-convergent as it quickly converges to the correct solution.

It is to be noted that we had considered three degrees of freedom at each node and then formulated various matrices as below:

{u}e12  x  1  =  [A]  12  x  12  {constant}12  x  1\{u{\}^{e_{}}}_{12\;x\;1}\;=\;\lbrack A\rbrack_{\;12\;x\;12}\;\{cons\tan t\}_{12\;x\;1}

{const}  =  [A]1{u}e6  x  1\displaystyle \Rightarrow\{const\}\;=\;\lbrack A\rbrack^{-1}\{u^{}{\}^e}_{6\;x\;1}

where {constant} = [c1 c2 c3 c4 c5 a0 a1 a2 a3 a4 a5 a6]' , and,

[A] is a 12 x 12 sized mat. calculated from the equations (13)to (15) at the nodal points and ,

{U} has displacement components {U} =[ u1  w1  φ1  u2  w2  φ2  u3  w3  φ3  u4  w4  φ4u_1\;w_1\;\varphi_1\;u_2\;w_2\;\varphi_2\;u_3\;w_3\;\varphi_3\;u_4\;w_4\;\varphi_4  ]' where the node numbering is as explained below :

    Schematic of model with node no.

    Then, we calculated the D(6 x 12) matrix which is simply matrix representation of deformations at given x location along the length of the beam.

    {U}6  x  1  =  [D]6  x  12  {const}\{U\}_{6\;x\;1}\;=\;\lbrack D\rbrack _{6\;x\;12}\; \{const\}  ​

    So, {U}6  x  1  =  [D]6  x  12  [A]112  x  12{u}e12  x  1\{U\}_{6\;x\;1}\;=\;\lbrack D\rbrack_{6\;x\;12\;}\lbrack A{\rbrack^{-1}}_{{}^{}12\;x\;12}^{}\{u{\}^e}_{12\;x\;1}

    The shape function matrix is given now as by [N]  =  [D][A]1\lbrack N\rbrack\;=\;\lbrack D\rbrack\lbrack A\rbrack^{-1}

    The strain components [ε]  =  []{U} [\varepsilon]\;=\;\lbrack \partial] \{U\}       

    where the [∂] is given by [∂] =  [x0zx0000x1000000x0zx0000x1]    \begin{bmatrix} \frac{\partial }{\partial x} &0&-z \frac{\partial }{\partial x} &0&0&0\\0& \frac{\partial }{\partial x} &-1&0&0&0\\0&0&0& \frac{\partial }{\partial x} &0& -z\frac{\partial }{\partial x} \\0&0&0&0& \frac{\partial }{\partial x} &-1\end{bmatrix}  

    Then, [ε]  =  [][N]{u}e    =[B]{u}e  \lbrack\varepsilon\rbrack\;=\;\lbrack\partial\rbrack\lbrack N\rbrack\{u\}^e\;\; = [B] \{u\}^e\;  ,

    where [B] = [∂D] [A]1\lbrack A\rbrack ^{-1} and

    [D]=  []  [D]\lbrack\partial D\rbrack=\;\lbrack\partial\rbrack\;\lbrack D\rbrack

    Now, stiffness mat.

    [K] = [B][C][B]  b  dx  dz  =  [A]T[D]T[C][D][A]1  b  dx  dz\begin{array}{l}\iint\lbrack B\rbrack'\lbrack C\rbrack\lbrack B\rbrack\;b\;dx\;dz\;\\=\;\iint\lbrack A\rbrack^{-T}\lbrack\partial D\rbrack^T\lbrack C\rbrack\lbrack\partial D\rbrack\lbrack A\rbrack^{-1}\;b\;dx\;dz\end{array}

    [C] is th material stiffness diagonal matrix whose diagonal terms are [E1G1E2G2] [E_1\; G_1\; E_2\; G_2] and similarly, we have mass matrix

    [M]  =    [N][S][N]  b  dx  dz\lbrack M\rbrack\;=\;\iint\;\lbrack N\rbrack'\lbrack S\rbrack\lbrack N\rbrack\;b\;dx\;dz

    [S] is diagonal density matrix of size 6 x 6 whose diagonal elements have ρ_t in the 1st 3 elements & ρ_b for the remaining 3 diagonal terms of the leading diagonal.

    Different codes were written in MATLAB to calculate the intermediate matrices, and also to check if the A matrix is non-singular. Symbolic differentiation and integration were used from MATLAB's symbolic toolbox.

    These integrals were performed using MATLAB and the results obtained were obeying the supposition of symmetricand non-negative real parts in the diagonal for both the stiffness and mass matrices.

    Time marching using Newmark-Beta method​

    This is the stage where dynamic analysis begins. So far we were creating the ground for the actual stud to begin, since this is where we could use our obtained mass and stiffness matrices to analyse the wave propagating through this model and then make an equal/equivalent idea of how it would be in the actual case of an adhesively bonded metal/composite joint.

    Here, we need to numerically integrate using the Newmark-Beta method, which has already been implemented in FORTRAN. In order to first understand the code, a learning of this language was carried out and a sample code generated which would be a part of the main code is shown below:

    1. c This is to solve system of eqns in NN unknowns
    2. c
    3. DIMENSION A(4,5)
    4. DIMENSION X(4)
    5. READ *,NN
    6. READ *,((A(I,J),J=1,NN+1),I=1,NN)
    7. DO 2 M=2,Nn
    8. DO 1 N=1,M-1
    9. IF(A(M,N).NE.0.)THEN
    10. IF(A(M-1,N).NE.0) THEN
    11. P =- A(M,N)/A(M-1,N)
    12. L=1
    13. ELSE
    14. P =- A(M,N)/A(M-2,N)
    15. L=2
    16. ENDIF
    17. pRint *, 'row', m-l
    18. DO 3 J=1,NN+1
    19. A(M,J) = A(M,J) + P*a(M-L,J)
    20. 3 CONTINUE
    21. ENDIF
    22. 1 CONTINUE
    23. pRint *,'U.T.',((A(I,J),J=1,NN+1),' \n ',I=1,NN)
    24. 2 CONTINUE
    25. pRint *,'U.T.',((A(I,J),J=1,NN+1),'\n ',I=1,NN)
    27. DO 4 I=NN,1,-1
    28. Q = 0.
    29. DO 5 J = NN,I+1,-1
    30. Q = Q+ X(J)*a(I,J)
    31. 5 CONTINUE
    32. X(I) = (A(I,NN+1) - Q) / A(I,I)
    33. 4 CONTINUE
    34. DO 6 I=1,NN
    35. PRINT *, 'X(',I,') = ',X(I)
    36. 6 CONTINUE
    37. STOP
    38. END
    Code in FORTRAN77 to solvesystem of equations

    The code solves a system of linear algebraic equations using :

    • Gauss-Elimination technique to get the Upper Trialngular matrix which has all elements below the leadimbg diagonal as 0.
    • Back substitution to solve the unkowns, begining at the last equation.

    It is to be noted that though the above code is valid in simple cases, pivoting has to be taken into account, as it then becomes robust.

    We next plan to implement the Newmark-Beta algorithm using MATLAB, which would no doubt, make the code look simpler and elegant but at the cost of using in-built functions in MATLAB, like a simple command A\b is equivalent to a dozens of FORTRAN lines.

    Conclusion and Recommendations


    The stiffness and mass matrices obtained matched with the theoretical reqirement that it should be symmetric with the diagonal terms having non-negative reap parts. This may serve as a check that our calculations may be right. After this the time response could be seen using only a single element, and later can be tried using more elements, say three or four.

    The future work may be carried out in these areas:

    • Experimental validation of the results
    • Setting up a new mathematical model of the physical situation

    The potential application in crack detection in adhesively bonded joints, may be useful in the aerospace industry. Other industries employing such mechanisms may also find it useful. As time elapses, the strength of the adhesive lessens and it may reflect on our model through the stiffness value of the spring KF. The number of elements may be increased from only 1 here with 4 nodes if the adhesive has different strengths along its length.

    Thus, the usefulness of this ground work may be developed through further work after which the maintenance cost can be lessened as computations can be done on a wide range of frequencies and this may, by the judicious use of the above method, increase the life cycles before going out of fail safe region, thus leading to an increased gap between two successive experimental maintenance checks, which nowadays, are being done routinely for aircrafts.


    Appendix A

    p_t D11t-B11t2A11t

    p_b D11b-B11b2A11b  

    s_t p_tγAt55\frac{p\_t}{\gamma {A^t}_{55}}

    Q_t KFγAt55\frac{KF}{\gamma {A^t}_{55}}

    O_t KF/p_t

    O_b KF/p_b

    β1\beta_1   -Q_t - Q_b

     β2    \beta_2   O_T - O_b

     β3   \beta_3    O_b

     β4  \beta_4   -O_b * s_t

    B11                                                    b(z22z12)E2(1ν2)    B_{11}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{b(z_2^2-z_1^2)E}{2(1-\nu^2)}\;\;        

     A11                                                    b(z2z1)E(1ν2)  \displaystyle   A_{11}\;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\frac{b(z_2-z_1) E}{(1-\nu^{2})}  

    ;  D11                                                    b(z23z13)E3(1ν2)     D_{11}\;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\frac{b(z^3_2-z^3_1) E}{3(1-\nu^{2})}      

     A11                                                    b(z2z1)E(1ν2)  \displaystyle   A_{11}\;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\frac{b(z_2-z_1) E}{(1-\nu^{2})}  
      A55t                                                    bG1(z2z1)  \displaystyle     A^t_{55}\;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;{bG_1(z_2-z_1) }   

    N.B. N. B.  The equations are identical for top and bottom beams, so, the equations have been written either in a general form or only one of the two. Subscripts/identifiers may be used appropriately.​


    I would like to thank the National Science Academies for organising this wonderful Summer Resesarch Fellowship Program 2018 and giving me the golden opportunity to work under Prof. S. Gopalakrishnan, for whom I cannot find words to express my sincere gratitude and thanks. Despite being busy, he always managed time for me and I could learn a lot during my stay here.

    He has given me such an amazing learning experience under the caring guidance of his senior PhD scholar, Shweta Paunikar Madam. She was continuously helping me to learn something new, gain as much knowledge as I can in an enjoyable, yet deep way. Her valuable research time was given to me for my improvement and thanking her here is too less in fromt of her guidance.

    Other lab mates of the Professor's Computational Wave Mechanics Lab, especially Mr. Sanjay Sharma, Mr. Manish Raut, Mr. Jay, Mr. Atul, Dr. Shivashankar, Mr. Ishaq, Mr. Rajshekhar and Mr. Venkata S M have also helped me a lot along with a senior PhD scholar, Mr. Sudipta from Rotocraft and Mechanics lab.

    We were provided with a wonderful platform for writing well-designed reports, which look professional. I would like to thank the AuthorCafe team for their support and they were readily available to solve any technical difficulties.

    I cannot forget my parents and all others who direcly/indirectly helped me in this project, including my co-interns, whose questions taught me how to explain this subject to a person with varying levels of understanding, starting from novice, to someone who is deep into the Mathematics, and who also taught me various other topics of their studies. In this I would also like to thank the faculty interns, including Dr Biren Sir and others, whose motivation were crucial, who gave me ideas which matched with my subject of project as well. There may be many others, including those mentioned above without names, whose names may not have been written, but I would like to acknowledge their efforts towards the upliftment of my knowledge, which led to an increase in the level of the study.


    • Srinivasan Gopalakrishnan, (2016), Wave Propagation in Materials and Structures

    • Robert M. Jones, (1999), Mechanics of Composite Materials, 2nd Indian Reprint  -Textbook

    • J. N. Reddy, ( 2005 , )An Introduction to the Finite Element Method -Textbook

    • Klaus-Jürgen Bathe, (1996),  Finite Element Procedures- Textbook

    • Robert D. Cook,David S. Malkus, Michael E. Plesha, Robert J. WItt, (2007) Concepts and Applications of FInite Element Analysis- Textbook

    • James F. Doyle,  (1997), Wave Propagation in Structure - Textbook.

    • V. Rajaraman, (1997), Computer Programming in FORTRAN 77- Textbook,


    • Fig 1: Idea: Shweta Mam's explanation; but built in:MS Paint 
    Written, reviewed, revised, proofed and published with