Wave Motion in an Elastically Connected Double Beam using Superconvergent Finite Element Formulation
Guided by:
Abstract
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 beamrod elements, capable of deforming in axial (rodlike), transverse (Euler Bernoulli beamlike) 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.
Abbreviations
FEM  Finite Element Method 
FRP  Fiber Reinforced Polymer 
ODE  Ordinary Differential Equation 
SFEM  Spectral Finite Element Method 
FEA  Finite Element Analysis 
u  axial displacement in beam 
`w  transverse displacement 
φ  rotation of plane normal u (axial) and w(transvers)displacements of the beam 
INTRODUCTION
Background/Rationale
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 yieldingintension 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 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 EulerBernoulli beam, say has a differential equation of the form:
This in general can be solved if EI(x) is not a nonlinear function in x; but if say it is nonlinear, 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 offdiagonal 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 faroff node so elements is very less (corresponding to those in the corners of the opright and bottomleft regions of the matrix).
 Nonnegative 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 rightward 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 crosssection leads to the necessity of a shear correction factor to give results with better accuracy. In our case, we have chosen a rectangular crosssectioned beam which has a shear correction factor γ close to 5/6.
LITERATURE REVIEW
Information
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 KlausJü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 NewmarkBeta 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.
METHODOLOGY
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 ẍ = $\frac{d^{2^{}}x}{dt^2}$ and ÿ=. $\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 nonlinear.
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 beamrod elements, connected by continuous elastic springs of stiffness K=KF, are shown below:
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:
Equations in the new notations:
These equations as shown above, may be solved in the same way as the Timoshenko beam
Differentiating equations (5) and (7) and rearranging them, and then substituting them in equations(6) and (8) we get:
$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:
where a's denote the 3 constants.
From (4) and (6)
Differentiating (7) & using the above 3 equations, yields, on simplification :
Now, this may easily yield $\varphi^b$as below, from which $\varphi^t$found out .
By substitution in equations above, 'w', the component of transverse displacements, normal to the longitudinal axis of the beam may be found to be :
This can be substituted back to the original governing equations, to get the axial components of the displacements u for each beam to be :
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 Superconvergent 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{\}^{e_{}}}_{12\;x\;1}\;=\;\lbrack A\rbrack_{\;12\;x\;12}\;\{cons\tan t\}_{12\;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} =[ $u_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 :
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}\;=\;\lbrack D\rbrack _{6\;x\;12}\; \{const\}$
So, $\{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 $\lbrack N\rbrack\;=\;\lbrack D\rbrack\lbrack A\rbrack^{1}$
The strain components $[\varepsilon]\;=\;\lbrack \partial] \{U\}$
where the [∂] is given by [∂] = $\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, $\lbrack\varepsilon\rbrack\;=\;\lbrack\partial\rbrack\lbrack N\rbrack\{u\}^e\;\; = [B] \{u\}^e\;$ ,
where [B] = [∂D] $\lbrack A\rbrack ^{1}$ and
$\lbrack\partial D\rbrack=\;\lbrack\partial\rbrack\;\lbrack D\rbrack$
Now, stiffness mat.
[K] = $\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 $[E_1\; G_1\; E_2\; G_2]$ and similarly, we have mass matrix
$\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 nonsingular. 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 nonnegative real parts in the diagonal for both the stiffness and mass matrices.
Time marching using NewmarkBeta 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 NewmarkBeta 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:
c This is to solve system of eqns in NN unknowns
c
DIMENSION A(4,5)
DIMENSION X(4)
READ *,NN
READ *,((A(I,J),J=1,NN+1),I=1,NN)
DO 2 M=2,Nn
DO 1 N=1,M1
IF(A(M,N).NE.0.)THEN
IF(A(M1,N).NE.0) THEN
P = A(M,N)/A(M1,N)
L=1
ELSE
P = A(M,N)/A(M2,N)
L=2
ENDIF
pRint *, 'row', ml
DO 3 J=1,NN+1
A(M,J) = A(M,J) + P*a(ML,J)
3 CONTINUE
ENDIF
1 CONTINUE
pRint *,'U.T.',((A(I,J),J=1,NN+1),' \n ',I=1,NN)
2 CONTINUE
pRint *,'U.T.',((A(I,J),J=1,NN+1),'\n ',I=1,NN)
C BACK PUTTING/SUBSTITUTION
DO 4 I=NN,1,1
Q = 0.
DO 5 J = NN,I+1,1
Q = Q+ X(J)*a(I,J)
5 CONTINUE
X(I) = (A(I,NN+1)  Q) / A(I,I)
4 CONTINUE
DO 6 I=1,NN
PRINT *, 'X(',I,') = ',X(I)
6 CONTINUE
STOP
END
The code solves a system of linear algebraic equations using :
 GaussElimination 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 NewmarkBeta algorithm using MATLAB, which would no doubt, make the code look simpler and elegant but at the cost of using inbuilt functions in MATLAB, like a simple command A\b is equivalent to a dozens of FORTRAN lines.
Conclusion and Recommendations
Conclusion
The stiffness and mass matrices obtained matched with the theoretical reqirement that it should be symmetric with the diagonal terms having nonnegative 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.
APPENDICES
Appendix A
p_t ${D}_{11}^{t}\frac{{{B}_{11}^{t}}^{2}}{{A}_{11}^{t}}$
p_b ${D}_{11}^{b}\frac{{{B}_{11}^{b}}^{2}}{{A}_{11}^{b}}$
s_t $\frac{p\_t}{\gamma {A^t}_{55}}$
Q_t $\frac{KF}{\gamma {A^t}_{55}}$
O_t KF/p_t
O_b KF/p_b
$\beta_1$ Q_t  Q_b
$\beta_2$ O_T  O_b
$\beta_3$ O_b
$\beta_4$ O_b * s_t
$B_{11}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{b(z_2^2z_1^2)E}{2(1\nu^2)}\;\;$
; $D_{11}\;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\;\frac{b(z^3_2z^3_1) E}{3(1\nu^{2})}$
$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.
ACKNOWLEDGEMENTS
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 welldesigned 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 cointerns, 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.
References

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

KlausJü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,
Source

Fig 1: Idea: Shweta Mam's explanation; but built in:MS Paint
Post your comments
Please try again.