\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small Ninth MSU-UAB Conference on Differential Equations and Computational Simulations. \emph{Electronic Journal of Differential Equations}, Conference 20 (2013), pp. 133--149.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \setcounter{page}{133} \title[\hfilneg EJDE-2013/Conf/20/ \hfil Computational study] {Computational study of a dynamic contact problem} \author[J. Patel, J. Turi \hfil EJDE-2013/conf/20 \hfilneg] {Jigarkumar Patel, Janos Turi} % in alphabetical order \address{Jigarkumar Patel \newline Department of Mathematical Sciences\\ University of Texas at Dallas\\ Richardson, TX 75080, USA} \email{jsp061000@utdallas.edu} \address{Janos Turi \newline Department of Mathematical Sciences\\ University of Texas at Dallas\\ Richardson, TX 75080, USA} \email{turi@utdallas.edu} \thanks{Published October 31, 2013.} \subjclass[2000]{34N05, 65P99, 49M99, 90-08, 90C05, 37M05, 35L87, \hfill\break\indent 93C20, 93C55} \keywords{Vibrations of cracked beams; elastic structures with defects; \hfill\break\indent linear complementarity problem with FEM} \begin{abstract} In this article, we describe a computational framework to study the influence of a normal crack on the dynamics of a cantilever beam; i.e., changes in its natural frequency, amplitude and period of vibration, etc. Due to the opening and closing of the crack during beam vibrations, unilateral contact boundary conditions are assumed at the crack location. In the numerical implementation the contact conditions lead to the consideration of a linear complementarity problem. An effective solution strategy for this problem using a modification of the simplex method is presented. Numerical experiments are included. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{definition}[theorem]{Definition} \allowdisplaybreaks \section{Introduction} It is a standard assumption that new structures are ideal, in other words they do not contain any defects. As per their usage over time, imperfections tend to develop at various locations. Cracks represent a severe form of defects in elastic structures and their effects on system behaviour has aroused considerable interest in the last few decades. As a particular example, the dynamic response of a beam with a normal surface crack has been studied by many researchers (e.g. see \cite{dim96,sfe01,CD98,PD89}). In reality, the most likely regions of structures that contain cracks are joints and corners, but the above papers did not address this situation. In this article, we use a variational framework (see also \cite{KO88}) to consider the dynamics of a two dimensional elastic cantilever beam with a normal surface crack located at the supporting wall (see also \cite{dt95}). The wall is considered to be a rigid object and non penetration unilateral contact conditions (e.g. see \cite{bs89,CJM05,glo99,ls84}) are assumed at the crack location during beam vibrations. Using classical variational analysis the problem can be reduced to an associated hyperbolic two point boundary-value problem with contact boundary conditions. For numerical simulations, we discretize the variational inequality and the unilateral contact conditions and obtain an associated linear complementarity problem \cite{CPS92}. Our primary goal is to introduce an effective method for finding the unique optimum (minimum) solution of this problem. Numerical results comparing vibrations of a cracked beam and an ideal beam indicate changes in frequency, amplitude and modes of vibration (i.e. Natural frequency and modes of vibration are considerably effected by the presence of a normal crack.) The rest of the article is organized as follows: In Section $1$ we formulate the problem and introduce space and time discretizations for it. Formulation of boundary conditions and admissibility of the numerical solution of the variational inequality satisfied by a cantilever beam is defined in Section $2$. In Section $3$ we describe an effective way to solve the linear complementarity problem. Some numerical results comparing the frequency of an ideal beam vibrations and a cracked beam vibrations are discussed in Section $4$ prior to the concluding comments in Section $5$. \subsection{Problem formulation} Consider a two dimensional cantilever beam with a normal surface crack at the supporting wall (see Figure \ref{fig:f1}.) We assume that only normal stresses are supported at the crack location in case of contact. \begin{figure}[ht] \begin{center} \includegraphics[width=0.6\textwidth]{fig1} % fig2 \end{center} \caption{Two dimensional cantilever beam with crack at the supporting wall} \label{fig:f1} \end{figure} Let $x=(x_1,x_2)$ be a two dimensional space variable. The quantities $l$ and $\theta$ represent length and thickness of the beam. The functions $u(t,x)$, $\varepsilon(t,x)$, and $\sigma(t,x)$ denote the displacement, linearized strain, and linearized stress, respectively at time $t$ of the space element with Lagrangian coordinates $x$. Also $\Omega = [0,l] \times [0,h]$ is the domain occupied by a cantilever beam in reference configuration. Let $\Gamma$ be the boundary of a cantilever beam. $\Gamma_{c}$ denotes the boundary of the normal crack. The elasticity modulus of the beam, $g(x)$ is a fourth order tensor which is symmetric and positive definite. All physical parameters are assumed to be independent of time and continuously differentiable with respect to $x$. Given some body force $f(x)$ the resulting deformation field $u(x)$ minimizes the potential energy in the static problem. In this case the potential energy of the system is given by \begin{equation} \label{e1.1} E_{\rm pot}=\frac {1}{2} \int_\Omega \frac{\partial u_i}{\partial x_{j}}(x) g_{ijkl}(x) \frac{\partial u_{k}}{\partial x_{l}}(x) dx - \int_\Omega f_i(x)u_i(x)dx. \end{equation} Here, $W^{1,2}(\Omega, R^2)$ is the natural state space for the deformation field. The displacement, $u$, is restricted to the closed convex set (see also \cite{dt98}) $\tilde{J}$ by Dirichlet type boundary conditions. $\tilde{J}$ is called admissibility set and it is given by \begin{equation} \label{e1.2} \begin{split} \tilde{J} = \big\{& u \in W^{1,2}(\Omega, R^2) : u(x,y) = 0, \text{ for $x = 0$ and $y \in [0,d]$} \\ & u(x,y) \geq 0 \text{ for } (x,y) \in \Omega\backslash\{ x=0,\, y \in [0,d] \} \big\} \end{split} \end{equation} The set $\tilde{J}$ contains the elements which are connected to or are in contact with the wall. The phase of the beam to the right of a crack cannot penetrate the rigid wall during vibration. Using Hamilton's principle of least action, the displacement field $u$ of the cantilever beam must satisfy the minimization problem: \begin{equation} \label{e1.3} \text{Find $u \in \tilde{J}$ such that } E_{\rm pot}(u) \leq E_{\rm pot}(v),\; \forall v \in \tilde{J} \end{equation} or equivalently satisfy the variational inequality \begin{equation} \label{e1.4} \frac {1}{2} \int_\Omega \frac{\partial u_i}{\partial x_{j}}(x) g_{ijkl}(x) \Big( \frac{\partial v}{\partial x_{l}}(x) - \frac{\partial u_{k}}{\partial x_{l}}(x)\Big) dx - \int_\Omega f_i(x)(v_i(x) - u_i(x))dx \geq 0. \end{equation} In the set up for the dynamic problem the displacement, $u$ depends also on time $t$. In this case the static body force $f$ is replaced by the inertial forces. \begin{equation} \label{e1.5} f(t,x) = - \rho \frac{\partial^{2} u}{\partial t^{2}}(t,x) \end{equation} Find $u(t,.) \in \tilde{J}$ for all $t \geq 0$ and all $v \in \tilde{J}$ satisfying the following variational inequality. \begin{equation} \label{e1.6} \begin{aligned} &\int_\Omega \frac{\partial u_i}{\partial x_{j}}(t,x) g_{ijkl}(x) \Big( \frac{\partial v_{k}}{\partial {l}}(x) - \frac{\partial u_{k}}{\partial x_{l}}(t,x) \Big) dx \\ &+ \int_\Omega \frac{\partial^{2} u}{\partial t^{2}}(t,x)\rho(x)(v_i(x) -u_i(t,x))dx \geq 0. \end{aligned} \end{equation} satisfying the initial and boundary conditions \begin{equation} \label{e1.7} u_i(0,x) = u_{0i}(x), \quad \frac{\partial {u_i}}{\partial t}(0,x)=u_{1i}(x) \quad \text{on } \Gamma, \end{equation} Here $ u_{0i} \in \tilde{J}$ and $ \frac{\partial u_i}{\partial t} \in L^{2}(\Omega, R^{2})$ for all $i$. Note that throughout this section we use the word admissible to represent elements in the set $\tilde{J}$. Discretization of the space variables $x_1$ and $x_2$ convert the displacement field into the form \begin{equation} \label{e1.8} u(t,x_1,x_2) = \sum_{n=1}^{N} w^n(t,x_1)\varphi^n(x_2) \end{equation} Here $\varphi^1$,$\varphi^2$,\dots ,$\varphi^N$ are linearly independent shape functions and $w: [0,\infty) \times [0,l] \to R^N$ is $n \times 1$ vector valued function. We want to find an admissible function $ w: [0,\infty) \times [0,l] \to R^N$ which satisfies the initial and boundary conditions and such that for all $t \geq 0 $ and for all admissible $v$ the following variational inequality holds \begin{equation} \label{e1.9} \begin{aligned} &\int_0^h \int_{0}^{l} \frac{\partial}{\partial x_{j}} \Big( \sum_{m=1}^N w^m(t,x_1) \varphi_i^m(x_2) \Big) g_{ijkl}(x_1,x_2) \\ &\times \frac{\partial}{\partial x_{k}} \Big( \sum_{n=1}^N (v^n(x_1)-w^n(t,x_1))\varphi^n_k(x_2)\Big) dx_2 dx_1 \\ &+ \int_0^h \int_{0}^{l} \Big( \sum_{m=1}^N \frac{\partial^{2} u}{\partial t^{2}}(t,x_1) \varphi_i^m(x_2)\Big) \rho(x_1,x_2) \\ &\times \Big( \sum_{n=1}^N (v^n(x_1) -w^n(t,x_1))\varphi_i^n(x_2)\Big) dx_2 dx_1 \geq 0 \end{aligned} \end{equation} Straightforward but tedious calculations convert the above integral inequality into the form \begin{equation} \label{e1.10} \begin{aligned} &\int_0^h \Big[ \frac{\partial w^{T}}{\partial x_{1}}(t,x_1)R(x_1) \Big( \frac{\partial v}{\partial x_{1}}(x_1) - \frac{\partial w}{\partial x_{1}}(t,x_1)\Big) \\ & + w^T(t,x_1)P^T(x_1)\Big( \frac{\partial v}{\partial x_{1}}(x_1) -\frac{\partial w}{\partial x_{1}}(t,x_1)\Big) \\ & + \frac{\partial w^{T}}{\partial x_{1}}(t,x_1)P(x_1)(v(x_1)-w(t,x_1)) + w^T(t,x_1)Q(x_1)(v(x_1)-w(t,x_1)) \\ &+ \frac{\partial^{2} w^{T}}{\partial t^{2}}(t,x_1)M(x_1)(v(x_1)-w(t,x_1)) \Big] dx_1 \geq 0, \end{aligned} \end{equation} where $R$, $P$, $Q$ and $M$ are $N \times N$ matrices given as follows. Note that $R$ and $M$ are positive definite and $P$ and $Q$ are positive semi definite matrices, respectively. \begin{equation} \label{e1.11} \begin{gathered} M(x_1)=(m^{mn}(x_1)) = \int_{0}^{l} \Big( \sum_{m,n=1}^N \varphi^m_i(x_2) \rho(x_1,x_2) \varphi^n_i(x_2)\Big) dx_2 , \\ R(x_1)=(r^{mn}(x_1)) = \int_{0}^{l} \Big( \sum_{m,n=1}^N \varphi_i^m(x_2) g_{i1k1}(x_1,x_2) \varphi_k^n(x_2)\Big) dx_2, \\ P(x_1)=(p^{mn}(x_1)) = \int_{0}^{l} \Big(\sum_{m,n=1}^N \varphi_i^m(x_2) g_{i1 k\beta}(x_1,x_2) \frac{\partial \varphi_{k}^n}{\partial \beta}(x_2)\Big) dx_2,\\ Q(x_1)=(q^{mn}(x_1))= \int_{0}^{l} \Big( \sum_{m,n=1}^N \frac{\partial \varphi_i^m}{\partial \alpha} (x_2) g_{i\alpha k\beta}(x_1,x_2) \frac{\partial \varphi_{k}^n}{\partial \beta} (x_2)\Big) dx_2. \end{gathered} \end{equation} To convert \eqref{e1.10} into a partial differential equation, let $z \in D([0,l],R^n)$ be a test function and let $z = v(x_1) - w(t,x_1)$. Integration by parts with respect to $x_1$ yields \begin{equation} \label{e1.12} \begin{aligned} &\int_0^h \Big[ \frac{\partial w^{T}}{\partial x_{1}}R(t,x_1)(x_1) + w^T(t,x_1)P^T(x_1) \\ &+ \frac{\partial w^{T}}{\partial x_{1}}(t,x_1)P(x_1) + w^T(t,x_1)Q(x_1) + \frac{\partial^{2} w^{T}}{\partial t^{2}}(t,x_1)M(x_1)]z \Big] dx_1 \geq 0 \end{aligned} \end{equation} Since the above equality holds for any test function so it holds for $-z$ also. We obtain it in the sense of distributions \begin{equation} \label{e1.13} M(x_1)\frac{\partial^{2} w^{T}}{\partial t^{2}}(t,x_1) = \frac{\partial}{\partial x_1}(R(x_1)\frac{\partial w^{T}}{\partial x_{1}} + P(x_1)w) - P^{T}(x_1)\frac{\partial w}{\partial x_{1}}(t,x_1) - Q(x_1)w(t,x_1) \end{equation} Now for an arbitrary $v \in \tilde{J}$ integrate again by parts, we have only boundary terms, \begin{equation} \label{e1.14} \Big(\frac{\partial w}{\partial x_1}(t,x_1)R(x_1) + w^{T}(t,x_1)P^{T}(x_1)\Big) (v(x)-w(x_1))\mid_{0}^{h} \geq 0 \end{equation} At the free end $v(l)- w(l)$ can take any value in $R^{N}$. \begin{equation} \label{e1.15} \big(\frac{\partial w}{\partial x_1}(l)R(l) + w(l)P^{T}(l)\big) = 0. \end{equation} At contact, $v(0)$ can take any value in $\tilde{J}$. \noindent\textbf{Problem:} Find admissible function $w: [0,\infty) \times [0,l] \to R^N$ for all $t \geq 0$ and for all admissible $v$ such that \begin{equation} \label{e1.16} \begin{gathered} \begin{aligned} M(x_1)\frac{\partial^{2} w^{T}}{\partial t^{2}}(t,x_1) &= \frac{\partial}{\partial x_1}(R(x_1)\frac{\partial w^{T}}{\partial x_{1}} + P(x_1)w) - P^{T}(x_1)\frac{\partial w}{\partial x_{1}}(t,x_1) \\ &\quad - Q(x_1)w(t,x_1) \end{aligned} \\ w_i(0,x) = w_{0i}(x), \quad \frac{\partial {w_i}}{\partial t}(0,x)=w_{1i}(x) \quad\text{on } \Gamma \\ v(l) - w(l) \quad\text{is free}\\ v(0) \in \tilde{J} \end{gathered} \end{equation} The first line in \eqref{e1.16} is the partial differential equation satisfied by cantilever beam while rest of the lines specify the initial and boundary conditions for cantilever beam with a crack at the supporting wall. Note that the partial differential equation in \eqref{e1.16} can be written in the equivalent form. \begin{equation} \label{e1.17} M \frac{\partial^{2} w}{\partial t^{2}}(t,x_1) = R \frac{\partial^{2} w}{\partial x_{1}^{2}} + [ \frac{\partial R}{\partial x_{1}}] \frac{\partial w}{\partial x_{1}} - [ \frac{\partial P}{\partial x_{1}} - Q ]w \end{equation} \subsection{Time discretization} In this section, we discretize the variational inequality with respect to the time variable. Replacing $ \frac{\partial^{2} w}{\partial t^{2}}$ using second order symmetric finite differences in \eqref{e1.17} \begin{equation} \label{e1.18} \frac{\partial^{2} w}{\partial t^{2}} = \frac{w(t+k) - 2w(t) + w(t-k)}{k^2} \end{equation} leads to the equation \begin{equation} \label{e1.19} M\big( w(t+k, x_1) \big) = k^2 R \frac{\partial^{2} w}{\partial x_{1}^{2}} + k^2 [ \frac{\partial R}{\partial x_{1}} ] \frac{\partial w}{\partial x_{1}} + [ k^2 \frac{\partial P}{\partial x_{1}} - k^2 Q + 2M ]w + M w(t-k, x_1). \end{equation} Here $w(t-k, x_1)$ represent the previous location of each particle of the beam. Similarly,$w = w(t, x_1)$ and $w(t+k,x_1)$ represent current and future locations of each particle, respectively. Using the initial and boundary conditions, the location of each particle can be calculated using \eqref{e1.19}. Since expression on the right side in \eqref{e1.19} can be calculated easily, we use following notation for further discussion. \begin{equation} \label{e1.20} - b = k^2 R \frac{\partial^{2} w}{\partial x_{1}^{2}} + k^2 [ \frac{\partial R}{\partial x_{1}}] \frac{\partial w}{\partial x_{1}} + [ k^2 \frac{\partial P}{\partial x_{1}} - k^2 Q + 2M ]w + R w(t-k, x_1) \end{equation} Equation \eqref{e1.19} takes the form \begin{equation} \label{e1.21} M\big( w(t+k, x_1) \big) + b = 0 \end{equation} \section{Numerical implementation} In this section we model the contact boundary condition for the cantilever beam with normal crack located right at the supporting wall. Body forces influence the dynamics of the beam which can be explained by studying the effect of the point load at the free end of the beam. \subsection{Modeling the boundary conditions} In this section, we model the boundary conditions for two dimensional discretized ideal and cracked beams (e.g. see \cite{SSB03} and \cite{KB96}). For an ideal beam, we simply solve linear system represented in \eqref{e1.21} and obtain the displacement field at a particular time step. Since $M$ is positive definite matrix this system admits a unique solution. In fact, the solution of system represented by \eqref{e1.21} for ideal beam is given by, \begin{equation} \label{e2.1} w(t+k, x_1) = M^{-1}(-b) \end{equation} Hence finding the displacement field for the dynamic problem corresponding to an ideal beam is rather straightforward. On the other hand, dynamic problem corresponding to a beam with crack at the supporting wall is challenging because of an additional boundary condition that is satisfied at the crack location during vibrations. \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig2a} %fg3a.eps \includegraphics[width=0.45\textwidth]{fig2b} %fg3b.eps \end{center} \caption{(a) Discretized ideal beam. (b) Discretized cracked beam} \label{fig:f2} \end{figure} To model the boundary condition at the crack location consider Figure \ref{fig:f2}(a). In the reference figure, heights of the beams are discretized into $N_{t}$ vertical nodes and the lengths of beams are discretized into $N_{l}$ horizontal nodes. The discretized ideal beam can be viewed in terms of two disjoint node sets $J_1$ and $J_2$ as shown in the reference Figure \ref{fig:f2}. We discuss now the modeling of the contact boundary conditions for two dimensional beams with a normal crack located right at the supporting wall. The most convenient way to apply contact boundary conditions is to divide the whole discretized beam into several point sets according to the crack location. Figure \ref{fig:f2}$(b)$ represents different point sets $J_1$, $J_2$ and $J_3$ of the two dimensional beam when the normal crack is located right at the wall. In reference Figure \ref{fig:f2}(b) the length of the crack is $N_{f}$. \begin{equation} \label{e2.2} \begin{gathered} J_1=\{(i,j)/ i=0, j \leq N_t \}, \\ J_2=\{(i,j)/ i>0, j > 0 \}, \end{gathered} \end{equation} Admissibility condition for the deformation field for an ideal beam is presented in form of the following definition. \begin{definition} \label{def1} \rm $w(i,j)$ is admissible if $$ w(i,j) = \begin{cases} 0 & \text{for } (i,j) \in J_1\\ \text{free} & \text{for } (i,j) \in J_2. \end{cases} $$ Here $w(i,j)$ is $free$ means that the displacement of the $(i,j)^{th}$ node can take any value corresponding to the applied body force. In our case, force is provided in such a way that we get observable vibration. \end{definition} To define a similar type of admissibility condition for the deformation field of a beam with a crack at the supporting wall, we divide the beam region into three disjoint node sets as shown in the reference Figure \ref{fig:f2}(b). \begin{equation} \label{e2.3} \begin{gathered} J_1=\{(i,j)/ i=0, j \leq N_x \}, \\ J_2=\{(i,j)/ i=0, j > N_x \}, \\ J_3=\{(i,j)/ (i,j)\notin \cup_{i=1}^{2} J_i\} \end{gathered} \end{equation} Set $J_1$ includes all those nodes which are connected to the wall and set $J_2$ includes all those nodes which are in contact with the wall. Set $J_3$ contains all other nodes of a beam in the reference configuration. Definition of admissibility for the deformation field of a discretized beam in case of a normal crack located at the supporting wall can be modeled as follows. \begin{definition}\label{def2.2}\rm $w(i,j)$ is admissible if $$ w(i,j) = \begin{cases} 0 & \text{for } (i,j) \in J_1\\ \geq 0 & \text{for } (i,j) \in J_2\\ \text{free} & \text{for } (i,j) \in J_3. \end{cases} $$ \end{definition} According to the above definition every node of $J_1$ should have displacement zero at any time step. Because of the non penetration condition, the phase of the beam cannot penetrate the wall. In other words, the displacement field for every node of $J_2$ is either positive or zero at any time step. If every node in $J_2$ has positive displacement at a particular time step then it means that at that time step the crack is completely open. If some nodes of $J_2$ have zero displacement and others have positive displacements, then the crack is partially open in the respective time step. If all nodes of $J_2$ have zero displacement at a particular time step, then at that time step the crack is completely closed. Nodes in $J_3$ will have displacement fields according to the applied force as well as the boundary conditions which lead to crack openings and closings. \subsubsection{Contact boundary conditions at the crack} Our problem is to find the displacement field, $w$, of a beam with a crack located at the supporting wall. The displacement field in this case depends not only on the applied force but on the boundary condition as well. Thus we would like to find $w$ such that \begin{equation} \label{e2.4} [Mw+b] = \begin{cases} 0 & \text{for } (i,j) \in J_1\\ \geq 0 & \text{for } (i,j) \in J_2\\ \text{free} & \text{for } (i,j) \in J_3. \end{cases} \end{equation} is satisfied. \begin{table}[ht] \caption{Table representing Equations \eqref{e2.4}} \begin{center} \begin{small} \begin{tabular}{|c|c|c|c|} \hline & $J_{2}$ ($w\geq 0$) $\quad$ $J_{3}$ ($w$ is free vector) & $-p$ vector & $-b$ = vector \\ \hline $\begin{bmatrix} J_{2} \\J_{3} \end{bmatrix}$ & $M = \begin{bmatrix} m_{11} & m_{12} & . & . & . & m_{1n} \\ m_{21} & m_{22} & . & . & . &m_{2n}\\. & . & .& . & . & .\\ . & . & .& . & . & .\\. & . & .& . & . & .\\ . & . & .& . & . & .\\. & . & .& . & . & . \\ m_{n1} & m_{n2} & . & . & . & m_{nn} \end{bmatrix}$ & $\begin{bmatrix} -1 \\.\\.\\-1\\0\\.\\.\\0 \end{bmatrix}$ & $\begin{bmatrix} -b_1 \\-b_2\\.\\.\\.\\.\\.\\-b_n \end{bmatrix}$ \\ \hline \end{tabular} \end{small} \end{center} \label{table1} \end{table} The displacement field of the entire beam with crack can be understood by studying the displacement of the crack. At the crack location unilateral contact conditions are satisfied whose discretization give rise to a linear complementarity problem. Let $p \geq 0$ be positive vector such that, \begin{equation} \label{e2.5} [Mw+b] = p \geq 0 \end{equation} Then the system represented by \eqref{e2.4} can be viewed in terms of matrices and vectors as shown in Table $(1)$. Here $p = [1,1,\dots ,1,0,\dots .0]^T \geq 0$ means that the initial condition is chosen in such a way that the crack is open. The displacement field for the system given in Table $(1)$ can be obtained by solving the given linear system. Being a positive definite matrix, $M$ is invertible and then there exists a displacement field $w$ for each time step. \begin{table}[ht] \caption{Displacement of every node for the static problem} \begin{center} \begin{small} \begin{tabular}{|p{0.5cm}|c|c|c|} \hline & $J_{2}$ $J_{3}$ & $-p$ & $-b$ \\ \hline $\begin{bmatrix} J_{2} \\ J_{3} \end{bmatrix}$ & $I$ & $ \begin{bmatrix} m_{11} & . & . & . & m_{1n} \\m_{21} & . & . & . &m_{2n}\\ . & .& . & . & .\\. & .& . & . & .\\. & .& . & . & .\\ . & .& . & . & .\\ . & .& . & . & . \\ m_{n1} & . & . & . & m_{nn} \end{bmatrix}^{-1}$ $\begin{bmatrix} -1 \\.\\.\\-1\\0\\.\\.\\0 \end{bmatrix}$ & $ \begin{bmatrix} m_{11} & . & . & . & m_{1n} \\ m_{21} & . & . & . &m_{2n}\\. & .& . & . & .\\. & .& . & . & .\\ . & .& . & . & .\\ . & .& . & . & .\\ . & .& . & . & . \\ m_{n1} & . & . & . & m_{nn} \end{bmatrix}^{-1}$ $\begin{bmatrix} -b_1 \\-b_2\\.\\.\\.\\.\\.\\-b_n \end{bmatrix}$ \\ \hline \end{tabular} \end{small} \end{center} \label{table2} \end{table} So the static problem related to this issue is easy to solve since it only involves the solution of a system of linear equations. Finding the displacement field for a dynamic problem requires a little more effort due to the fact that during vibrations the crack opens and closes. In such a case to get the displacement field for the beam one has to solve a related linear complementarity problem satisfied at the crack location. To find the deformation field of the dynamic problem, orthogonality should be assumed naturally to get an optimum solution. Also the natural state space for the deformation field is the Sobolev space $W^{1,2}(\Omega, R^{2})$. Consequently, at the crack location, the deformation field should satisfy \begin{equation} \label{e2.6} \begin{gathered} b+Mw = p\\ w \geq 0, \quad p \geq 0\\ w \perp p. \end{gathered} \end{equation} This equation is called linear complementarity problem. \section{A method for solving the linear complementarity problem} In this section, we introduce a very effective method to find the unique minimum solution of the associated linear complementarity problem (see \cite{CPS92,DL76}). We want to find $w$ and $p$ such that \eqref{e2.6} is satisfied. \begin{equation} \label{e3.1} \begin{gathered} b+Mw = p\\ w \geq 0, \quad p \geq 0\\ w \perp p. \end{gathered} \end{equation} The symmetry of $M$ allows us to convert it into the matrix where all diagonal entries are $1$ using the following similarity transformation: $ \tilde{M} = G^{-1}MG$ such that $\tilde{m_{ii}} = 1$ for $i = 1,2,\dots ,N$, where \begin{equation} \label{e3.2} G = \begin{pmatrix} \frac{1}{\sqrt m_{11}}&0&.&.&.&0\\ 0&\frac{1}{\sqrt m_{22}}&.&.&.&0\\ .&.&.&.&.&.\\ .&.&.&.&.&.\\ 0&0&.&.&.&\frac{1}{\sqrt m_{NN}} \end{pmatrix} \end{equation} Note that $G$ is invertible, symmetric and positive definite. \begin{equation} \label{e3.3} \tilde{M} = G^{-1}MG \Longleftrightarrow M = G\tilde{M}G^{-1} \end{equation} Using \eqref{e3.3}$ in \eqref{e2.6}$ we obtain \begin{equation} \label{e3.4} \begin{gathered} b+G\tilde{M}G^{-1}w = p\\ w \geq 0, \quad p \geq 0\\ w \perp p. \end{gathered} \end{equation} Multiplying \eqref{e3.4} by $G^{-1}$ converts the problem into the following form \begin{equation} \label{e3.5} \begin{gathered} G^{-1}b+G^{-1}GMG^{-1}w = G^{-1}p\\ G^{-1}w \geq 0, \quad G^{-1}p \geq 0\\ G^{-1}w \perp G^{-1}p. \end{gathered} \end{equation} We adopt following notation, $ \tilde{w} = G^{-1}w$, $ \tilde{b} = G^{-1}b$ and $\tilde{p} = G^{-1}p$ and then \eqref{e3.5} takes the form \begin{equation} \label{e3.6} \begin{gathered} \tilde{M}\tilde{w}+\tilde{b} = \tilde{p}\\ \tilde{w} \geq 0, \quad \tilde{p} \geq 0\\ \tilde{w} \perp \tilde{p}. \end{gathered} \end{equation} We use the simplex method to solve the associated linear complementarity problem. If $\tilde{p} = (\tilde{p}_{1},\tilde{p}_{2},\dots ,\tilde{p}_{n})^{T}$ has all elements strictly positive during every time step during the vibration of the beam, then due to the orthogonality condition every element of the displacement vector $w$ must be zero during each time step of the vibration. In such a case, we conclude that the beam is ideal. This observation may be used as an identification technique of possible cracks in the structure. If every element of $w$ is strictly positive during every time step during beam vibrations then due to the orthogonality condition every element of the vector $p$ must be zero during every time step of the vibration. This indicates that the crack stays open all the time during the vibration, which is not a realistic case. Hence some elements in both vector should be non negative and zero. Existence of a crack guarantees that there exist some elements $\tilde{w}_i$ of $\tilde{w}$ and $\tilde{p}_i$ of $\tilde{p}$ which are strictly positive depending on if the crack is open, closed or partially open. The components, $\tilde{b}_i$, of the vector $\tilde{b}$ depend on the derivatives of the displacement and on the history of the previous location according to \eqref{e1.20}. These components could be positive or negative. Again existence of the crack guarantees some negative components of the vector $\tilde{b}$. We replace those negative elements $\tilde{b}_i$ by $-\tilde{b}_i = \tilde{j}_i$ for $i=1,2,\dots ,n$. Thus $\tilde{j}$ is an $n \times 1$ vector whose components are all positive and greater than or equal to the corresponding component of $\tilde{b}$. Now our problem is to find $\tilde{w}$ and $\tilde{p}$ satisfying \begin{equation} \label{e3.7} \begin{gathered} \tilde{M}\tilde{w}-\tilde{p} \leq \tilde{j}\\ \tilde{w} \geq 0, \quad \tilde{p} \geq 0\\ \tilde{w} \perp \tilde{p}. \end{gathered} \end{equation} Introducing a slack variable to the above linear inequality we obtain \begin{equation} \label{e3.8} \begin{gathered} \tilde{M}\tilde{w}-\tilde{p}+\tilde{k} = \tilde{j}\\ \tilde{w} \geq 0, \quad \tilde{p} \geq 0\\ \tilde{w} \perp \tilde{p}. \end{gathered} \end{equation} Here $\tilde{k}$ is an $n \times 1 $ vector and its $i^{th}$ entry is non zero if $\tilde{b}_i \leq 0$. In order to formulate a linear programming problem, we need an objective function. The coefficient vector in the objective function is an $n \times 1$ vector, say $\tilde{l}$ whose $i^{th}$ component is $1$ if the corresponding $i^{th}$ slack variable is non zero. If $i^{th}$ slack variable is zero then the corresponding component of $\tilde{l}_i$ is $0$. Using this information, coefficients of the objective row can be generated in terms of the vector given by $ \tilde{O}_{R} = -\tilde{l}\tilde{M}$. Using the simplex method, we find the solution of the linear programming problem which minimizes the objective function. Moreover the solution of the linear complementarity problem satisfies the non negativity and orthogonality conditions. Problem: \begin{equation} \label{e3.9} \begin{gathered} \text{Minimize } z = \tilde{O}_{R}x \quad \text{where } x = [x_1\;x_2\;\dots\; x_n]^{T}\\ \text{subject to constraints: } \tilde{M}\tilde{w}-\tilde{p}+\tilde{k} = \tilde{j} \\ \tilde{w} \geq 0, \quad \tilde{p} \geq 0,\quad \tilde{w} \perp \tilde{p}. \end{gathered} \end{equation} During each pivoting step, we have to keep track of basic and non basic artificial variables. Due to this reason, we define the column and row indices $C_{I} = [-1,-2,\dots ,-n]$ and $R_{I} = [1,2,\dots ,n]$, respectively. To apply a pivoting step, we arrange all vectors and matrices in the table $T$ as follows \begin{equation} \label{e3.10} T =\begin{pmatrix} 0 & C_{I} & 0\\ (R_{I})^T & \tilde{M} & \tilde{j}\\ 0 & \tilde{O_R} & 0\\ 0 & \tilde{l} & 0\\ \end{pmatrix} \end{equation} We find the pivot column by finding $\min \big(\tilde{l}_i/\tilde{o}_i\big)$ for $i = 1,2,\dots ,n$. Similarly, we find the pivot row by finding $\min\big(\tilde{j}_i/a_{ik} \big)$ where the $a_{ik}$'s are elements of $k^{th}$ pivot column. The common element of the pivot column and pivot row is called the pivot element. Using the row reduction operation, we make the pivot element 1 and the rest of the entries in the pivot column equal to $0$. To keep track of the basic artificial variables, we interchange respective row and column indices. Repeat the pivot step until all artificial variables becomes non-basic. After all pivoting steps are completed $T$ takes the form \begin{scriptsize} \begin{equation*} %\label{e} T = \begin{pmatrix} 0&(-1)\text{ or } 1&(-2)\text{ or }2&(-3)\text{ or }3&(-4)\text{ or }4&.&.&.&(-n)\text{ or }n&0\\ 1\text{ or }(-1)&.&.&.&.&.&.&.&.&\tilde{j}_{c1}\\ 2\text{ or }(-2)&.&.&.&.&.&.&.&.&\tilde{j}_{c2}\\ 3\text{ or }(-3)&.&.&.&.&.&.&.&.&\tilde{j}_{c3}\\ 4\text{ or }(-4)&.&.&.&.&.&.&.&.&\tilde{j}_{c4}\\ .&.&.&.&.&.&.&.&.&.\\ .&.&.&.&.&.&.&.&.&.\\ .&.&.&.&.&.&.&.&.&.\\ n\text{ or }(-n)&0&0&0&0&.&.&.&.&\tilde{j}_{cn}\\ 0&.&.&.&.&.&.&.&0&0\\ 0&0&0&0&0&.&.&.&0&0\\ \end{pmatrix} \end{equation*} \end{scriptsize} Here $[\tilde{j}_{ci}]$ for $i = 1,2,\dots ,n$ are the new elements obtained from the $[\tilde{j}_i]$s after all pivoting steps are completed. From the values of $[\tilde{j}_{ci}]$, we generate $\tilde{w}$ and $\tilde{p}$. \begin{equation} \label{e3.12} \tilde{w}_i = \begin{cases} \tilde{j}_{ci} & \text{for } R_I < 0 \\ 0 & \text{for } R_I > 0. \end{cases} \end{equation} Elements of the vector $\tilde{p}$ are equal to $\tilde{p_i} = \tilde{j_{ci}} - \tilde{w_i}$. This is a clever way of preserving orthogonality and non negativity of the solution of the linear complementarity problem. Using the simplex method we obtained vectors $\tilde{w}$ and $\tilde{p}$ satisfying the linear programming problem given by \eqref{e3.5}. The vectors $w$, $b$ and $p$ can be retrieved from $\tilde{w}$, $\tilde{b}$ and $\tilde{p}$ using the matrix $G$ as follows. \begin{equation} \label{e3.13} \begin{gathered} w = G\tilde{w}\\ b = G\tilde{b}\\ p = G\tilde{p}. \end{gathered} \end{equation} The simplex method is used effectively to solve the linear complementarity problem satisfied at the crack location. We state the following uniqueness theorem about the solutions of the linear complementarity problem. \begin{theorem}[Uniqueness Theorem] \label{thm3.3} Using the simplex method, the problem \begin{gather*} M_{n \times n}w_{n \times 1} + b_{n \times 1} = p_{n \times 1}\\ w_{n \times 1} \geq 0, \quad p_{n \times 1} \geq 0\\ w \perp p. \end{gather*} has a unique optimum (minimum) solution $w_{n\times 1}$ and $p_{n \times 1}$ where $M_{n \times n}$ is a symmetric positive definite matrix. \end{theorem} The proof of the this theorem follows immediately due to positive definiteness of matrix $M$, and the positivity and orthogonality of the solution. \section{Results} Using the finite element discretization, we compare the vibrations of a cantilever beam having the normal surface crack at the supporting wall with vibrations of an ideal beam. \emph{Timestep} vs \emph{amplitude} plots are presented for an ideal beam and for a cantilever beam with a normal crack located right at the supporting wall. The following parameters are used in the finite element discretization. \begin{itemize} \item Length of the beam $l = 12$ cm. \item Thickness of the beam $\theta$ = $2.2$ cm. \item Lame parameters $\lambda = 1$ and $\mu = 1$ \item Density $\rho = 0.1$ \item Damping factor $d = 0.002$. \item Initial force $1000$ unit \end{itemize} \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig3a} %fg4a.eps \includegraphics[width=0.45\textwidth]{fig3b} \\ %fg4b.eps \includegraphics[width=0.45\textwidth]{fig3c} %fg4c.eps \end{center} \caption{(a) Ideal beam vibration. (b) Beam vibration with crack length 4. (c) Beam vibration with crack length 6} \label{fig3} \end{figure} Figure \ref{fig3} (a), (b) and (c) represent positions of an ideal beam and cracked beam with $8$ vertical and $5$ horizontal nodes in finite element discretization. These results are taken by freezing each of these beams at various time steps. Length of the crack in Figure \ref{fig3}(b) is kept 4 while the crack length on Figure \ref{fig3}(c) is kept 6. The crack is completely open in Figure \ref{fig3}(b) and (c). The force applied in all of these cases is the same. \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig4} %{Merg.eps}}%<---angle here \end{center} \caption{Ideal beam vibration and cracked beam vibration.} \label{fig:Mr} \end{figure} Figure \ref{fig:Mr} represent the amplitudes of vibrations of an ideal and a cracked beam with $8$ vertical and $5$ horizontal nodes in finite element discretization. The amount of point load or force applied on the free end of the beam in both cases is the same. We run the simulation for $2000$ time steps for both cases. Length of the crack in Figure \ref{fig:Mr} is 4. In other words, four consecutive nodes on the supporting (fixed) end are kept loose. Dotted curve represents the vibration of a beam with crack size 4 and solid curve represents vibration of an ideal beam. Comparing both curves in Figure \ref{fig:Mr}, we notice that the amplitude in the negative direction increases and positive direction is decreasing due to presence of a crack at the supporting wall. Noticeable change in the period and amplitude can be observed in the reference Figure \ref{fig:Mr}. \begin{table}[ht] \begin{center} \begin{small} \begin{tabular}{|c|c|c|c|c|} \hline $1^{st}$ mode in $Y_-$ & $2^{nd}$ mode in $Y_-$ & $1^{st}$ mode in $Y_+$ & $2^{nd}$ mode in $Y_+$ & Period \\ \hline 0.0953 & 0.0789 & 0.0849 & 0.0744 & 1.85 \\ \hline \end{tabular} \end{small} \end{center} \caption{Amplitude in modes and Period for an ideal beam} \label{table3} \end{table} \begin{table}[ht] \begin{center} \begin{scriptsize} \begin{tabular}{|c|c|c|c|c|c|} \hline Crack size & $1^{st}$ mode in $Y_-$ & $2^{nd}$ mode in $Y_-$ & $1^{st}$ mode in $Y_+$ & $2^{nd}$ mode in $Y_+$ & Period \\ \hline 1 & 0.0979 & 0.0795 & 0.0848 & 0.0749 & 1.88 \\ \hline 2 & 0.1012 & 0.0833 & 0.0844 & 0.0752 & 1.93 \\ \hline 3 & 0.1049 & 0.0863 & 0.0828 & 0.0751 & 2 \\ \hline 4 & 0.11 & 0.0935 & 0.0819 & 0.0748 & 2.1 \\ \hline 5 & 0.116 & 0.1023 & 0.0813 & 0.0741 & 2.23 \\ \hline 6 & 0.1221 & 0.1157 & 0.0804 & 0.0734 & 2.4 \\ \hline 7 & 0.1544 & 0.1361 & 0.0797 & 0.0722 & 2.74 \\ \hline 8 & 0.209 & 0.1844 & 0.0791 & 0.071 & 3.43 \\ \hline \end{tabular} \end{scriptsize} \end{center} \caption{Amplitude in modes and Periods for a beam with crack at the supporting wall} \label{table4} \end{table} The direction of the applied force is considered as a positive direction. The direction opposite to the applied force direction is considered as a negative direction. In reference Tables \ref{table3} to \ref{table6}, trough represents the amplitude in the negative direction while peak represents the amplitude in the positive direction. Results of Tables \ref{table3} and \ref{table4} are obtained by taking $10$ vertical and $5$ horizontal nodes in the finite element discretization of beams without and with a crack. Table \ref{table3} represents amplitudes of an ideal beam vibrations for first two cycles, period and natural frequency for the ideal beam. Considerable change in frequency, amplitude and period can be observed with increasing crack size. Succeeding Figure \ref{fig:f5} represents amplitude versus crack size graph which indicates change in the amplitude in positive and negative direction. As crack size decreases from $8$ to $0$, amplitude in both modes in each direction converges to the amplitude of an ideal beam. Note that crack size $0$ refers to am amplitude of ideal beam vibration in respective mode. Solid curves in both sub figures represent amplitude in $1^{st}$ mode in each direction while dotted curves in both sub figures represent amplitude in $2^{nd}$ mode in each direction. \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig5a} % {CryN0.eps}} \includegraphics[width=0.45\textwidth]{fig5b} % {CryP.eps}} \end{center} \caption{Change in amplitude in first two modes in $Y_-$ and $Y_+$.} \label{fig:f5} \end{figure} Table \ref{table4} represents amplitudes of cracked beam vibrations for the first two cycle, periods and natural frequencies for different crack lengths. Notice that trough values in the first and second cycles are directly proportional to the crack size. In other words, increasing crack size lead to increase in trough values in respective cycles. Comparing first and second trough values of an ideal beam with a cracked beam, we notice an increment in respective values. This indicates, that the ideal beam constitutes least trough values in the first and second cycles compare to the cracked beam. One can consider other cycles as well. We do not include them since observable changes occur in the first two cycles for reasonably small but noticeable deflections. Comparing the first and second peak values of an ideal beam with a cracked beam, we notice a decrement in respective peak values. That indicates, that the ideal beam has higher peak values in respective cycles compared to the cracked beam. Change in the period and amplitude of the vibration are directly proportional to the crack size. All properties mentioned above hold true in case of finer finite element discretizations. To make a more concrete conclusion, we ran simulations for $20$ vertical and $10$ horizontal nodes for $5000$ time-steps and displayed the results for the ideal beam and the cracked beam with the crack located at the supporting wall. All parameter values are kept the same for simplicity. Tables \ref{table5} and \ref{table6} describe the amplitudes of beam vibrations for the first two cycles, period, natural frequency for the ideal beam, and the beam with crack located at the supporting wall respectively for more nodes. \begin{table}[ht] \begin{center} \begin{small} \begin{tabular}{|c|c|c|c|c|} \hline $1^{st}$ mode in $Y_-$ & $2^{nd}$ mode in $Y_-$ & $1^{st}$ mode in $Y_+$ & $2^{nd}$ mode in $Y_+$ & Period \\ \hline 0.598 & 0.4589 & 0.5269 & 0.4552 & 5.49 \\ \hline \end{tabular} \end{small} \end{center} \caption{Amplitude in modes and Period for an ideal beam} \label{table5} \end{table} \begin{table}[ht] \begin{center} \begin{scriptsize} \begin{tabular}{|c|c|c|c|c|c|} \hline Crack size & $1^{st}$ mode in $Y_-$ & $2^{nd}$ mode in $Y_-$ & $1^{st}$ mode in $Y_+$ & $2^{nd}$ mode in $Y_+$ & Period \\ \hline 1 & 0.6126 & 0.4588 & 0.5272 & 0.4557 & 5.51 \\ \hline 2 & 0.6144 & 0.4618 & 0.5261 & 0.4551 & 5.55 \\ \hline 3 & 0.62 & 0.4708 & 0.5252 & 0.4527 & 5.61 \\ \hline 4 & 0.6332 & 0.4785 & 0.5209 & 0.4509 & 5.68 \\ \hline 5 & 0.6385 & 0.4971 & 0.5149 & 0.4497 & 5.78 \\ \hline 6 & 0.6592 & 0.5093 & 0.5066 & 0.4451 & 5.89 \\ \hline 7 & 0.6673 & 0.5303 & 0.4924 & 0.4324 & 6.02 \\ \hline 8 & 0.6888 & 0.5564 & 0.486 & 0.4272 & 6.18 \\ \hline 9 & 0.7042 & 0.5734 & 0.4784 & 0.4257 & 6.36 \\ \hline 10 & 0.7245 & 0.5894 & 0.4676 & 0.4236 & 6.55 \\ \hline 11 & 0.7516 & 0.6071 & 0.4657 & 0.422 & 6.73 \\ \hline 12 & 0.772 & 0.6501 & 0.4575 & 0.4145 & 6.97 \\ \hline 13 & 0.8497 & 0.7284 & 0.4504 & 0.4106 & 7.29 \\ \hline 14 & 0.9513 & 0.741 & 0.4468 & 0.4073 & 7.9 \\ \hline 15 & 1.063 & 0.814 & 0.4460 & 0.3886 & 8.52 \\ \hline 16 & 1.196 & 0.9208 & 0.4419 & 0.3847 & 9.44 \\ \hline 17 & 1.392 & 1.068 & 0.4397 & 0.3785 & 10.89 \\ \hline 18 & 1.853 & 1.392 & 0.433 & 0.3545 & 14.23 \\ \hline \end{tabular} \end{scriptsize} \end{center} \caption{Amplitude in modes and Periods for a beam with crack at the supporting wall} \label{table6} \end{table} Succeeding Figure \ref{fig:f6} represents amplitude versus crack size graph for finer finite element discretization which also indicates change in the amplitude in positive and negative direction. Numerical convergence in amplitude of vibrations can be observed as crack size decreases from $18$ to $0$. Solid curves in both sub figures represent amplitude in $1^{st}$ mode in each direction while dotted curves in both sub figures represent amplitude in $2^{nd}$ mode in each direction. \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig6a} % CryN.eps \includegraphics[width=0.45\textwidth]{fig6b} % CryP1.eps \end{center} \caption{Change in amplitude in first two modes in $Y_-$ and $Y_+$.} \label{fig:f6} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.48\textwidth]{fig7a} % fg6a.eps \includegraphics[width=0.48\textwidth]{fig7b} \\ % fg6b.eps} \includegraphics[width=0.48\textwidth]{fig7c} % fg6c.eps} \includegraphics[width=0.48\textwidth]{fig7d} \\ % fg6d.eps} \includegraphics[width=0.48\textwidth]{fig7e} % fg6e.eps} \includegraphics[width=0.48\textwidth]{fig7f} \\ % fg6f.eps} \includegraphics[width=0.48\textwidth]{fig7g} % fg6g.eps} \end{center} \caption{Movement of nodes on the crack during vibrations.} \label{fig:f7} \end{figure} Results in Figure \ref{fig:f7} are obtained by dividing the beam into $5$ horizontal nodes and $8$ vertical nodes. Length of the crack is 6. Figure \ref{fig:f7} represents movement of nodes living on the crack where unilateral contact conditions are satisfied. In all sub-figures, the blue curve denotes variation in the $y$ coordinate and the red curve denotes variation in the $x$ coordinate of the top node of the crack. First sub-figure on the upper left corner refers to the movement of the top node of the crack. Similarly, the first figure in the upper right corner is the second node on the crack from the top. The next four sub-figures refer to displacements of the next four consecutive nodes. Since the crack length is 6, there are six nodes which are not connected to the supporting wall. Opening and closing of the crack is assumed in such a way that the first node closest to the crack tip comes in contact with the wall and then second node closest to the crack should come in contact, \dots etc. The top node of the crack comes in contact with the supporting wall at last. In other words, the crack is closed when the last node comes in contact with the supporting wall. First six sub-figures are due to the solution of the linear complementarity problem at the crack location. From all these sub-figures, it is clear that the non-penetration and unilateral conditions are satisfied at all crack nodes. The last sub-figure refers to the crack tip which is connected with the supporting wall. It does not move under the applied force so both the blue and red curves take the value zero for each time step. \subsection*{Conclusions} We discretize a standard two dimensional cantilever beam using finite element method. We model contact boundary conditions for a cantilever beam with normal crack located at the supporting wall. Discretization of the unilateral contact type boundary conditions satisfied at the normal crack location lead us to a linear complementarity problem satisfied at the crack location. We have developed an effective method to find the solution of the linear complementary problem. Using numerical results, we compared vibrations of a beam containing a normal crack at the supporting wall with an ideal beam, and noticed change in the amplitude, frequency and period of vibrations. Size of a crack is a major factor which play crucial role in change of amplitude, period and frequency during vibration of a cracked beam. Amplitude, frequency and period of the beam vibration are influenced due to the presence of crack. Convergence of this model with actual 3 dimensional beam model will be addressed in our future paper. \begin{thebibliography}{00} \bibitem{bs89} M. Schatzman, M. Bercovier; \emph{Numerical approximation of a wave equation with unilateral constraints}, Math. of Computation, 53, (1989), 55-79. \bibitem{CPS92} R. W. Cottle, J-S Pang, R. E. Stone; \emph{The Linear Complementarity Problem}, Academic Press, Boston, (1992). \bibitem{CJM05} C. Eck, J. Jarusek, M. Krbec; \emph{Unilateral Contact Problems; Variational Methods and Existence Theorems}, CRC Press, Boca Raton,(2005). \bibitem{dt95} W. Desch, J. Turi; \emph{Hyperbolic two-point contact problems and a differential inclusion}, SRC Technical Report, Karl Franzens University of Graz, 48, (1995). \bibitem{dt98} W. Desch, J. Turi; \emph{The stop operator related to a convex polyhedron}, Journal of Differential Equations, 157, (1999), 329-347. \bibitem{dim96} A. D. Dimarogonas; \emph{Vibration of cracked strucures: A state of the art review}, Engrg. Fracture Mechanics, 55, (1996), 831-857. \bibitem{SSB03} Z. Belhachmi, J. M. Sac-Epee, J. Sokolowski; \emph{Mixed Finite Element Method For Smooth Domain Formulation of Crack Problem}, SIAM J. Numer. Anal., Vol. 43, No. 3, (2003), pp. 1295-1320. \bibitem{glo99} Ch. Glocker; \emph{Formulation of spatial contact situations in rigid multibody systems}, Comput. Methods Appl. Mech. Engrg., 177, (1999), 199-214. \bibitem{KO88} N. Kikuchi, J. T. Oden; \emph{Contact Problems in Elasticity, A Study of Variational Inequalities and Finite Element Methods}, SIAM, Philadelphia, (1988). \bibitem{KB96} Y. W. Kwon, H. Bang; \emph{The Finite Element Method Using MATLAB}, CRC Press, Boca Raton, (1996). \bibitem{ls84} G. Lebeau, M. Schatzman; \emph{A wave problem in a half-space with a unilateral constraint}, J. Differential Equations, 53, (1984), 309-361. \bibitem{sfe01} J. K. Sinha, M. I. Friswell, S. Edwards; \emph{Simplified models for the location of cracks in beam structures using measured vibration data ,Journal of Sound and Vibration}, 251, (2001), 13-18. \bibitem{DL76} G. Duvaut, J. L. Lions; \emph{Inequalities in Mechanics and Physics}, Springer, New York, (1976). \bibitem{CD98} T. G. Chondros, A.D. Dimarogonas; \emph{Vibration of a Cracked Cantilever Beam}, Journal of Vibration and Acoustics, 742/vol 120, July 1998. \bibitem{PD89} N. Papaeconomou, A. Dimarogonas; \emph{Computational Mechanics}, Springer-Verlag, (1989) 5, 88-94. \end{thebibliography} \end{document}