text:
Nonlinear finite elements/Updated Lagrangian approach

Updated Lagrangian Approach
For the total Lagrangian approach, the discrete equations are formulated with respect to the reference configuration. For the updated Lagrangian approach, the discrete equations are formulated in the  current configuration, which is assumed to be the new  reference configuration.

The independent variables in the total Lagrangian approach are $$X$$ and $$t$$. In the updated Lagrangian they are $$x$$ and $$t$$ which are with respect to the new reference configuration.

The dependent variable in the total Lagrangian approach is the displacement $$u(X,t)$$. In the updated Lagrangian approach, the dependent variables are the Cauchy stress $$\sigma(x,t)$$ and the velocity $$v(x,t)$$.

Updated Lagrangian Stress and Strain Measures
The stress measure is the  Cauchy stress:

{ \sigma(x,t) = \cfrac{T}{A} ~. } $$ The strain measure is the  rate of deformation (also called velocity strain):

{ D(x,t) = \frac{\partial v}{\partial x} = v_{,x} ~. } $$ Note that the derivative is a spatial derivative. It can be shown that

\int_0^t D(x,\tau)~d\tau = \ln(F(x,t)) $$ where $$F$$ is the deformation gradient. So the integral of the rate of deformation in 1-D is similar to the logarithmic strain (also called natural strain).

Governing Equations in Updated Lagrangian Form
The  Updated Lagrangian governing equations are:

Conservation of Mass

{ \rho~J = \rho_0~. } $$ For the rod,

\rho~A~F = \rho_0~A_0 ~. $$ These are the same as those for the total Lagrangian approach.

Conservation of Momentum

{ \frac{\partial }{\partial x}(A~\sigma) + \rho~A~b = \rho~A_0~\cfrac{Dv}{Dt} ~. } $$ In this case $$A$$ may not be constant at along the length and further simplification cannot be done. In short form,

(A\sigma)_{,x} + \rho~A~b = \rho~A~\dot{v} ~. $$

Conservation of Energy
If we ignore heat flux and heat sources

{ \rho\frac{\partial w_{\text{int}}}{\partial t} = \sigma~D~. } $$ If we include heat flux and heat sources

{ \rho\frac{\partial w_{\text{int}}}{\partial t} = \sigma~D - \frac{\partial q}{\partial x} + \rho~s~. } $$ In short form:

\rho\dot{w}_{\text{int}} = \sigma~D - q_{,x} + \rho~s~. $$

Total Form
For a linear elastic material:

{ \sigma(x,t) = E^{\sigma D} \int_0^t D(x,\tau)~d\tau = E^{\sigma D} \ln(F(x,t)) ~. } $$ The superscript $$\sigma D$$ refers to the fact that this function relates $$\sigma$$ and $$D$$.

For small strains:

E^{\sigma D} = \text{Youngs modulus}~. $$

Rate Form
For a linear elastic material:

{ \dot{\sigma}(x,t) = E^{\sigma D}~D(x,t) ~. } $$

Initial and Boundary Conditions for Updated Lagrangian
For the updated Lagrangian approach, initial conditions are needed for the stress and the velocity. The initial displacement is assumed to be zero.

The initial conditions are:

{ \begin{align} \sigma(x,0) & = \sigma_0(x) & ~\text{for}~ & x \in [0, L] \\ v(x,0) & = v_0(x) & ~\text{for}~ & x \in [0, L]\\ u(x,0) & = 0 & ~\text{for}~ & x \in [0, L] \end{align} } $$

The  essential boundary conditions are

{ v(x,t) = \bar{v}(x,t) \qquad ~\text{for}~x \in \Gamma_v ~. } $$ The  traction boundary conditions are

{ n~\sigma(x,t) = t_x(x,t) \qquad ~\text{for}~x \in \Gamma_t ~. } $$ The  unit normal to the boundary in the current configuration is $$n$$. The tractions in the current configuration are related to those in the reference configuration by

{ t_x~A = t_x^0~A_0 ~. } $$

The  interior continuity or  jump condition is

{ \lfloor\sigma~A\rceil = 0 ~. } $$

Weak Form for Updated Lagrangian
The momentum equation is

(A~\sigma)_{,x} + \rho~A~b = \rho~A~\cfrac{Dv}{Dt} ~. $$ To get the weak form over an element, we multiply the equation by a weighting function and integrate over the  current length of the element (from $$x_a$$ to $$x_b$$).

\int_{x_a}^{x_b} w\left[(A~\sigma)_{,x} + \rho~A~b - \rho~A~\cfrac{Dv}{Dt}\right]~dx = 0 ~. $$ Integrate the first term by parts to get

\int_{x_a}^{x_b} w (A~\sigma)_{,x}~dx = \left[w~A~\sigma\right]_{x_a}^{x_b} - \int_{x_a}^{x_b} w_{,x}~(A~\sigma)~dx ~. $$ Plug the above into the weak form and rearrange to get

\int_{x_a}^{x_b} w_{,x}~(A~\sigma)~dx + \int_{x_a}^{x_b} w~\rho~A~\cfrac{Dv}{Dt}~dx = \int_{x_a}^{x_b} w~\rho~A~b~dx + \left[w~A~\sigma\right]_{x_a}^{x_b} ~. $$ If we think of the weighting function $$w$$ as a variation of $$v$$ that satisfies the kinematic admissibility conditions, we get the  principle of virtual power:

{ \int_{x_a}^{x_b} (\delta v)_{,x} (A~\sigma)~dx + \int_{x_a}^{x_b}\delta v~\rho~A~\cfrac{Dv}{Dt}~dx = \int_{x_a}^{x_b} \delta v~\rho~A~b~dx + \left[\delta v~A~\sigma\right]_{x_a}^{x_b} ~. } $$ Recall that

v_{,x} = D(x,t) \qquad\text{and}\qquad t_x = n~\sigma ~. $$ Therefore,

\delta v_{,x} = \delta D(x,t) $$ and

\left[\delta v~A~\sigma\right]_{x_a}^{x_b} = \left[\delta v~A~t_x\right]_{\Gamma_t} ~. $$ Hence we can alternatively write the weak form as

{ \int_{x_a}^{x_b}\delta D~A~\sigma~dx + \int_{x_a}^{x_b}\delta v~\rho~A~\cfrac{Dv}{Dt}~dx = \int_{x_a}^{x_b}\delta v~\rho~A~b~dx +\left[\delta v~A~t_x\right]_{\Gamma_t}~. } $$ Comparing with the energy equation, we see that the first term above is the internal virtual power. Using the substitutions $$A~dx = d\Omega$$ and $$Dv/DT = \dot{v}$$, we can also write the weak form as

{ \int_{\Omega}\delta D~\sigma~d\Omega + \int_{\Omega}\delta v~\rho~\dot{v}~d\Omega = \int_{\Omega}\delta v~\rho~b~d\Omega +\left[\delta v~A~t_x\right]_{\Gamma_t}~. } $$ The weak form may also be written in terms of the virtual powers as

{ \delta P = \delta P_{\text{int}} - \delta P_{\text{ext}} + \delta P_{\text{kin}} = 0 } $$ where,

{ \begin{align} \delta P_{\text{int}} & = \int_{x_a}^{x_b} (\delta v)_{,X}~(A~\sigma)~dx \\ \delta P_{\text{ext}} & = \int_{x_a}^{x_b} \delta v~\rho~A~b~dx + \left[\delta v~A~t_X\right]_{\Gamma_t} \\ \delta P_{\text{kin}} & = \int_{x_a}^{x_b} \delta v~\rho~A~\dot{v}~dx \end{align} } $$

Finite Element Discretization for Updated Lagrangian
The  trial velocity field is

{ v_h(x,t) = \sum_{j=1}^n N_j(x) v_j(t) } $$ In matrix form,

{ v_h(x,t) = \mathbf{N}(x)~\mathbf{v}(t) } \text{where} \mathbf{N} = \begin{bmatrix} N_1(X) & N_2(X) && N_n(X) \end{bmatrix} ~\text{and}~ \mathbf{v} = \begin{bmatrix} v_1(t) \\ v_2(t) \\ \vdots \\ v_n(t) \end{bmatrix} ~. $$ The resulting  acceleration field is the  material time derivative of the velocity

a_h(x,t) = \dot{v}_h(x,t) = \sum_{j=1}^n N_j(x) \dot{v}_j(t) ~. $$ Hence,

{ a_h(x,t) = \sum_{j=1}^n N_j(x) a_j(t) } ~\text{or}~ { a_h(X,t) = \mathbf{N}(x)~\mathbf{a}(t) } \text{where} \mathbf{a} = \begin{bmatrix} a_1(t) \\ a_2(t) \\ \vdots \\ a_n(t) \end{bmatrix} ~. $$ Note that the shape functions are functions of $$X$$ and not of $$x$$. We could transform them into function of $$x$$ using the inverse mapping

x = \varphi^{-1}(x, t) ~. $$ However, in that case the shape functions become functions of time and the procedure becomes more complicated.

The  test (weighting) function is

{ \delta v_h(x) = \sum_{i=1}^n N_i(x) \delta v_i ~. } \text{or} { \delta v_h(x) = \mathbf{N}~\delta \mathbf{v} } \text{where} \delta \mathbf{v} = \begin{bmatrix} \delta v_1 \\ \delta v_2 \\ \vdots \\ \delta v_n \end{bmatrix} ~. $$ The  derivatives of the test functions with respect to $$x$$ are

{ (\delta v_h)_{,x} = \sum_{i=1}^n N_{i,x} \delta v_i ~. } \text{or} { (\delta v_h)_{,x} = \mathbf{B}~\delta \mathbf{v} } \text{where} \mathbf{B} = \begin{bmatrix} N_{1,x} & N_{2,x} && N_{n,x} \end{bmatrix} ~. $$ It is convenient to use the isoparametric concept to compute the spatial derivatives. Let us reexamine what this approach involves.

Figure 1 shows a two-noded 1-D element in the reference and current configurations along with the parent element.

The  map between the Eulerian coordinates and the parent coordinates is

x(\xi, t) = (1-\xi)~x_1^e(t) + \xi~x_2^e(t)~\xi = N_1(\xi)~x_1^e(t) + N_2(\xi)~x_2^e(t)\qquad \text{or} { x(\xi, t) = \mathbf{N}(\xi)~\mathbf{x}^e(t) ~. } $$ At $$t = 0$$,

x(\xi, 0) = x(\xi) = \mathbf{N}(\xi)~\mathbf{X}^e ~. $$ Therefore the  displacement in the parent coordinates is

u(\xi, t) = x(\xi, t) - x(\xi) = \mathbf{N}(\xi)~\left[\mathbf{x}^e(t) - \mathbf{X}^e\right] \qquad \text{or} { u(\xi, t) = \mathbf{N}(\xi)~\mathbf{u}^e(t) ~. } $$ Similarly, the  velocity and its variation in the parent coordinates are given by

{ v(\xi, t) = \mathbf{N}(\xi)~\mathbf{v}^e(t) ~,\qquad \delta v(\xi, t) = \mathbf{N}(\xi)~\delta \mathbf{v}^e(t) ~. } $$ The  acceleration in the parent coordinates is given by

{ a(\xi, t) = \mathbf{N}(\xi)~\dot{\mathbf{v}}^e(t) ~. } $$ The  rate of deformation is given by

D(x, t) = v_{,x} (x, t) = \mathbf{N}_{,x}(X(x,t))~\mathbf{v}^e(t)~. $$ We can evaluate the derivative with respect to $$x$$ by simply using the map to the parent coordinate instead of mapping back to the reference coordinates. Using the chain rule (for the two noded element)

N_{1,\xi} = N_{1,x}~x_{,\xi} \qquad\text{and}\qquad N_{2,\xi} = N_{2,x}~x_{,\xi}~. $$ In matrix form,

{ \mathbf{N}_{,\xi} = \mathbf{N}_{,x}~x_{,\xi} } \qquad\implies\qquad { \mathbf{N}_{,x} = \mathbf{N}_{,\xi}~\left(x_{,\xi}\right)^{-1} =: \mathbf{B}(\xi)~. } $$ The  rate of deformation  in parent coordinates is then given by

{ D(\xi, t) = \mathbf{B}(\xi)~\mathbf{v}^e(t)~. } $$

To derive the finite element system of equations, we follow the usual approach of substituting the trial and test functions into the approximate weak form

\int_{x_a}^{x_b} (\delta v_h)_{,x} (A~\sigma)~dx + \int_{x_a}^{x_b}\delta v_h~\rho~A~\cfrac{Dv_h}{Dt}~dx = \int_{x_a}^{x_b} \delta v_h~\rho~A~b~dx + \left[\delta v_h~A~\sigma\right]_{x_a}^{x_b} ~. $$ Let us proceed term by term.

First LHS Term
The first term represents the  virtual internal power

\delta P_{\text{int}} = \int_{x_a}^{x_b} (\delta v_h)_{,x} (A~\sigma)~dx ~. $$ Plugging in the derivative of the test function, we get

\begin{align} \delta P_{\text{int}} & = \int_{x_a}^{x_b} \left[\sum_{i=1}^n N_{i,x}\delta v_i\right] (A~\sigma)~dx\\ & = \sum_{i=1}^n \delta v_i \left[\int_{x_a}^{x_b} N_{i,x}(A~\sigma)~dx\right]\\ & = \sum_{i=1}^n \delta v_i f_i^{\text{int}} ~. \end{align} $$ The matrix form of the expression for virtual internal work is

{ \delta P_{\text{int}} = \delta \mathbf{v}^T \mathbf{f}_{\text{int}} } \text{where} \mathbf{f}_{\text{int}} = \begin{bmatrix} f_1^{\text{int}} \\ f_2^{\text{int}} \\ \vdots \\ f_n^{\text{int}} \end{bmatrix} ~. $$ The internal force is

{ f_i^{\text{int}} = \int_{x_a}^{x_b} N_{i,x}(A~\sigma)~dx = \int_{-1}^{1} N_{i,\xi}(x_{,\xi})^{-1}(A~\sigma)~x_{,\xi}~d\xi = \int_{-1}^{1} N_{i,\xi}(A~\sigma)~d\xi ~. } $$ Note that the above simplification only occurs in 1-D. In matrix form,

{ \mathbf{f}_{\text{int}} = \int_{x_a}^{x_b} \mathbf{B}^T~(A~\sigma)~dx = \int_{-1}^{1} (\mathbf{N}_{,\xi})^T~(A~\sigma)~d\xi~. }\text{or} { \mathbf{f}_{\text{int}} = \int_{\Omega} \mathbf{B}^T~\sigma~d\Omega ~. } $$

 Remark:

Note that we can write

N_{,x} = N_{,X}~\cfrac{dX}{dx} \qquad \implies \qquad N_{,x}~dx = N_{,X}~dX ~. $$ If we use the above relation, we get

\mathbf{f}_{\text{int}} = \int_{x_a}^{x_b} (\mathbf{N}_{,x})^T~\sigma~A~dx = \int_{X_a}^{X_b} (\mathbf{N}_{,X})^T~\sigma~A~dX ~. $$ Using the relation

\sigma~A = P~A_0 $$ we get

{ \mathbf{f}_{\text{int}} = \int_{X_a}^{X_b} (\mathbf{N}_{,X})^T~P~A_0~dX ~. } $$ The above is the same as the expression we had for the internal force in the total Lagrangian formulation.

Second LHS Term
The second term represents the  virtual kinetic power

\delta P_{\text{kin}} = \int_{x_a}^{x_b}\delta v_h~\rho~A~\cfrac{Dv_h}{Dt}~dx ~. $$ Plugging in the test function, we get

\begin{align} \delta P_{\text{kin}} & = \int_{x_a}^{x_b}\left[\sum_{i=1}^n \delta v_i N_i\right] \rho~A~\cfrac{Dv_h}{Dt}~dx\\ & = \sum_{i=1}^n \delta v_i \left[\int_{x_a}^{x_b} \rho~A~N_i~\cfrac{Dv_h}{Dt}~dx\right]\\ & = \sum_{i=1}^n \delta v_i f_i^{\text{kin}} ~. \end{align} $$ The inertial (kinetic) force is

{ f^{\text{kin}}_i = \int_{x_a}^{x_b}\rho~A~N_i~\cfrac{Dv_h}{Dt}~dx ~. } $$ Now, plugging in the trial function into the expression for the inertial force, we get

\begin{align} f^{\text{kin}}_i & = \int_{x_a}^{x_b}\rho~A~N_i \left[\sum_{j=1}^n \cfrac{Dv_j}{Dt} N_j\right]~dx \\ & = \sum_{j=1}^n \left[\int_{x_a}^{x_b}\rho~A~N_i~N_j~dx\right] \cfrac{Dv_j}{Dt} \\ & = \sum_{j=1}^n M_{ij} a_j \end{align} $$ where

{ M_{ij} := \int_{x_a}^{x_b}\rho~A~N_i~N_j~dx \qquad \leftarrow \quad \text{Consistent Mass Matrix Coefficient.} } $$ and

a_j := \cfrac{Dv_j}{Dt} = \dot{v}_j \qquad\leftarrow \quad\text{Acceleration.} $$ In matrix form,

{ \mathbf{f}_{\text{kin}} = \mathbf{M} \mathbf{a} } \text{where} \mathbf{a} = \begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} \dot{v}_1 \\ \dot{v}_2 \\ \vdots \\ \dot{v}_n \end{bmatrix} \text{and} \mathbf{f}_{\text{kin}} = \begin{bmatrix} f_1^{\text{kin}} \\ f_2^{\text{kin}} \\ \vdots \\ f_n^{\text{kin}} \end{bmatrix} ~. $$ The consistent mass matrix in matrix notation is

{ \mathbf{M} =\int_{x_a}^{x_b}\rho~A~\mathbf{N}^T~\mathbf{N}~dx = \int_{\Omega}\rho~\mathbf{N}^T~\mathbf{N}~d\Omega ~. } $$ Plugging the expression for the inertial force into the expression for virtual kinetic power, we get

\begin{align} \delta P_{\text{kin}} & = \sum_{i=1}^n \delta v_i \left[\sum_{j=1}^n M_{ij} a_j\right] \\ & = \sum_{i=1}^n \sum_{j=1}^n \delta v_i M_{ij} a_j \end{align} $$ In matrix form,

{ \delta P_{\text{kin}} = (\delta \mathbf{v})^T \mathbf{M} \mathbf{a} = (\delta \mathbf{v})^T~\mathbf{f}_{\text{kin}} ~. } $$

 Remark:

Note that since the integration domain is a function of time, the mass matrix for the updated Lagrangian formulation is also a function of time. However, from the conservation of mass, we have

\rho_0~A_0~dX = \rho~A~dx ~. $$ Therefore, we can write the mass matrix as

{ \mathbf{M} =\int_{x_a}^{x_b}\rho~A~\mathbf{N}^T~\mathbf{N}~dx = \int_{X_a}^{X_b}\rho_0~A_0~\mathbf{N}^T~\mathbf{N}~dX = \int_{\Omega_0}\rho_0~\mathbf{N}^T~\mathbf{N}~d\Omega_0 ~. } $$ Hence the mass matrices are the same for both formulations.

RHS Terms
The right hand side terms represent the  virtual external power

\delta P_{\text{ext}} = \int_{x_a}^{x_b}\delta v_h~\rho~A~b~dx + \left[\delta v_h~A~t_x\right]_{\Gamma_t}~. $$ Plugging in the test function into the above expression gives

\begin{align} \delta P_{\text{ext}} & = \int_{x_a}^{x_b}\left[\sum_{i=1}^n N_i~\delta v_i\right]\rho~A~b~dx + \left[\left(\sum_{i=1}^n N_i~\delta v_i\right) A~t_x\right]_{\Gamma_t} \\ & = \sum_{i=1}^n \delta v_i\left[\int_{x_a}^{x_b} N_i~\rho~A~b~dx\right] + \left[\sum_{i=1}^n \delta v_i~N_i~A~t_X\right]_{\Gamma_t} \\ & = \sum_{i=1}^n \delta v_i\left[\int_{x_a}^{x_b} N_i~\rho~A~b~dx + \left[N_i~A~t_x\right]_{\Gamma_t}\right] \\ & = \sum_{i=1}^n \delta v_i~f^{\text{ext}}_i ~. \end{align} $$ In matrix notation,

{ \delta P_{\text{ext}} = \delta \mathbf{v}^T~\mathbf{f}_{\text{ext}} } \text{where} \mathbf{f}_{\text{ext}} = \begin{bmatrix} f_1^{\text{ext}} \\ f_2^{\text{ext}} \\ \vdots \\ f_n^{\text{ext}} \end{bmatrix} ~. $$ The external force is given by

{ f_i^{\text{ext}} =\int_{x_a}^{x_b} N_i~\rho~A~b~dx + \left[N_i~A~t_x\right]_{\Gamma_t}~. } $$ In matrix notation,

{ \mathbf{f}_{\text{ext}} = \int_{x_a}^{x_b} \mathbf{N}^T~\rho~A~b~dx + \left[\mathbf{N}^T~A~t_x\right]_{\Gamma_t} }\text{or} { \mathbf{f}_{\text{ext}} = \int_{\Omega} \rho~\mathbf{N}^T~b~d\Omega + \left[\mathbf{N}^T~A~t_x\right]_{\Gamma_t}~. } $$

 Remark:

Using the conservation of mass

\rho_0~A_0~dx = \rho~A~dx $$ and the relations between the traction sin the reference and the current configurations

A_0~t_x^0 = A~t_x $$ we can transform the integral in the expression for the external force to one over the reference coordinates as follows:

\mathbf{f}_{\text{ext}} = \int_{x_a}^{x_b} \mathbf{N}^T~\rho_0~A_0~b~dx + \left[\mathbf{N}^T~A_0~t_x^0\right]_{\Gamma_t}~. $$ The above is the same as the expression we had for the external force in the total Lagrangian formulation.

Discrete Equations for Updated Lagrangian
The finite element equations in updated Lagrangian form are:

{ \begin{align} v(x,t) & = \sum_i N_i(x)~v_i^e(t) = \mathbf{N}~\mathbf{v}^e \\ D(x,t) & = \sum_i N_{i,x}~v_i^e(t) = \mathbf{B}~\mathbf{v}^e \\ \mathbf{f}^e_{\text{int}} & = \int_{\Omega^e} \mathbf{B}^T~\sigma~d\Omega \\ \mathbf{f}^e_{\text{ext}} & = \int_{\Omega^e} \rho~\mathbf{N}^T~b~d\Omega + \left[\mathbf{N}^T~A~t_x \right]_{\Gamma_t^e} \\ \mathbf{M}^e & = \int_{\Omega_0^e}~\rho_0~\mathbf{N}^T~\mathbf{N}~d\Omega_0 \\ \mathbf{M}~\ddot{\mathbf{u}} & = \mathbf{f}_{\text{ext}} - \mathbf{f}_{\text{int}} ~. \end{align} } $$


license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

title: File:Config1D.png


text:
Nonlinear finite elements/Effect of mesh distortion

An Example: Effect of Mesh Distortion
Consider the three-noded quadratic displacement element shown in Figure 1.

The shape functions for the parent element are
 * $$\begin{align}

N_1(\xi) & = \frac{1}{2}\xi(\xi-1) \\ N_2(\xi) & = 1-\xi^2 \\ N_3(\xi) & = \frac{1}{2}\xi(\xi+1) \end{align}$$ In matrix form,

\mathbf{N}(\xi) = \begin{bmatrix} \frac{1}{2}\xi(\xi-1) & 1-\xi^2 & \frac{1}{2}\xi(\xi+1) \end{bmatrix} ~. $$ The trial functions (in terms of the parent coordinates) are
 * $$\begin{align}

u_h(\xi, t) & = \mathbf{N}(\xi)~\mathbf{u}(t) \\ v_h(\xi, t) & = \mathbf{N}(\xi)~\mathbf{v}(t) \\ a_h(\xi, t) & = \mathbf{N}(\xi)~\mathbf{a}(t) ~. \end{align}$$ The mapping from the Eulerian coordinates to the parent element coordinates is

x(\xi, t) = N_1(\xi)~x_1(t) + N_2(\xi)~x_2(t) + N_3(\xi)~x_3(t) ~. $$ In matrix form,

x(\xi, t) = \mathbf{N}(\xi)~\mathbf{x}(t)~. $$ Therefore, the derivative with respect to $$\xi$$ is

x_{,\xi} = N_{1,\xi}~x_1(t) + N_{2,\xi}~x_2(t) + N_{3,\xi}~x_3(t) = \frac{1}{2}(2\xi-1)~x_1 - 2\xi~x_2 + \frac{1}{2}(2\xi+1)~x_3 ~. $$ In matrix form,

x_{,\xi} = \mathbf{N}_{,\xi}~\mathbf{x}(t) = \begin{bmatrix} \frac{1}{2}(2\xi-1) & - 2\xi & \frac{1}{2}(2\xi+1) \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}~. $$ Now, the $$\mathbf{B}$$ matrix is given by

\mathbf{B} = \mathbf{N}_{,\xi} (x_{,\xi})^{-1} ~. $$ Therefore,

\mathbf{B} = \cfrac{1}{2 x_{,\xi}} \begin{bmatrix} 2\xi-1 & - 4\xi & 2\xi+1 \end{bmatrix} ~. $$ The rate of deformation is then given by

{ D = \mathbf{N}_{,x}~\mathbf{v}(t) = \mathbf{B}~\mathbf{v}(t) = \cfrac{1}{2 x_{,\xi}} \begin{bmatrix} 2\xi-1 & - 4\xi & 2\xi+1 \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \\ v_3 \end{bmatrix} ~. } $$ The stress can be calculated using the relation

\sigma = E^{\sigma D}~D ~. $$ The internal forces are given by

\mathbf{f}_{\text{int}} = \int_{x_1}^{x_3} \mathbf{B}^T~\sigma~A~dx ~. $$ Plugging in the expression for $$\mathbf{B}$$, and changing the limits of integration, we get

{ \mathbf{f}_{\text{int}} = \int_{-1}^{1} \cfrac{1}{2 x_{,\xi}} \begin{bmatrix} 2\xi-1 \\ - 4\xi \\ 2\xi+1 \end{bmatrix} ~\sigma~A~x_{,\xi}~d\xi = \frac{1}{2}\int_{-1}^{1}\begin{bmatrix} 2\xi-1 \\ - 4\xi \\ 2\xi+1 \end{bmatrix}~\sigma~A~d\xi ~. } $$ If node $$2$$ is midway between node $$1$$ and node $$3$$,

x_2 = \frac{1}{2}(x_1 + x_3) ~. $$ Then we have,

x_{,\xi} = \frac{1}{2}(2\xi-1)~x_1 - \xi~(x_1+x_3) + \frac{1}{2}(2\xi+1)~x_3 = \frac{1}{2}(x_3-x_1) ~. $$ The rate of deformation becomes

D = \cfrac{1}{x_3-x_1} \begin{bmatrix} 2\xi-1 & - 4\xi & 2\xi+1 \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \\ v_3 \end{bmatrix} $$ which is a linear function of $$\xi$$.

However, if node $$2$$ moves away from the middle during deformation, then $$x_{,\xi}$$ is no longer constant and can become zero or negative. Under such situations the rate of deformation is either infinite or the element inverts upon itself since the isoparametric map is no longer one-to-one.

Let us consider the case where $$x_{,\xi}$$ is zero. In that case, the Jacobian becomes

J = \cfrac{A}{A_0} x_{,X} = \cfrac{A}{A_0} x_{,\xi}~(X_{,\xi})^{-1} = 0 ~. $$ Similarly, when $$x_{,\xi}$$ is negative, $$J$$ is negative. This implies that the conservation of mass is violated.

To find the location of $$x_2$$ when this happens, we set the relation for $$x_{,\xi}$$ to zero. Then we get,

\frac{1}{2}(2\xi-1)~x_1 - \xi~(x_1+x_3) + \frac{1}{2}(2\xi+1)~x_3 = 0 \qquad \implies \qquad x_2 = \cfrac{x_1+x_3}{2} + \cfrac{x_3-x_1}{4\xi}~. $$ If $$x_{,\xi} = 0$$ at $$\xi = 1$$, then

x_2 = \cfrac{x_1+3x_3}{4}~. $$ This means that as node $$2$$ gets closer and closer to node $$3$$, the rate of deformation become infinite at node $$3$$ and then negative.

If $$x_{,\xi} = 0$$ at $$\xi = -1$$, then

x_2 = \cfrac{3x_1+x_3}{4}~. $$ This means that as node $$2$$ gets closer and closer to node $$1$$, the rate of deformation becomes infinite at node $$1$$ and then negative.

These effects of mesh distortion can be severe in multiple dimensions. That is the reason that linear elements are preferred in large deformation simulations.


title: File:UpdateLag3Node.png


text:
Nonlinear finite elements/Solution procedure

Solution procedure
The finite element system of equations is

\mathbf{M}~\ddot{\mathbf{u}} = \mathbf{f}_{\text{ext}} - \mathbf{f}_{\text{int}} = \mathbf{f} ~. $$

Let's assume that the mass matrix is diagonal. Let us also assume that we would like to solve this problem using an explicit method that uses central differencing.

Let the time step size be $$\Delta t$$. Let us also assume that the time step is constant. Then, at time step $$n$$, the time is $$t_n = n\Delta t$$.

Let

\mathbf{u}_n = \mathbf{u}(t_n) ~. $$ Let us take a half step to compute the velocity at the middle of the timestep. Then,

\mathbf{v}^{n+1/2} = \cfrac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\Delta t} ~. $$ We will use this half-step velocity to compute the acceleration

\mathbf{a}^n = \ddot{\mathbf{u}}^n = \cfrac{\mathbf{v}^{n+1/2} - \mathbf{v}^{n-1/2}}{\Delta t} ~. $$ Now,

\ddot{\mathbf{u}}^n = \mathbf{M}^{-1}~\mathbf{f}^n ~. $$ Therefore,

\mathbf{M}^{-1}~\mathbf{f}^n =\cfrac{\mathbf{v}^{n+1/2} - \mathbf{v}^{n-1/2}}{\Delta t} \qquad \implies \qquad \mathbf{v}^{n+1/2} = \mathbf{v}^{n-1/2} + \Delta t~\mathbf{M}^{-1}~\mathbf{f}^n ~. $$

Algorithm

 * 1) Initialize:
 * 2) Set $$t = 0$$ and $$n = 0$$.
 * 3) Set the initial velocities ($$\mathbf{v}^0$$) at the nodes.
 * 4) Set the initial accelerations ($$\mathbf{a}^0$$) at the nodes.
 * 5) Set the initial stresses ($$\sigma^0$$) at the element Gauss points.
 * 6) Compute the lumped mass $$\mathbf{M}$$ at the nodes.
 * 7) Set the displacements and velocities at the first time step:
 * 8) Displacement: $$ \mathbf{u}^0 = 0 ~.$$
 * 9) Velocity: $$ \mathbf{v}^{1/2} = \mathbf{v}^0 + \cfrac{\Delta t}{2}~\mathbf{a}^0 ~. $$
 * 10) Apply essential BCs: $$ [\mathbf{v}^{1/2}]_{\Gamma_v} = \bar{\mathbf{v}}(\mathbf{x}, \Delta t/2) ~. $$
 * 11) Update the nodal displacements: $$ \mathbf{u}^1 = \mathbf{u}^0 + \Delta t~\mathbf{v}^{1/2} ~. $$
 * 12) Set $$n=1$$ and $$t = \Delta t$$.
 * 13) Loop through the following steps until $$t = t_{\text{max}}$$.
 * 14) For each element (at the Gauss points):
 * 15) Compute the strain measure: $$ \varepsilon^n = u_{,X} \text{or} D^n = v_{,x} ~. $$
 * 16) Compute the stress: $$ P = E^{PF}~\varepsilon^n \text{or} \dot{\sigma} = E^{\sigma D}~D^n ~.$$
 * 17) For each node:
 * 18) Compute the internal force $$(f_i^n)_{\text{int}}$$.
 * 19) Compute the external force $$(f_i^n)_{\text{ext}}$$.
 * 20) Compute the total force $$f_i^n$$.
 * 21) Compute the acceleration $$a_i^n$$: $$ a_i^n = M_{ii}^{-1}~f_i^n ~. $$
 * 22) Update the velocity $$v_i^{n+1/2}$$. $$ v_i^{n+1/2} = v_i^{n-1/2} + \Delta t~a_i^n ~. $$
 * 23) Apply the essential boundary conditions.
 * 24) Update the displacement $$u_i^{n+1}$$. $$ u_i^{n+1} = u_i^n + \Delta t~v_i^{n+1/2}~. $$
 * 25) Update the counters: $$n = n+1$$, $$t = t + \Delta t$$.

Stability of the explicit algorithm
If the time step $$\Delta t$$ is too large, the algorithm may not be stable and may give unreasonable results. To provide a check on the time step size, we use the CFL (Courant-Friedrichs-Lewy) condition to determine $$\Delta t$$. This condition (in 1-D) states that

{ \Delta t_{\text{crit}} = \cfrac{l_0}{c_0} } $$

where $$l_0$$ is the  initial  length of the element, and $$c_0$$ is speed of sound in the material (wavespeed) given by

{ c_0 = \cfrac{E^{PF}}{\rho_0} ~. } $$


text:
Nonlinear finite elements/Natural vibration

Special case : natural vibrations
Recall that the finite element system of equations has the form

\mathbf{M}~\ddot{\mathbf{u}} = \mathbf{f}_{\text{ext}} - \mathbf{f}_{\text{int}}~. $$ We could also have written this equation as

\mathbf{M}~\ddot{\mathbf{u}} + \mathbf{K}~\mathbf{u} = \mathbf{f} ~. $$ For natural vibrations, the forces and the displacements are assumed to be periodic in time, i.e.,

\mathbf{u} = \mathbf{u}^0~\exp(i\omega t) ~. $$ and

\mathbf{f} = \mathbf{f}^0~\exp(i\omega t) ~. $$ Then, the accelerations take the form

\ddot{\mathbf{u}} = (i\omega)^2~\mathbf{u}^0~\exp(i\omega t) = -\omega^2~\mathbf{u}^0~\exp(i\omega t) ~. $$ Plugging these into the FE system of equations, we get

[-\omega^2~\exp(i\omega t)]\mathbf{M}~\mathbf{u}^0 + \exp(i\omega t)~\mathbf{K}~\mathbf{u}^0 = \exp(i\omega t)~\mathbf{f}^0 ~. $$ After simplification, we get

{ \left(-\omega^2\mathbf{M} + \mathbf{K}\right)~\mathbf{u}^0 = \mathbf{f}^0 ~. } $$ If there is no forcing, the right hand side is zero, and we get the finite element system of equations for  free vibrations

{ \left(-\omega^2\mathbf{M} + \mathbf{K}\right)~\mathbf{u}^0 =0 ~. } $$ The above equation is similar to the  eigenvalue problem of the form

\mathbf{A}~\mathbf{x} = \lambda~\mathbf{x} \qquad \equiv \qquad \left(\mathbf{A} - \lambda\mathbf{I}\right)\mathbf{x} = 0 ~. $$ Since the right hand side is zero, the finite element system of equations has a solution only if

{ \det(\mathbf{K} -\omega^2\mathbf{M}) =0 ~. } $$ For a two noded element,

\mathbf{K} = \begin{bmatrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{bmatrix} ~\text{and}~ \mathbf{M} = \begin{bmatrix} M_{11} & M_{12} \\ M_{21} & M_{22} \end{bmatrix} ~. $$ Therefore,

\mathbf{K} - \omega^2~\mathbf{M} = \begin{bmatrix} K_{11} -\omega^2~M_{11} & K_{12} -\omega^2~M_{12} \\ K_{21} -\omega^2~M_{21} & K_{22} -\omega^2~M_{22} \end{bmatrix}~. $$ The determinant is

\det(\mathbf{K} - \omega^2~\mathbf{M}) = (K_{11} -\omega^2~M_{11})(K_{22} -\omega^2~M_{22}) - (K_{12} -\omega^2~M_{12})(K_{21} -\omega^2~M_{21})~. $$ This gives us a quadratic equation in $$\omega^2$$ which can be solved to find the natural frequencies of the element.


text:
Nonlinear finite elements/Homework 1

Problem 1: Mathematical Model Development
Consider the motion of a pendulum (see Figure 1).

Assume that:

 * 1) The rod and the mass at the end are rigid.
 * 2) The mass of the rod is negligible compared to the mass of the bob.
 * 3) There is no friction at the pivot.

Find the following:

 * 1) The equation of motion of the bob (angular displacement as a function of time). State all other assumptions.
 * 2) Why is the equation of motion "nonlinear"?
 * 3) Assume that $$\theta$$ is "small". Derive the equation of motion for this situation.
 * 4) Why is the small $$\theta$$ equation of motion "linear"?
 * 5) Derive a finite element formulation for the linear equation of motion.
 * 6) Use ANSYS or some other tool of your choice (including your own code) to solve the linear equation of motion via finite elements. Compare with the exact solution.  You can use any reasonable values for time ($$t$$), mass ($$m$$), and length ($$l$$).
 * 7) What happens when you try to solve the nonlinear equation of motion with finite elements? (You will need to formulate the FE formulation to see the difference.)

Problem 2: Numerical Exercise
In the learning resources section, you have read about the process of formulating a FE model for an axial bar under a distributed body load. Ansys 9 shows you how to solve the same problem using ANSYS. Work the example out on your own.

Next, consider a mini centrifuge of the type manufactured by Kisker Biotech (see Figure 2).

This type of mini centrifuge is widely used in the fields of genetic research and microbiology. The centrifuge is rate to run at 15,000 RPM. The allowable temperature range is -20 C to 40 C.  The nominal radius is 140 mm. Assume other geometric quantities on the basis of the figure.

We want to extend the temperature range to -80 C - 100 C and the allowable RPM to 25,000.

Assume that the current material is a light metal (any other material of your choice will also do). The only requirement is that it should be cheap and stiff. Ignore bending and the weight of the sample tubes.

Find the thinnest beam of the material of your choice that can stand the axial load that is generated due to rotation at the new speed. Use ANSYS or your own code to solve the problem. State all simplifying assumptions.

We just need a ballpark value. A 1-D simulation will do.

title: Nonlinear finite elements/Homework 1


text:
Nonlinear finite elements/Homework 2

Problem 1: Classification of PDEs
Two commonly encountered PDEs in engineering are

u_{,xx} = \alpha~u_{,t} $$ and

u_{,xxxx} = \alpha~u_{,tt} $$

 * 1) What physical processes do the two PDEs model?
 * 2) Write out the PDEs in elementary partial differential notation (that is, using symbols of the form $$\frac{\partial u}{\partial x}$$).
 * 3) Determine whether the PDEs are elliptic, hyperbolic, or parabolic. Show how you arrived at your classification.

Problem 2: Models of Physical Problems
The differential equation which can be written as

\cfrac{d}{dx}\left(a(x)\cfrac{du(x)}{dx}\right) + c(x)~u(x) = f(x) \qquad \text{for} \qquad x\in(O,L) ~. $$ arises in a number of fields of physics and engineering.

Write out the ODEs for at least five of these physical situations and identify the meanings of the variables.

Problem 3: Weighted Residual Methods
Give examples from the scientific literature, with full references, where the collocation, subdomain, and least squares weighted residual methods have been used.

Explain why the standard finite element method has not been used in these examples.

One example for each method will suffice.

Problem 4: Least Squares Method
The residual over an element (for the general second-order ODE) is defined as

R^e(x,c^e_1,c^e_2,\dots,c^e_n) := \cfrac{d}{dx}\left(a~\cfrac{du^e_h}{dx}\right) + c~u^e_h - f(x) $$ where

u^e_h(x) = \sum_{j=1}^n c^e_j \varphi^e_j(x) ~. $$

In the least squares method, the residual $$R^e$$ is minimized in the following sense:

\frac{\partial }{\partial c_i}\int_{x_a}^{x_b} [R^e(x,c^e_1,c^e_2,\dots,c^e_n)]^2~dx = 0, \qquad i=1,2,\dots,n ~. $$

 * 1) What is the weighting function used in the least squares approach?
 * 2) Develop a least-squares finite element model for the problem. You don't have to find the solution to a specific problem - just the components of $$\mathbf{K}$$, $$\mathbf{u}$$, and $$\mathbf{f}$$ for an element will do.
 * 3) Discuss some functions ($$\varphi^e_j$$) that can be used to approximate $$u$$.

Problem 5: Heat Transfer
The differential equation that is used to model one-dimensional heat transfer is

-\cfrac{d}{dx}\left(k~\cfrac{dT}{dx}\right) = q \qquad \text{for} \qquad x \in (0,L) $$ where $$T$$ is the temperature, $$k$$ is the thermal conductivity, and $$q$$ is the heat generation.

Let us assume that the boundary conditions are
 * $$\begin{align}

T(0) &= T_0 \\ \left.\left[ k\cfrac{dT}{dx}+\beta~(T - T_{\infty})+\hat{q}\right] \right|_{x=L} & = 0 \end{align}$$ where $$T_0$$ is a prescribed temperature, $$\beta$$ is a convective heat  transfer coefficient, $$T_{\infty}$$ is the ambient temperature, and $$\hat{q}$$ is a prescribed heat flux.

 * 1) Formulate the finite element model for this problem over one element, that is, find expressions for the stiffness matrix and the load vector.
 * 2) Use ANSYS (or any other software tool) to solve the problem of heat conduction in an insulated rod. Compare your solution with the exact solution (if available).  Try at least three refinements of the mesh to confirm that your solution converges.
 * 3) Write your own code to solve the problem in part 2 (use two elements). Compare your solution with the ANSYS (or other software) solution.

Use the following data:
 * $$\begin{align}

q         & = 0 \\ \hat{q}   & = 0 \\ L         & = 0.1~\text{m} \\ k         & = 0.01~\text{W/(m-K)} \\ \beta     & = 25~\text{W/(m}^2\text{-K)} \\ T_0       & = 323 ~\text{K} \\ T_{\infty} & = 278 ~\text{K} \end{align}$$

You can look at the 1-D heat conduction problem on the Introduction to Ansys page to see how ANSYS deals with heat conduction.

title: Nonlinear finite elements/Homework 2


text:
Nonlinear finite elements/Homework 3

Problem: The Patch Test
The patch test is one of the first checks used to determine the accuracy of a finite element.

Your first task is to read the following paper on the patch test.

Taylor, R.L., Simo, J.C., Zienkiewicz, O.C., and Chan, A.C.H, 1986, The patch test - a condition for assessing FEM convergence,  International Journal for Numerical Methods in Engineering, 22, pp. 39-62.

After reading the paper, please answer the following questions:

 * 1) What is the patch test and why is it used?
 * 2) What are the consistency requirement and the stability condition?
 * 3) What do the authors mean by convergence?
 * 4) List the 2-D structural elements available in ANSYS and describe what shape functions and degrees of freedom each uses.
 * 5) Check whether each of these elements passes the patch test. Show details of your calculations.
 * 6) Which parts of the paper did you not understand?

title: Nonlinear finite elements/Homework 3


text:
Nonlinear finite elements/Homework 4

Problem 1: Axially Loaded Bar with Nonlinear Modulus
In the learning resource section we have discussed the solution of the problem of a nonlinear axially loaded bar with a body load. A Matlab code for solving the problem with finite elements can also be found on the relevant page.

In the code, we apply a tensile load at the end of the bar and use linear shape functions. Therefore, each element has two nodes.

In this problem, you are required to modify the Matlab code to use quadratic elements. Each element will have three nodes ($$1$$, $$2$$, and $$3$$). Assume that all elements have the same length $$h$$ and that node $$2$$ is at the middle of the element. The shape functions for the nodes are

\begin{align} N_1(x) & = \cfrac{2(x_2 - x)(x_3 - x)}{h^2} \\ N_2(x) & = \cfrac{4(x - x_1)(x_3 - x)}{h^2} \\ N_3(x) & = \cfrac{2(x - x_1)(x - x_2)}{h^2} \end{align} $$ Apply a compressive load $$R$$ at the free end of the bar.

 * 1) Plot the shape functions for an element. Show that they add up to 1 at all points within the element.
 * 2) Plot the displacement, stress, strain, and the error norm when you use 10 elements to discretize the bar. See the figure in the learning resource section for how you should plot these data.  Discuss and comment on your observations.
 * 3) Do the same for 100 elements. What are your observations on the effect of mesh refinement on the solution?
 * Now, increase the applied load. At what load does the solution fail to converge?Explain why the solution fails to converge at this load.

You could alternatively write your own Maple code or use ANSYS for these calculations depending on your interests.

Problem 2: Nonlinear Thermal Conduction
Consider the model chip system shown in Figure 1. Assume that heat is being generated in the silicon carbide region of the chip at the rate of 100 W/cm$$^2$$. The cooling channels contain a fluid at a constant 250 K and (supposedly) can dissipate heat at the rate of 800 W/cm$$2$$. The surroundings are at room temperature (300 K).

The thermal conductivity of the silicon carbide and the aluminum nitride are also given in Figure 1. The thermal conductivity of the epoxy resin is 1.5 W/m-K.

 * 1) First run a linear simulation of the problem using room temperature thermal conductivities of the materials.  This will act as a check of the model and setup and will help you determine any additional constraints that you will need.  Do the cooling channels really dissipate heat at 800 W/cm$$2$$?
 * 2) Fit curves to the thermal conductivity versus temperature data for Silicon Carbide and Aluminum Nitride.
 * 3) Perform a nonlinear simulation using these temperature dependent properties to find the equilibrium temperature inside the silicon carbide region. State all simplifying assumptions that you have made.

title: Nonlinear finite elements/Homework 4


text:
Nonlinear finite elements/Homework 5

Problem 1: Nonlinear Beam Bending
The differential equations governing the bending of straight beams are

{ \begin{align} \cfrac{dN_{xx}}{dx} + f(x) & = 0 \\ \cfrac{dV}{dx} + q(x) & = 0 \\ \cfrac{dM_{xx}}{dx} -V + N_{xx}\cfrac{dw_0}{dx} & = 0 \end{align} } $$ 1) Show that the weak forms of these equations can be written as
 * $$\begin{align}

\int_{x_a}^{x_b} \cfrac{dv_1}{dx}N_{xx}~dx & = \int_{x_a}^{x_b} v_1 f~dx + \left.v_1 N_{xx}\right|_{x_a}^{x_b}\\ \int_{x_a}^{x_b} \left[ \cfrac{dv_2}{dx}\left(\cfrac{dw_0}{dx}N_{xx}\right) - \cfrac{d^2v_2}{dx^2}M_{xx}\right]~dx & = \int_{x_a}^{x_b} v_2 q~dx + \left.v_2\left(\cfrac{dw_0}{dx}N_{xx}+\cfrac{dM_{xx}}{dx}\right) \right|_{x_a}^{x_b} - \left.\cfrac{dv_2}{dx}M_{xx} \right|_{x_a}^{x_b} \end{align}$$

2) The  von Karman strains are related to the displacements by
 * $$\begin{align}

\varepsilon_{xx} & = \varepsilon^0_{xx} + z\varepsilon^1_{xx} \\ \varepsilon^0_{xx} & = \cfrac{du_0}{dx} + \frac{1}{2}\left(\cfrac{dw_0}{dx}\right)^2 \\ \varepsilon^1_{xx} & = -\cfrac{d^2w_0}{dx^2} \end{align}$$ The stress and moment resultants are defined as
 * $$\begin{align}

N_{xx} & = \int_A \sigma_{xx}~ dA \\ M_{xx} & = \int_A z\sigma_{xx}~ dA \end{align}$$ For a linear elastic material, the stiffnesses of the beam in extension and bending are defined as
 * $$\begin{align}

A_{xx} & = \int_A E~dA\qquad \leftarrow \qquad \text{extensional stiffness}\\ B_{xx} & = \int_A zE~dA\qquad \leftarrow \qquad \text{extensional-bending stiffness}\\ D_{xx} & = \int_A z^2E~dA\qquad \leftarrow \qquad \text{bending stiffness} \end{align}$$ where $$E$$ is the Young's modulus of the material.

Derive expressions for $$N_{xx}$$ and $$M_{xx}$$ in terms of the displacements $$u_0$$ and $$w_0$$ and the extensional and bending stiffnesses of the beam assuming a linear elastic material.

3) Express the weak forms in terms of the displacements and the extensional and bending stiffnesses.

4) Assume that the approximate solutions for the axial displacement and the transverse deflection over a two noded element are given by
 * $$\begin{align}

u_0(x) & =u_1 \psi_1(x) + u_2 \psi_2(x) \\ w_0(x) & =w_1 \phi_1(x) + \theta_1 \phi_2(x) + w_2 \phi_3(x) + \theta_2 \phi_4(x) \end{align}$$ where $$\theta = -(dw_0/dx)$$.

Compute the element stiffness matrix for the element.

5) Show the alternate procedure by which the element stiffness matrix can be made symmetric.

6)Derive the element tangent stiffness matrix for the element.

Problem 2: Numerical simulation
Consider the cantilever beam shown in Figure 1. The beam is made of 6061-T6 aluminum alloy.The Young's modulus of the beam is 69 GPa, the Poisson's ratio is 0.33, the tensile yield strength is 275 MPa.

The length of the beam is 5 m, the width and the thickness are both 10 cm. A point load of 10 kN is applied at the free end as shown in the figure.

 * 1) Perform a linear three-dimensional finite element analysis of the beam using 20 elements along the length and 3 elements across the thickness. Then,
 * 2) Plot the axial stress of the centroidal line as a function of distance from the fixed end.
 * 3) Plot the axial displacement of the centroidal as a function of distance from the fixed end.
 * 4) Plot the transverse displacement of the centroidal line as a function of distance from the fixed end.
 * 5) Perform a  nonlinear three-dimensional finite element analysis of the beam using 20 elements along the length and 3 elements across the thickness. Then,
 * 6) Plot the axial stress of the centroidal line as a function of distance from the fixed end.
 * 7) Plot the axial displacement of the centroidal as a function of distance from the fixed end.
 * 8) Plot the transverse displacement of the centroidal line as a function of distance from the fixed end.
 * 9) Perform a linear finite element analysis of the beam using 4 beam elements. Then,
 * 10) Plot the axial stress of the centroidal line as a function of distance from the fixed end.
 * 11) Plot the axial displacement of the centroidal as a function of distance from the fixed end.
 * 12) Plot the transverse displacement of the centroidal line as a function of distance from the fixed end.
 * 13) Perform a  nonlinear finite element analysis of the beam using 4 beam elements. Then,
 * 14) Plot the axial stress of the centroidal line as a function of distance from the fixed end.
 * 15) Plot the axial displacement of the centroidal as a function of distance from the fixed end.
 * 16) Plot the transverse displacement of the centroidal lineas a function of distance from the fixed end.
 * 17) Comment on your results.

title: Nonlinear finite elements/Homework 5


text:
Nonlinear finite elements/Homework 6

Problem 1: Continuum-based Nonlinear Beam Bending
You have seen that the derivation of the finite element equations for beams can get quite complicated.

One alternative is to use the continuum-based approach given in Belytschko, Liu, and Moran. Your task for this problem is to read the portion on beams from the chapter on Beams and Shells from that book. After you have completed that task, answer the following questions.

 * 1) What is the continuum-based approach for the finite element analysis of beams? Draw figures to elucidate.
 * 2) Why is the continuum-based approach preferred in nonlinear finite element analysis?
 * 3) The rate of deformation is defined as

\boldsymbol{D} = \text{sym}(\boldsymbol{\nabla} \mathbf{v}) = \frac{1}{2}\left[\boldsymbol{\nabla} \mathbf{v} + (\boldsymbol{\nabla} \mathbf{v})^T\right] $$ where $$\mathbf{v}$$ is the velocity.In index notation, we write

D_{ij} = \text{sym}(v_{i,j}) = \frac{1}{2}\left(v_{i,j} + v_{j,i}\right) = \frac{1}{2}\left(\frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right), \qquad i,j=1,2,3~. $$ Given the above definition, derive equations (9.2.1) through (9.2.7) from the book.
 * 4) Why is the Euler-Bernoulli beam theory $$C^1$$? Why is the Timoshenko beam theory $$C^0$$?
 * 5) Why is a $$C^0$$ beam theory preferable?
 * 6) What are the major assumptions in the Timoshenko beam formulation described in the handout (the authors call these the modified Mindlin-Reissner assumptions)?
 * 7) Derive the equation (9.3.13) of the book starting from equation (9.3.3).

Problem 1a: Curved laminated beam
To enforce the plane stress assumption, we have to rotate the rate of deformation so that its components are aligned with a local coordinate system aligned with the laminae. This exercise is meant to improve your understanding of the process.

Consider the curved composite beam shown in Figure 1.

Assume that the beam has been shaped into an arc of a circle. The material of the beam is a transversely isotropic fiber composite material with the fibers running along the length of the beam. The rate constitutive relation of the material is given by

\cfrac{D}{Dt} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{31} \\ \sigma_{12} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{13} & 0 & 0 & 0 \\ C_{13} & C_{13} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & (C_{11}-C_{12})/2 \end{bmatrix} \begin{bmatrix} D_{11} \\ D_{22} \\ D_{33} \\ D_{23} \\ D_{31} \\ D_{12} \end{bmatrix}~. $$

 * 1) Construct a laminar coordinate system at the blue point.  (Use equations 9.3.16 and 9.3.17 from the book).  Assume that the blue point is at the center of element 5.
 * 2) Assume that the velocities at nodes 5 and 6 are $$\mathbf{v}_5$$ and $$\mathbf{v}_6$$, respectively.  Compute the velocity gradient at the blue point.
 * 3) Compute the rate of deformation at the blue point.
 * 4) Rotate the rate of deformation so that its components are with respect to the laminar coordinate system.
 * 5) Compute the stress rate at the blue point.  This stress rate is expressed in terms of the laminar coordinate system.
 * 6) Apply the plane stress condition and solve for $$D_{yy}$$.
 * 7) Express the stress rate and the modified rate of deformation in the global coordinate system.

You may use numbers instead of symbols to make the above process more concrete in your mind.

Problem 1b: ( Optional: Extra Credit)
Show how we obtain equation 9.3.22 of the book. Recall that, in matrix notation, the nodal internal force is given by

\mathbf{f}^{\text{int}} = \int_{\Omega} \mathbf{B}^T \boldsymbol{\sigma}~d\Omega $$ where the $$\mathbf{B}$$ matrix relates the rate of deformation to the nodal velocities ($$\mathbf{v}$$) by

\mathbf{D} = \mathbf{B} \mathbf{v} ~. $$

Problem 1c: ( Optional: Extra Credit)
As you have noticed, not all the steps on the finite element solution process have been described by the authors. Write down the extra steps that you think you will need to complete the solution process.

Problem 2: Exploring commercial software
Track down the different types of beam and shell elements that ANSYS and LS-DYNA provide.For each element that you find, answer the following
 * 1) What is the name of the element?
 * 2) How many nodes does the element have?
 * 3) How many displacement type (displacements, rotations) degrees of freedom are there at each node?  List them.
 * 4) How many force type (forces, moments) boundary conditions can be applied at each node?List them.
 * 5) What is the form of the shape functions used by the element?  (Alternatively, you may write down the approximate solutions assumed by the element.)
 * 6) Is the beam element a continuum-based one or one derived from the equations of beam theory?
 * 7) Does the element allow for full and reduced integration?

title: Nonlinear finite elements/Homework 6


text:
Nonlinear finite elements/Homework 7

Problem 1: Index Notation

 * Determine whether the following expressions are valid in index notation. If valid, identify the free indices and the dummy indices.
 * $$\begin{align}

A_{ms} & = b_m (c_r - d_r) \\ A_{ms} & = b_m (c_s - d_s) \\ t_i& = \sigma_{ji} n_j \\ t_i& = \sigma_{ji} n_i \\ x_i x_i & = r^3 \\ B_{ij} c_j & = 3 \end{align}$$
 * Show the following
 * $$\begin{align}

\delta_{ii} &= 3 \\ e_{ijk}e_{pqk} &= \delta_{ip}\delta_{jq} - \delta_{iq}\delta_{jp} \\ \delta_{ij}e_{ijk} & = 0 \\ e_{qrs} d_q d_s & = 0 \end{align}$$
 * The elasticity tensor is given by

\boldsymbol{\mathsf{C}} = \lambda~ \boldsymbol{\mathit{1}}\otimes\boldsymbol{\mathit{1}} + 2\mu~\boldsymbol{\mathsf{I}} $$ where $$\lambda, \mu$$ are Lame constants, $$\boldsymbol{\mathit{1}}$$ is the second order identity tensor, and $$\boldsymbol{\mathsf{I}}$$ is the fourth-order symmetric identity tensor.The two identity tensors are defined as
 * $$\begin{align}

\boldsymbol{\mathit{1}} & = \delta_{ij}~\mathbf{e}_i\otimes\mathbf{e}_j \\ \boldsymbol{\mathsf{I}} & = \frac{1}{2}[\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}]~ \mathbf{e}_i\otimes\mathbf{e}_j\otimes\mathbf{e}_k\otimes\mathbf{e}_l \end{align}$$ The stress-strain relation is

\boldsymbol{\sigma} = \boldsymbol{\mathsf{C}} : \boldsymbol{\varepsilon} $$ Show that the stress-strain relation can be written in index notation as

\sigma_{ij} = 2\mu\varepsilon_{ij} + \lambda\varepsilon_{kk}\delta_{ij}~. $$ Write the stress-strain relations in expanded form.

Problem 2: Rotating Vectors and Tensors
Let ($$\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3$$) be an orthonormal basis.Let $$\boldsymbol{A}$$ be a second order tensor and $$\mathbf{u}$$ be a vector with components
 * $$\begin{align}

\boldsymbol{A} &= 5~\mathbf{e}_1\otimes\mathbf{e}_1 - 4~\mathbf{e}_2\otimes\mathbf{e}_1 + 2~\mathbf{e}_3\otimes\mathbf{e}_3\\ \mathbf{u} &= -2~\mathbf{e}_1 + 3~\mathbf{e}_3 \end{align}$$

the basis ($$\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3$$).
 * Write out $$\boldsymbol{A}$$ and $$\mathbf{u}$$ in matrix notation.
 * Find the components of the vector $$\mathbf{v} = \boldsymbol{A}\bullet\mathbf{u}$$ in the basis ($$\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3$$).
 * Find the components of the vector $$\mathbf{w} = \mathbf{v}\times\mathbf{u}$$ in
 * Find the components of the tensor $$\boldsymbol{C} = \mathbf{v}\otimes\mathbf{w}$$ in the orthonormal basis.
 * Rotate the basis clockwise by 30 degrees around the $$\mathbf{e}_3$$ direction. Find the components of $$\mathbf{u}$$, $$\mathbf{v}$$, $$\mathbf{w}$$, and $$\boldsymbol{C}$$ in the rotated basis.

Problem 3: More Beams
Use ANSYS to solve the following problems.We want to see how the beam element BEAM188 that ANSYS provides behaves. If you use some other tool choose the equivalent element.

Beam188 is based on Timoshenko beam theory. The element uses linear shape functions for all degrees of freedom.The beam element refers to the following papers:
 * Simo, J.C. and Vu-Quoc, L., 1986, "A three-dimensional finite strain rod model. Part II: Computational Aspects," Computer Methods in Applied Mechanics and Engineering,  58, pp. 79-116.
 * Ibrahimbegovic, A., 1995, "On finite element implementation of geometrically nonlinear Reissner's beam theory: Three-dimensional curved beam elements," Computer Methods in Applied Mechanics and Engineering,  122, pp. 11-26.

Part A:
Consider a beam of length $$L$$ = 100 in., cross-section 1 in. $$\times$$ 1 in., and subjected to a uniformly distributed transverse load $$q_0$$ lbf/in. Model one half of the beam using symmetry considerations.

Hinged-Hinged Beam
The boundary conditions are

w_0(0) = u_0(L/2) = \varphi_x(L/2) = 0 ~. $$ Compute a plot for this case using Beam188 elements. What do you observe?

Clamped-Clamped Beam
The boundary conditions are

u_0(0) = w_0(0) = \varphi_x(0) = u_0(L/2) = \varphi_{x = L/2} = 0. $$ Compute a plot for this case using Beam188 elements. Comment on your plot.

You will have the save your results at the end of each load step to get the data you need.

Part B:
For this problem you will try to reproduce some of the results given in Simo and Vu Quoc (1986) and Ibrahimgbegovic (1995) using Beam188 elements.

 * 1) Simulate the unrolling of a cantilever beam from Section 4.1.1 of Ibrahimbegovic (1995) and compare your results with the results shown in the paper.
 * 2) Simulate the clamped-hinged deep circular arch from Example 7.3 of Simo and Vu Quoc (1986) and compare you results with the results shown in the paper.
 * 3) Simulate the buckling of a hinged right-angle frame under both fixed and follower loads from Example 7.4 of Simo and Vu Quoc (1986) and compare your results with those shown in the paper.

title: Nonlinear finite elements/Homework 7


text:
Nonlinear finite elements/Homework 8

Problem 1: Plate Bending
Consider the flat plate shown in Figure 1.

Transverse normals are defined as the straight lines perpendicular to the midsurface. If we assume that the plate is thin, the Kirchhoff-Love hypothesis holds, i.e.,
 * Transverse normals remain straight.
 * Transverse normals remain normal to the midsurface.
 * Transverse normals are inextensible.

Answer the following questions.

Part 1
The displacement field that satisfies the Kirchhoff-Love hypothesis is given by
 * $$\begin{align}

u_1(x_1,x_2,x_3) & = u^0_1(x_1,x_2) - x_3~u^0_{3,1} \\ u_2(x_1,x_2,x_3) & = u^0_2(x_1,x_2) - x_3~u^0_{3,2} \\ u_3(x_1,x_2,x_3) & = u^0_3(x_1,x_2) \end{align}$$ where $$u^0_i$$ are the displacements of the midsurface.

Write the above equations (using the standard partial derivative notation with $$\partial$$~s) after making the following substitutions.
 * $$\text{(2)} \qquad

\begin{align} & x := x_1 ~; y := x_2 ~; z := x_3 \\ & u := u_1 ~; v := u_2 ~; w := u_3 \\ & u_0 := u^0_1 ~; v_0 := u^0_2 ~; w_0 := u^0_3 ~. \end{align} $$

Part 2
Assume that the plate undergoes infinitesimal strains and rotations. Then the strain-displacement relations are

\boldsymbol{\varepsilon} = \frac{1}{2}\left[\boldsymbol{\nabla} \mathbf{u} + (\boldsymbol{\nabla} \mathbf{u})^T\right]~. $$ Show that, in terms of the midsurface displacements, the strain-displacement relations can be written in index notation as
 * $$\begin{align}

\varepsilon_{\alpha\beta} & = \frac{1}{2} \left(u^0_{\alpha,\beta} + u^0_{\beta,\alpha}\right) - x_3~u^0_{3,\alpha\beta} \\ \varepsilon_{\alpha 3} & = 0 \\ \varepsilon_{33} & = 0 \end{align}$$ where $$\alpha = 1,2$$ and $$\beta = 1,2$$.

Part 3
The internal virtual work in the plate is given by (using index notation)

\delta W_{\text{int}} = \int_{\Omega_0}\left(\int_{-h}^{h} \sigma_{\alpha\beta}~\delta \varepsilon_{\alpha\beta}~dx_3\right) d\Omega_0 ~. $$ Note that if a component of the strain tensor is zero, then the variation of that strain component is also zero.

Show that the internal virtual work can be written as
 * $$\text{(3)} \qquad

\delta W_{\text{int}} = \int_{\Omega_0}\left( N_{\alpha\beta}\delta \varepsilon_{\alpha\beta}^0 + M_{\alpha\beta}\delta \varepsilon_{\alpha\beta}^1\right)~d\Omega_0 $$ where
 * $$\begin{align}

N_{\alpha\beta} & = \int_{-h}^h \sigma_{\alpha\beta}~dx_3 \\ M_{\alpha\beta} & = \int_{-h}^h x_3~\sigma_{\alpha\beta}~dx_3 \\ \varepsilon_{\alpha\beta}^0 & = \frac{1}{2} \left(u^0_{\alpha,\beta} + u^0_{\beta,\alpha}\right) \\ \varepsilon_{\alpha\beta}^1 & = - u^0_{3,\alpha\beta} ~. \end{align}$$ Write down the expanded form of equation (3).

Part 4
Reddy's book uses the notation of equation (2) to derive the following  equilibrium equations for the plate.
 * $$\begin{align}

\frac{\partial N_{xx}}{\partial x} + \frac{\partial N_{xy}}{\partial y} & = 0 \\ \frac{\partial N_{xy}}{\partial x} + \frac{\partial N_{yy}}{\partial y} & = 0 \\ \frac{\partial^2 M_{xx}}{\partial x^2} + 2\frac{\partial^2 M_{xy}}{\partial x\partial y} + \frac{\partial^2 M_{yy}}{\partial y^2} + \mathcal{N}(u_0,v_0,w_0) + q & = 0 \end{align}$$ where

\mathcal{N}(u_0,v_0,w_0) = \frac{\partial }{\partial x} \left(N_{xx}\frac{\partial w_0}{\partial x} + N_{xy}\frac{\partial w_0}{\partial y}\right) + \frac{\partial }{\partial y}\left(N_{xy}\frac{\partial w_0}{\partial x} + N_{yy}\frac{\partial w_0}{\partial y} \right) ~. $$ Write the above equations in index notation. Then, show that in index-free  Gibbs notation, the equilibrium equations for the plate can be written as
 * $$\begin{align}

\boldsymbol{\nabla} \bullet \boldsymbol{N} & = \mathbf{0} \text{(4)} \qquad \\ \boldsymbol{\nabla} \bullet (\boldsymbol{\nabla}\bullet\boldsymbol{M}) + \boldsymbol{\nabla} \bullet (\boldsymbol{N}\bullet\boldsymbol{\nabla u^0_3)} + q & = 0 \text{(5)} \qquad \end{align}$$ where

\begin{align} \boldsymbol{N} & = N_{\alpha\beta}~\mathbf{e}_{\alpha}\otimes\mathbf{e}_{\beta} ~;& \boldsymbol{M} & = M_{\alpha\beta}~\mathbf{e}_{\alpha}\otimes\mathbf{e}_{\beta}\\ \boldsymbol{\nabla} \varphi & = \varphi_{,\alpha} ~\mathbf{e}_{\alpha} ~;& \boldsymbol{\nabla} \mathbf{v} & = v_{\alpha,\beta} ~\mathbf{e}_{\alpha}\otimes\mathbf{e}_{\beta} ~;& \boldsymbol{\nabla} \boldsymbol{S} & = S_{\alpha\beta,\gamma} ~\mathbf{e}_{\alpha}\otimes\mathbf{e}_{\beta}\otimes\mathbf{e}_{\gamma} \\ \boldsymbol{\nabla} \bullet \mathbf{v} & = v_{\alpha,\alpha} ~;& \boldsymbol{\nabla} \bullet \boldsymbol{S} & = S_{\alpha\beta,\beta}~\mathbf{e}_{\alpha} ~;& \boldsymbol{\nabla} \bullet \boldsymbol{T} & = T_{\alpha\beta\gamma,\gamma} ~\mathbf{e}_{\alpha}\otimes\mathbf{e}_{\beta} \end{align} $$ and $$\alpha=1,2$$, $$\beta=1,2$$, and $$\gamma=1,2$$.$$\varphi$$ is a scalar, $$\mathbf{v}$$ is a vector, $$\boldsymbol{S}$$ is a second-order tensor, and $$\boldsymbol{T}$$ is a third-order tensor.

Part 5
Derive the symmetric weak forms of the equilibrium equations starting with equations (4) and (5). Use the coordinate-free Gibbs notation in your derivation.

You will find the following identities useful in your derivation
 * $$\begin{align}

\boldsymbol{\nabla} \bullet (\varphi\mathbf{v}) & = \varphi (\boldsymbol{\nabla} \bullet \mathbf{v}) + \mathbf{v}\bullet(\boldsymbol{\nabla} \varphi)\\ \boldsymbol{\nabla} \bullet (\boldsymbol{S}^T\bullet\mathbf{v}) & = \mathbf{v}\bullet(\boldsymbol{\nabla} \bullet \boldsymbol{S}) + \boldsymbol{S}:(\boldsymbol{\nabla} \mathbf{v}) \end{align}$$

Part 6

 * Write down the natural boundary conditions in Gibbs notation.
 * What are the essential boundary conditions that can be applied to the plate?

Problem 2: Transient Heat Conduction
Consider the model problem shown in Figure 2.

A 4340 steel casing is used to contain a solid high energy (HE) material. The HE material reacts on the inner (star-shaped) surface and generates gases.Since the gases are confined, the temperature of the gases increases according to the ideal gas law ($$PV = nRT$$).

We will try to simulate the transient heat transfer and the associated thermal stresses in the HE material and the steel casing.

Since the actual problem is quite complex and involves chemical reactions, equations of state, and fluid-structure interaction, we will simplify the problem.

The assumptions we will make are:

 * 1) The burn-surface does not change position.
 * 2) The temperature is the same at all points on the burn surface.
 * 3) The temperature at the burn surface increases exponentially with time according to the relation

T = T_0 + A~\exp{(B~t)} $$
 * where $$T_0$$ is the initial temperature (300 K), $$T$$ is the temperature at time $$t$$, $$A = 2$$ K, and $$B = 3$$ K/s.

The inner radius of the steel casing is 5 cm. The outer radius of the steel casing is 6 cm.

The inner and outer radii of the star shaped region are 4 cm and 2 cm, respectively. You can make your own assumptions about the shape of the curves in the star shaped region.

Use standard values for the density, Young's modulus, Poisson's ratio, coefficient of thermal expansion, and thermal conductivity for the steel. You can use Matweb to find the properties you need.

Assume that the high energy material has the following properties:
 * Density = 1200 kg/m$$^3$$~.
 * Young's modulus = 1 GPa~.
 * Poisson's ratio = 0.49~.
 * Coefficient of thermal expansion = 1.2$$\times$$10$$^{-4}$$ /K~.
 * Thermal Conductivity = 5 W/m-K (this is 10 times higher than theactual value).
 * Specific Heat Capacity = 1000 J/kg-K.

 * 1) Assume that the initial temperature of the steel and the HE material is 300 K.Also assume that there is no heat flux on the external surface of the steel casing. Plot the temperature contours and stresses in the two materials at $$t = 1$$ s and $$t = 2$$ s.
 * 2) Assume that the initial temperature of the steel and the HE material is 300 K.Also assume that the external surface of the steel casing is maintained at 300 K. Plot the temperature contours and stresses in the two materials at $$t = 1$$ s and $$t = 2$$ s.

You should look at examples given in VM 33 and VM 174 in the ANSYS Verification manual for ideas on how to solve this problem.

Problem 3: Explicit Finite Elements
Figure 3(a) shows a beam that is built-in at one and loaded at the other. Figure 3(b) shows a plate that is built-in at one end and loaded at the other. The load that is applied is transient and has the form shown in Figure 3(c).

Assume that the beam and the plate are made of 4340 steel. The beam has dimensions 12 in. $$\times$$ 0.25 in. $$\times$$ 0.25 in. The plate has dimensions 12 in. $$\times$$ 0.25 in. $$\times$$ 6 in.

Simulate the problem for 10 secs. using LS-DYNA or some other explicit FE code. Use beam elements for the beam and shell elements for the plate.Use an element size of 2 in. for the beam and 2 in. $$\times$$ 2 in. for the shell.

 * 1) Plot the displacement of the end of the beam as a function of time.
 * 2) Plot the axial stress at two points on the top of the beam as a function of time.
 * 3) Plot the displacement of the end of the plate as a function of time.
 * 4) Plot the axial stress at two points on the top of the plate as a function of time.

After you run the LS-DYNA, the program will create a binary file 'd3plot' and ASCII file 'nodout'. To view the result:

 * 1) Run LS-Prepost and open the binary file 'd3plot'.
 * 2) To view the deflection at free end: click ASCII button.
 * 3) In the window below it, select 'nodout', then click LOAD.
 * 4) Available nodes will be shown in the next window below.
 * 5) Select the node.  For shell element, those are 2 nodes at the free end.
 * 6) Select 'y-displacement' for beam, and 'z-displacement' for shell, then click PLOT.

You should also try the animation tool to get a feel for how the beam and shell deform over time.

title: Nonlinear finite elements/Homework 8


text:
Nonlinear finite elements/Homework 9

Problem 1: Total Lagrangian
Consider the tapered two-node element shown in Figure 1. The displacement field in the element is linear.

The reference (initial) cross-sectional area is

A_0 = (1-\xi)~ A_{01} + \xi~ A_{02} ~. $$ Assume that the nominal (engineering) stress is also linear in the element, i.e.,

P = (1-\xi)~ P_1 + \xi~ P_2 ~. $$

 * 1) Using the total Lagrangian formulation, develop expressions for the internal nodal forces.
 * 2) What are the internal nodal forces if the reference area and the nominal stress are constant over the element?
 * 3) Assume that the body force is constant. Develop expressions for the external nodal forces for that case.
 * 4) What are the external nodal forces if the reference area and the nominal stress are constant over the element?
 * 5) Develop an expression for the consistent mass matrix for the element.
 * 6) Obtain the lumped (diagonal) mass matrix using the row-sum technqiue.
 * 7) Obtain the lumped (diagonal) mass matrix using the row-sum technqiue.
 * 8) Find the natural frequencies of a single element with consistent mass by solving the eigenvalue problem

\mathbf{K}~\mathbf{u} = \omega^2~\mathbf{M}~\mathbf{u} $$
 * with

\mathbf{K} = \cfrac{E(A_{01} + A_{02}}{2 l_0} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} $$
 * where $$E$$ is the Young's modulus and $$l_0$$ is the initial length of the element.

Problem 2: Updated Lagrangian
Consider the tapered two-node element shown in Figure 1.

The current cross-sectional area is

A = (1-\xi)~ A_1 + \xi~ A_2 ~. $$ Assume that the Cauchy stress is also linear in the element, i.e.,

\sigma = (1-\xi)~ \sigma_1 + \xi~ \sigma_2 ~. $$

 * 1) Using the updated Lagrangian formulation, develop expressions for the internal nodal forces.
 * 2) Assume that the body force is constant.Develop expressions for the external nodal forces for that case.

Problem 3: Modal Analysis

 * 1) Consider the axially loaded bar in problem VM 59 of the ANSYS Verification manual. Assume that the bar is made of Tungsten carbide.
 * 2) Find the fundamental natural frequency of the bar.
 * 3) Find the first three modal frequencies for a load of 40,000 lbf.
 * 4) Consider the stretched circular membrane in problem VM 55 of the ANSYS Verification manual. Assume that the membrane is made of OFHC (Oxygen-free High Conductivity) copper.
 * 5) Find the fundamental natural frequency of the bar.
 * 6) Find the first five modal frequencies for a load of 10,000 lbf.

title: Nonlinear finite elements/Homework 9


text:
Nonlinear finite elements/Homework 10

Problem 1: Kinematics and Stress Rates
Figure 1 shows a linear three-noded triangular element in the reference configuration.

The motion of the nodes is given by:
 * $$\begin{align}

\text{Node}~1:&& x_1(t) & = 0& \qquad y_1(t) & = 0 \\ \text{Node}~2:&& x_2(t) & = 2(1 + at)\cos\cfrac{\pi t}{2} &\qquad y_2(t) & = 2(1 + at)\sin\cfrac{\pi t}{2} \\ \text{Node}~3:&& x_3(t) & = -(1 + bt)\sin\cfrac{\pi t}{2} &\qquad y_3(t) & =(1 + bt)\cos\cfrac{\pi t}{2} \end{align}$$ The configuration ($$x,y$$) of the element at time $$t$$ is given by

x(\boldsymbol{\xi}, t) = \sum_{i=1}^3 x_i(t)~N_i(\boldsymbol{\xi}) ~;\qquad y(\boldsymbol{\xi}, t) = \sum_{i=1}^3 y_i(t)~N_i(\boldsymbol{\xi}) $$

 * 1) Write down expressions for $$N_1$$, $$N_2$$, and $$N_3$$ in terms of the initial configuration ($$X,Y$$) ?
 * 2) Derive expressions for the deformation gradient and the Jacobian determinant for the element as functions of time.
 * 3) What are the values of $$a$$ and $$b$$ for which the motion is isochoric?
 * 4) For which values of $$a$$ and $$b$$ do we get invalid motions?
 * 5) Derive the expression for the Green (Lagrangian) strain tensor for the element as a function of time.
 * 6) Derive an expression for the velocity gradient as a function of time.
 * 7) Compute the rate of deformation tensor and the spin tensor.
 * 8) Assume that $$a=1$$ and $$b=2$$. Sketch the undeformed configuration and the deformed configuration at $$t=1$$ and $$t=2$$. Draw both the deformed and undeformed configurations on the same plot and label.
 * 9) Compute the polar decomposition of the deformation gradient with the above values of $$a$$ and $$b$$.
 * 10) Assume an isotropic, hypoelastic constitutive equation for the material of the element. Compute the material time derivative of the Cauchy stress using (a) the Jaumann rate and (b) the Truesdell rate.

Problem 2: Hyperelastic Pinched Cylinder Problem
Read the following paper on shells:

Buchter, N., Ramm, E., and Roehl, D., 1994, "Three-dimensional extension of non-linear shell formulation based on the enhanced assumed strain concept," '' Int. J. Numer. Meth. Engng.'',  37, pp. 2551-2568.

Answer the following questions:

 * 1) What do the authors mean by "enhanced assumed strain" ?
 * 2) Example 8.2 (and Figures 3 and 4 and Table III) in the paper discusses the simulation of a hyperelastic cylinder. Perform the same simulation using ANSYS for a shell thickness of 0.2 cm.  Use shell elements and the Neo-Hookean hyperelastic material model that ANSYS provides.
 * 3) Compare the total load needed to achieve an edge displacement of 16 cm with the results given in Table III. Comment on your results.

Problem 3: Elastic-Plastic Punch Indentation
Read the following paper on elastic-viscoplastic FEA:

Rouainia, M. and Peric, D., 1998, "A computational model for elasto-viscoplastic solids at finite strain with reference to thin shell applications," '' Int. J. Numer. Meth. Engng.'',  42, pp. 289-311.

Answer the following questions:

 * 1) Example 5.4 of the paper shows a simulation of the deformation of a thin sheet by a square punch. Perform the same simulation for 6061-T6 aluminum. Assume linear isotropic hardening and no rate dependence.
 * 2) Draw a plot of the punch force vs. punch travel and compare your result with the results shown in Figure 13 of the paper (qualitative comparison only).

title: Nonlinear finite elements/Homework 10


text:
Nonlinear finite elements/Homework 11

Problem 1: Small Strain Elastic-Plastic Behavior
For small strains, the strain tensor is given by

\boldsymbol{\varepsilon} = \frac{1}{2}\left[\boldsymbol{\nabla} \mathbf{u} + (\boldsymbol{\nabla} \mathbf{u})^T\right] \qquad\text{or}\qquad \varepsilon_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}) ~. $$ In classical (small strain) rate-independent plasticity we start off with an additive decomposition of the strain tensor

\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^e + \boldsymbol{\varepsilon}^p \qquad\text{or}\qquad \varepsilon_{ij} = \varepsilon_{ij}^e + \varepsilon_{ij}^p ~. $$ Assuming linear elasticity, we have the following elastic stress-strain law

\boldsymbol{\sigma} = \boldsymbol{\mathsf{C}} : \boldsymbol{\varepsilon}^e \qquad\text{or}\qquad \sigma_{ij} = C_{ijkl}\varepsilon_{kl}^e ~. $$ Let us assume that the $$J_2$$ theory applies during plastic deformation of the material. Hence, the material obeys an associated flow rule
 * $$\text{(1)} \qquad

\dot{\boldsymbol{\varepsilon}}^p = \dot{\gamma}\frac{\partial f(\boldsymbol{\sigma},\alpha,T)}{\partial \boldsymbol{\sigma}} $$ where $$\dot{\gamma}$$ is the plastic flow rate, $$f$$ is the yield function, $$T$$ is the temperature, and $$\alpha$$ is an internal variable.

Answer the following questions. Show your derivations in a clear and step-by-step manner.

Part 1
Let $$\alpha$$ be the equivalent plastic strain, defined as

\alpha := \sqrt{\cfrac{2}{3}} \lVert\boldsymbol{\varepsilon}^p\rVert \qquad \text{where} \qquad \lVert\mathbf{a}\rVert = \sqrt{\mathbf{a}:\mathbf{a}} ~. $$ Express the time derivative of $$\alpha$$ in terms of $$\dot{\gamma}$$ and $$\partial f/\partial \boldsymbol{\sigma}$$. This is the evolution law for $$\alpha$$.

Part 2
For an adiabatic process, the rate of change of temperature can be written as

\dot{T} = \cfrac{\chi}{\rho C_p} \boldsymbol{\sigma}:\dot{\boldsymbol{\varepsilon}}^p $$ where $$\chi$$ is the Taylor-Quinney coefficient, $$\rho$$ is the density, and $$C_p$$ is the specific heat. Express $$\dot{T}$$ in terms of $$\dot{\gamma}$$ and $$\partial f/\partial \boldsymbol{\sigma}$$. This is the evolution law for $$T$$.

Part 3
Write down the rate form of the elastic stress-strain law. Assume that deformations are small so that objectivity of the rates is not a concern.

Part 4
The consistency condition during plastic flow requires that

\dot{f}(\boldsymbol{\sigma}, \alpha, T) = 0 ~. $$ Write down an expression for the time derivative of $$f(\boldsymbol{\sigma},\alpha,T)$$ using the chain rule.

Part 5
Use the consistency condition and the expressions you have derived in the previous parts to derive an expression for $$\dot{\gamma}$$ in terms of $$\partial f/\partial\boldsymbol{\sigma}$$, $$\partial f/\partial\alpha$$, $$\partial f/\partial T$$, $$\boldsymbol{\mathsf{C}}$$, and $$\dot{\boldsymbol{\varepsilon}}$$.

Part 6
The continuum elastic-plastic tangent modulus is defined by the following relation

\dot{\boldsymbol{\sigma}} = \boldsymbol{\mathsf{C}}^{\text{ep}}:\dot{\boldsymbol{\varepsilon}} ~. $$ Derive an expression for the elastic plastic tangent modulus using the results you have derived in the previous parts.

Part 7
The $$J_2$$ theory of plasticity also states that the material satisfies the von Mises yield condition

f(\boldsymbol{\sigma},\alpha,T) := \sqrt{\cfrac{3}{2}}\lVert\mathbf{s}\rVert_{} - \sigma_y(\alpha,T) $$ where $$\mathbf{s}$$ is the deviatoric part of the stress $$\boldsymbol{\sigma}$$. Derive an expression for $$\partial f/\partial\boldsymbol{\sigma}$$ in terms of the normal to the yield surface

\mathbf{n} = \cfrac{\mathbf{s}}{\lVert\mathbf{s}\rVert_{}} ~. $$

Part 8
The yield stress $$\sigma_y$$ is given by the Johnson-Cook model

\sigma_y(\alpha,T) = \left[\sigma_0 + B \alpha^n\right] \left[1 - \left(\cfrac{T - T_0}{T_m -T_0}\right)\right] $$ where $$\sigma_0$$ is the initial yield stress, $$B, n$$ are constants, $$T_0$$ is a reference temperature, and $$T_m$$ is the melt temperature. Derive expressions for $$\partial f/\partial \alpha$$, and $$\partial f/\partial T$$ for the von Mises yield condition with the Johnson-Cook flow stress model.

Part 9
Assume that the elastic response of the material is linear, i.e.,

\boldsymbol{\mathsf{C}} = \lambda\boldsymbol{\mathit{1}}\otimes\boldsymbol{\mathit{1}} + 2\mu\boldsymbol{I} ~. $$ Derive the expression for the elastic-plastic tangent modulus for a von Mises yield condition with Johnson-Cook flow stress for a linear elastic material using the expressions that you have derived in the previous parts.

Part 10
Discretize the equations for $$\dot{\mathbf{e}}^p$$ (equation 1), $$\dot{\alpha}$$ (from part 1), and $$\dot{T}$$ (from part 2) using Forward Euler. Use the following notation in your discretized equations:
 * $$\begin{align}

\Delta\gamma & := \dot{\gamma}_n~\Delta t \\ \mathbf{n}_{n} & := \cfrac{\mathbf{s}_n}{\lVert\mathbf{s}_n\rVert_{}} \end{align}$$ where $$\Delta t$$ is the time step, $$\dot{\gamma}_n$$ is the value of $$\dot{\gamma}$$ at $$t = t_n$$, $$\mathbf{s}_n$$ is the value of $$\mathbf{s}$$ at $$t = t_n$$.

Part 11
The Kuhn-Tucker loading-unloading conditions are

\dot{\gamma} \ge 0 ~, f(\boldsymbol{\sigma}, \alpha, T) \le 0 ~; \dot{\gamma} f(\boldsymbol{\sigma},\alpha,T) = 0 ~. $$ Write down a discrete form of the Kuhn-Tucker conditions.

Part 12
In the radial return algorithm, we define a trial elastic state as
 * $$\text{(2)} \qquad

\boldsymbol{\sigma}_{n+1}^{\text{trial}} = \boldsymbol{\sigma}_n + \boldsymbol{\mathsf{C}}:(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_n) $$ where $$\boldsymbol{\sigma}_n,\boldsymbol{\varepsilon}_n$$ are the stress and strain at $$t = t_n$$ and $$\boldsymbol{\sigma}_{n+1},\boldsymbol{\varepsilon}_{n+1}$$ are the values at $$t = t_{n+1}$$. Show that, if the elastic response of the material is linear, equation (2) can be written as
 * $$\text{(3)} \qquad

\boldsymbol{\sigma}_{n+1}^{\text{trial}} = \left[\lambda~\text{tr}(\boldsymbol{\varepsilon}_{n+1})~\boldsymbol{\mathit{1}} + 2\mu~\boldsymbol{\varepsilon}_{n+1}\right] - \left[\lambda~\text{tr}(\boldsymbol{\varepsilon}_n^p)~\boldsymbol{\mathit{1}} + 2\mu~\boldsymbol{\varepsilon}_n^p\right] ~. $$ Hint: Start by showing that

\boldsymbol{\sigma}_{n+1}^{\text{trial}} = \boldsymbol{\mathsf{C}}:(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_n^p) ~. $$

Part 13
Starting from equation (3) show that

\mathbf{s}_{n+1}^{\text{trial}} = \mathbf{s}_n + 2\mu~(\mathbf{e}_{n+1}-\mathbf{e}_n) $$ where $$\mathbf{s}$$ is the deviatoric part of $$\boldsymbol{\sigma}$$ and $$\mathbf{e}$$ is the deviatoric part of $$\boldsymbol{\varepsilon}$$.

Part 14
Show that

\mathbf{s}_{n+1} = \mathbf{s}_{n+1}^{\text{trial}} - 2\mu~\Delta\gamma~\mathbf{n}_n ~. $$ Hint: The stress $$\boldsymbol{\sigma}_{n+1}$$ is given by

\boldsymbol{\sigma}_{n+1} = \boldsymbol{\mathsf{C}}:\boldsymbol{\varepsilon}_{n+1}^e ~. $$ Express this equation in terms of $$\boldsymbol{\varepsilon}_{n+1}$$ and $$\boldsymbol{\varepsilon}_{n+1}^p$$. Then use the discretized equation for $$\dot{\boldsymbol{\varepsilon}}^p$$ (part 10) and the relation for $$\boldsymbol{\mathsf{C}}$$ for isotropic elasticity. Finally compute the deviatoric stress terms after showing that

\text{tr}(\mathbf{s}) = 0 ~. $$

Part 15
The discretized form of the Kuhn-Tucker conditions in conjunction with the consistency condition gives us

f(\boldsymbol{\sigma}_{n+1},\alpha_{n+1},T_{n+1}) = 0 ~. $$ Use this condition and the relations you have derived in the previous sections to arrive at a nonlinear equation in $$\Delta\gamma$$ that can be solved using Newton iterations.

Part 16
Let the nonlinear equation be $$g(\Delta\gamma)$$. Recall that the Newton method requires that we iterate using the formula

\Delta\gamma_{r+1} = \Delta\gamma_r - \cfrac{g(\Delta\gamma_r)} {\cfrac{dg(\Delta\gamma_r)}{d\Delta\gamma}} $$ where $$r$$ is the Newton iteration number. Derive an expression for the derivative of $$g$$ that is required in the above formula.

(You can use  Computational Inelasticity by J.C. Simo and T.J.R. Hughes for pointers.)

Problem 2: Billet Upset Forging
Consider the isothermal upset forging of the cylindrical billet shown in Figure 2.

Assume that the dies are rigid. Also assume that sticking friction is in effect between the billet and the die faces when they are in contact.

The billet has an initial radius of 10 mm and its initial height is 30 mm. The shear modulus of the material is 384.6 MPa, the bulk modulus of the material is 833.3 MPa, the initial yield stress is 1 MPa and the linear hardening modulus is 3 MPa.

Model a quarter of the cylinder using symmetry boundary conditions.

Apply a compressive force of 1 kN to the die.

 * 1) Plot the final shape of the billet. Compare your results  with those shown in Simo and Hughes (Fig. 9.8, p. 325).
 * 2) Plot a curve of the die force (kN) versus the die stroke (mm).   Compare your results with those shown in Simo and Hughes.  Do  you observe any volumetric locking?

(Use an implicit software to solve these problems.)

Problem 3: Taylor Impact Tests
Consider the impact of a cylindrical Taylor impact specimen on a rigid target. The undeformed and deformed profiles of the specimen are shown in in Figure 3.

The initial length of the specimen is 30 mm. The initial diameter is 6 mm. The initial velocity is 188 m/s. The initial temperature is 718 K.

The material of the specimen is OFHC copper. The properties of the bar are (in SI units):

The plastic deformation of the specimen is described by the Johnson-Cook model and $$J_2$$ plasticity. The Johnson-Cook model parameters are (in SI units):

 * Use LS-DYNA to simulate the Taylor impact test. Assume that there is no friction between the anvil and the specimen. Plot the final deformed shape of the specimen and compare that with the experimentally determined shape given in the table below. What differences do you observe and why?
 * {| border="2" cellpadding="4" cellspacing="0"

!Point ! x (mm) ! y (mm) |}
 * 1
 * 0.000000
 * 0.000000
 * 2
 * 5.436409
 * 0.000000
 * 3
 * 4.711554
 * 0.852540
 * 4
 * 4.611804
 * 2.040725
 * 5
 * 4.581879
 * 3.228910
 * 6
 * 4.615129
 * 4.141866
 * 7
 * 4.585204
 * 4.980980
 * 8
 * 4.448878
 * 6.175878
 * 9
 * 4.312552
 * 7.223092
 * 10
 * 4.073150
 * 8.344149
 * 11
 * 3.870324
 * 9.465205
 * 12
 * 3.597672
 * 10.868203
 * 13
 * 3.388196
 * 11.707317
 * 14
 * 3.218620
 * 12.902215
 * 15
 * 3.152120
 * 13.949429
 * 16
 * 2.982544
 * 15.070486
 * 17
 * 2.952618
 * 16.674871
 * 18
 * 0.000000
 * 16.674871
 * 19
 * 0.000000
 * 0.000000
 * 15
 * 3.152120
 * 13.949429
 * 16
 * 2.982544
 * 15.070486
 * 17
 * 2.952618
 * 16.674871
 * 18
 * 0.000000
 * 16.674871
 * 19
 * 0.000000
 * 0.000000
 * 19
 * 0.000000
 * 0.000000
 * 0.000000
 * You can generate a mesh in ANSYS and tranfer it to LS-DYNA if you want.

 * Show whether energy is conserved during your simulation.

title: Nonlinear finite elements/Homework 11


text: File:Billet.png

title: File:Billet.png


text: File:TaylorTest.png

title: File:TaylorTest.png


text: File:NFE HW10Prob1.png

title: File:NFE HW10Prob1.png


text: File:NFE HW9Prob1.png

title: File:NFE HW9Prob1.png


text: File:KirchhoffPlate.png

title: File:KirchhoffPlate.png


text: File:NFE HW8Prob2.png

title: File:NFE HW8Prob2.png


text: File:NFE HW8Prob3.png

title: File:NFE HW8Prob3.png


text: File:NFE HW8Prob3 1.png

title: File:NFE HW8Prob3 1.png


text: File:CompositeShell.png

title: File:CompositeShell.png


text: File:CantileverEndLoad.png

title: File:CantileverEndLoad.png


text: File:ChipHeat.png

title: File:ChipHeat.png


text: File:NFE Homework1 1.png

title: File:NFE Homework1 1.png


text: File:NFE Homework1 2a.png

title: File:NFE Homework1 2a.png


text:
Nonlinear finite elements/Quiz 1

Quiz
''Answer the following questions. You have 15 minutes.''

Heat conduction in an isotropic material with a constant thermal conductivity and no internal heat sources is described by Laplace's equation

\nabla^2 T = 0 ~\qquad~ \text{or,} \qquad \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} = 0 ~. $$

 * 1) Derive a symmetric weak form for the Laplace equation in 1-D (an insulated rod).
 * 2) What are the expressions for the components of the finite element stiffness matrix ($$K_{ij}$$) and the load vector ($$f_i$$) for this 1-D problem?
 * 3) Assume that the one of ends of the rod is maintained at a temperature of $$T_1$$ (which is nonzero) and the other end has a prescribed heat flux of $$Q_2$$. If we discretize the rod into two elements, what does the reduced finite element system of equations look like?  You do not have to work out the terms of the stiffness matrix - just use generic labels.
 * Now, assume that the thermal conductivity of the material varies with temperature. What form does the governing equation take? (We will call this the modified problem.)
 * 1) List the steps needed to solve the modified problem using finite elements.

title: Nonlinear finite elements/Quiz 1


text:
Introduction to vibrations

Introduction
Vibrations (or oscillations) occur primarily in solid objects (or assemblies of objects) as a result of forces acting upon the objects from outside, or from forces originating inside the objects.

Example
An example of an outside force is wind acting upon a building or bridge. An example of an internal force is the cylinders in a gasoline car engine.

Damping vibrations
Soft materials tend to damp out vibrations quickly while harder materials maintain vibrations longer. An object submerged in a dense fluid (liquid or gas) will also have any vibrations quickly dampened by the fluid.

Examples
A rubber band (an "elastic" in the UK) will quickly stop vibrating, after being plucked, since the energy is rapidly dissipated as heat and sound. Glass, however, being quite hard, may be shattered by sounds at the resonant frequency, as the glass is incapable of rapidly damping the sound energy. Thus, the acoustic energy increases until fracture occurs.

Home experiment
Strike a bell, then submerge the bell in water. The bell's vibrations and sound should end almost immediately, due to the damping effect of the water. You may observe waves generated in the water as the vibrations pass from the bell to the water.

Frequency and wavelength
Every vibration has both a frequency (number of vibrations per time period) and a wavelength (time, or sometimes distance, between vibrations). These two values are inversely proportional to one another.

Natural frequency and resonance
Every object "wants" to vibrate at a certain frequency, called the natural or resonance frequency. This frequency is dependent on the object's material, mass, and configuration. Some objects will also resonate at integer multiples of the natural frequency (twice, three times, etc.).

Example
Each bell has a specific frequency at which it will vibrate, and thus gives off a certain note when struck.

Home experiment
Try striking an empty wine glass with your fingernail. You should hear an audio vibration at a certain pitch. This pitch is the resonant frequency of the glass. Now partially fill the glass with water and repeat the experiment. The resonant frequency will have changed for the new configuration.

Cancelling vibrations
Vibrations may be thought of as a wave phenomenon, meaning that each vibration wave has an associated crest and trough. If two vibrations of the same wavelength/frequency are superimposed (combined), they will either add to each other, if they are in phase (the crests and troughs match), or they will tend to cancel each other, if out of phase (the crests from one wave match the troughs from another). In fact, it is theoretically possible to completely cancel all vibrations by using two waves of identical magnitude and frequency/wavelength.

Examples
Noise canceling headphones are available which electronically shift sounds half a wavelength behind the input sounds and thus cancel the noise. Some mufflers employ a mechanical system to accomplish the same result. Car engines almost always have an even number of cylinders, so that half the cylinders are up when the other half is down, in order to minimize vibration.

Speed of propagation of vibrations
Many vibrations propagate (spread) at, or near, the speed of sound in the given material (the speed of sound, unlike the speed of light, is not constant). Some vibrations propagate several times faster however, via shock waves, and other vibrations propagate much more slowly.

Example
The Earth undergoes a very slow oscillation, known as precession, which causes it to "wobble" about the poles every 25,800 years. The most noticeable effect is that the North Star changes over tens of thousands of years. This extremely slow vibration is caused by tidal forces exerted by the Moon and Sun.

Uses for vibrations

 * Diagnosing mechanical systems: Each type of wear, failure, or loss of lubrication in a mechanical system changes it's configuration and/or acoustic damping properties and thus it's resonant frequency. Experienced operators can recognize differences in sounds made by a machine and know that those differences indicate that specific maintenance operations must be performed on the machine.

 * Acoustic lubrication: Since vibrating objects are only in contact with each other part of the time, they exhibit reduced friction between them. For optimal effect, the frequency of the vibration must match the resonant frequency of the objects.

Dangers of vibrations

 * Vibrations may lead to the failure of a component or entire mechanical system. During an earthquake, for example, vibrations may lead to structural failure of buildings, causing loss of life.  Wind may also induce vibrations in improperly designed buildings, such as the Tacoma Narrows Bridge, which collapsed due to resonant frequency vibrations in only moderate winds.

title: Introduction to vibrations


text:
Wikiversity:Project of the...

This page is intended to allow for the deciding on a project of the week/month/quarter, etc.

Weeks, months, or?
Here please discuss how often the "Project of the..." should be changed. What say you?
 * How about a month? A week seems too short, and a quarter seems to long. --Remi 04:58, 16 August 2007 (UTC)

Project to focus on
Please discuss here what projects you may be interested in utilizing for Wikiversity wide collaboration. Projects might include organizing, a certain course, a certain specific project, etc.

title: Wikiversity:Project of the...


text:
Category:Pages unlisted anywhere

There are a lot of these, whether you have found them or not.

See: Special:Lonelypages Category:Wikiversity

Pages can be added to this category using.

title: Category:Pages unlisted anywhere


text: Template:Orphan

title: Template:Orphan


text:
Wikiversity:Candidates for Custodianship/RM77

Hello everybody. Before all, I'm sorry for the english mistakes I could make, French is my native language. I wish I could have the Custodian rights on EN.WV. I recap my reasons, it's not for personal accomplishment :) : Indeed, I'm a bureaucrat on FR.WV, and a regular contributor on it since december (opening of French version) with about 4K edits. I'm also sysop on FR.WB, and a regular user of FR.WP.
 * There is an increasing amount of vandalism which is let daily on the main namespace (i.e : Drilling Engineering)
 * One more admin is worth it on a relatively new project like Wikiversity
 * This could help the collaboration between FR.WV and EN.WV.

Questions and discussion

 * Hi RM77. You only have 7 edits on this project... perhaps you might want to wait a while (and get familiar with simple stuff like the delete template when patrolling for vandalism). --SB_Johnny | PA! 14:17, 16 August 2007 (UTC)
 * I note your comment about increasing vandalism. I wasn't aware that we have had increasing vandalism recently - indeed, the opposite. If I've been missing a lot of vandalism, I'd like to know about it - perhaps you can point me to the regularly vandalised pages? McCormack 15:58, 21 August 2007 (UTC)
 * See here and here. Keep in mind that I've done a couple hardblocks too, using CU data from wikibooks and the CU logs for Wikipedia... it's still fairly low-level, but definitely a growing problem. You probably don't notice it much because we're blessed with a fairly large and active cadre of Custodians, so many hands are making light work for all of us. --SB_Johnny | PA! 17:39, 21 August 2007 (UTC)
 * Hmm... I note that 39% of this month's deletions are mine. I think I notice... We interpret the volume and severity differently. -- McCormack 17:52, 21 August 2007 (UTC)
 * It's less a question of volume than percentage of volume. Special:Statistics shows only about 4500 pages likely to be content. If you look at that and then look back on those logs (as well as reversions), it's actually not so tiny an issue. Also keep in mind that all wm projects tend to have less vandalism during late summer in the Northern Hemisphere, because the kids are mostly outdoors, rather than bored to death in the house. --SB_Johnny | PA! 18:03, 21 August 2007 (UTC)
 * Yes - I agree with all of that. -- McCormack 18:32, 21 August 2007 (UTC)
 * Hi RM77. Though custodianship is "no big deal", I think you really need to learn how a project works before being helpful as a custodian. I know you're a very good user and I trust you based on your work on the French-language Wikiversity, but I agree you should wait before requesting custodian tools here. On a sidenote, I'd love to see more collaboration between fr.wikiversity and en.wikiversity, but you don't really need custodian tools to do that, do you? Cheers, guillom 13:16, 29 August 2007 (UTC)
 * I agree, you should know the project better. E.g. I am custodian at the German wikiversity. So, one would think then it is about the same. But there are differences between these two Wikiversities (as recommended below, join the IRC channel and you can acclimate faster). Erkan Yilmaz ( evaluate me!, discussion ) 17:45, 29 August 2007 (UTC)

 * The major method by which en-Wikiversity tries to lower the barrier to custodianship is by allowing existing custodians to mentor probationary custodians. If you can get one existing custodian to trust you and be your mentor, then you get a month with the admin tools and can then demonstrate to the rest of the community that you can be trusted. One way to get to know some custodians is to participate at #wikiversity-en. --JWSchmidt 14:08, 29 August 2007 (UTC)
 * I would recommend as JWSchmidt said. He is my mentor I got to know him at the IRC channel and can ask for things and get them cleared fast. Visit us at the IRC channel. Erkan Yilmaz ( evaluate me!,  discussion ) 17:45, 29 August 2007 (UTC)

❌ User withdrawn. --Goldenburg111 22:06, 30 March 2014 (UTC)
 * His last edit was on 16th August. Should we write him an email to awake his interest again ? LOL Erkan Yilmaz ( evaluate me!, discussion ) 17:46, 29 August 2007 (UTC)
 * No mail will be needed :) I had been busy last days.
 * I agree with you all, I think I asked Custodianship too early. Guillom, I'll make all my possible to divide my time between FR.WP and FR.WV, so that's why the collaboration between FR.WV and EN.WV is not that much advanced. I swear I make my possible :)
 * If being a mentor is one of the way to get the Custodians tools, I will try it. I'm quite used with admin tools...
 * I go around here when I can. Forgive me for my late answer :) RM77 <=> We talk. 20:07, 29 August 2007 (UTC)


text:
User talk:RM77

You may better choose this one :)

speedy deletion
What are you trying to do with edits like this edit?

Please do not mark pages for speedy deletion unless you read the page, try to understand it and provide a reason for the deletion. We only speedy delete content that is causing harm to Wikiversity. If a page is not causing harm, there is no need to be in a hurry. Saying "orphan talk page" does not indicate that you have tried to figure out why it is an orphan and if the contents of the page should be kept. Do not delete the contents of pages that you mark for deletion.--JWSchmidt 17:13, 16 August 2007 (UTC)

title: User talk:RM77


text: User:Fica Blok38

title: User:Fica Blok38


text: Category:User sr

title: Category:User sr


text: Template:User sr

title: Template:User sr


text:
User talk:Fica Blok38

Hello Fica Blok38, and welcome to Wikiversity! If you need help, feel free to visit my talk page, or contact us and ask questions. When leaving comments for others to read, remember to sign and date; it helps everyone keep track of who is writing messages. The signature icon in the edit window makes it simple. To help you get started:
 * Take a guided tour and learn how to edit;
 * Explore our learning projects;
 * Browse our portals, schools, and research pages; and
 * Read and help develop our community policies.

Also, don't forget to experiment with the links to your left. Be bold, and see you around Wikiversity! --HappyCamper 14:31, 16 August 2007 (UTC)

title: User talk:Fica Blok38


text:
User talk:Pusey

Hello Pusey, and welcome to Wikiversity! If you need help, feel free to visit my talk page, or contact us and ask questions. When leaving comments for others to read, remember to sign and date; it helps everyone keep track of who is writing messages. The signature icon in the edit window makes it simple. To help you get started:
 * Take a guided tour and learn how to edit;
 * Explore our learning projects;
 * Browse our portals, schools, and research pages; and
 * Read and help develop our community policies.

Also, don't forget to experiment with the links to your left. Be bold, and see you around Wikiversity! --HappyCamper 14:32, 16 August 2007 (UTC)

title: User talk:Pusey


text:
User talk:Vonjy TODISOA

Hello Vonjy TODISOA, and welcome to Wikiversity! If you need help, feel free to visit my talk page, or contact us and ask questions. When leaving comments for others to read, remember to sign and date; it helps everyone keep track of who is writing messages. The signature icon in the edit window makes it simple. To help you get started:
 * Take a guided tour and learn how to edit;
 * Explore our learning projects;
 * Browse our portals, schools, and research pages; and
 * Read and help develop our community policies.

Also, don't forget to experiment with the links to your left. Be bold, and see you around Wikiversity! --HappyCamper 14:31, 16 August 2007 (UTC)

title: User talk:Vonjy TODISOA


text:
User talk:Denny12~enwikiversity

Hello Denny12, and welcome to Wikiversity! If you need help, feel free to visit my talk page, or contact us and ask questions. When leaving comments for others to read, remember to sign and date; it helps everyone keep track of who is writing messages. The signature icon in the edit window makes it simple. To help you get started:
 * Take a guided tour and learn how to edit;
 * Explore our learning projects;
 * Browse our portals, schools, and research pages; and
 * Read and help develop our community policies.

Also, don't forget to experiment with the links to your left. Be bold, and see you around Wikiversity! --HappyCamper 14:31, 16 August 2007 (UTC)

Your account will be renamed
Hello,

The developer team at Wikimedia is making some changes to how accounts work, as part of our on-going efforts to provide new and better tools for our users like cross-wiki notifications. These changes will mean you have the same account name everywhere. This will let us give you new features that will help you edit and discuss better, and allow more flexible user permissions for tools. One of the side-effects of this is that user accounts will now have to be unique across all 900 Wikimedia wikis. See the announcement for more information.

Unfortunately, your account clashes with another account also called Denny12. To make sure that both of you can use all Wikimedia projects in future, we have reserved the name Denny12~enwikiversity that only you will have. If you like it, you don't have to do anything. If you do not like it, you can pick out a different name.

Your account will still work as before, and you will be credited for all your edits made so far, but you will have to use the new account name when you log in.

Sorry for the inconvenience.

Yours, Keegan Peterzell Community Liaison, Wikimedia Foundation 23:31, 17 March 2015 (UTC)

Renamed
 This account has been renamed as part of single-user login finalisation. If you own this account you can |log in using your previous username and password for more information. If you do not like this account's new name, you can choose your own using this form after logging in: . -- Keegan (WMF) (talk) 06:06, 19 April 2015 (UTC)

title: User talk:Denny12~enwikiversity


text:
User:Val.samonis

INTRO: Educated internationally, Val Samonis (PhD, CPC) gained over two decades of cross-functional research, teaching/training, knowledge management, and advisory experience combining theoretical and hands-on empirical approaches to business, technology, and education frontiers globally. Expert in the Internet-based and other high-tech methods of knowledge management, teaching, and problem solving, he is an outside-of-the-box thinker and people-oriented intercultural communicator (fluent in English, German, Russian, Polish, Lithuanian; good knowledge of other languages) with many important contacts in emerging markets, Europe, and North America.

CORE EXPERIENCE: Dr. Val Samonis has worked with top business and education leaders and reformers in many countries, e.g. Harvard President Larry Summers, Nobel Laureates in Economics and Finance (Tobin, Arrow, Solow, Leontieff, Klein), Advisory to The British Commonwealth of Nations’ “Commonwealth of Learning” on modern online and blended university education, The Royal Roads University President’s Globalization Expert Commission, Advisory to the Academic Dean, Lansbridge University, Polish Deputy Prime Minister L. Balcerowicz, Czech Deputy Prime Minister P. Mertlik, other top experts globally; lectured internationally on trade, investment, knowledge economy, corporate governance and enterprise restructuring, competitiveness, and other business reform and technology issues in the global economy; and has been extensively published (list of his publications is over 30 pages long). He managed and/or worked in international research and advisory ("blue ribbon") teams sponsored by the Base of the Pyramid Initiative (Cornell), The African Capacity Building Foundation (G-7 countries and global institutions), The Hudson Institute, World Bank, The OECD-World Bank Private Sector Advisory Group on Corporate Governance, CASE Warsaw, Soros Foundations, The Center for European Integration Studies (ZEI Bonn), Joint Committee on Corporate Governance (Canada), ACCC/CIDA (Canada), and a number of governments, e.g. in Baltic States, Poland, Canada. Val served as an advisor to the Czech Government, the Lithuanian Parliament, international organizations (e.g. UNCTAD, WTO), and multinational corporations (e.g. CARE, Medley Global Advisors, Andersen Consulting). Dr. Samonis has been teaching online, onsite, and his DMI(SM) dual-mode comprehensive and specialized courses in Global Management/Business, Frontiers in Emerging Markets (technological and geographical), Corporate Governance and Finance, Global Knowledge Economy, as well as supervising graduate (MBA, MA, PhD) students at the University of Toronto (Canada), University of Maryland (USA), Hawaii Pacific University USA), Royal Roads University (Canada), Center for University Studies (USA & Mexico), Warsaw School of Economics (Poland), University of Tasmania (Australia), The Center for European Integration Studies (Germany), Vilnius University (Lithuania), and other institutions internationally.

title: User:Val.samonis


text:
User:A.K.Karthikeyan

Hello, I am a beta in everthing that's science, technology and logic. Hoping to get an alpha in something.

title: User:A.K.Karthikeyan


text: User:Multipleidentitynumberthree

title: User:Multipleidentitynumberthree


text: User talk:Multipleidentitynumberthree

title: User talk:Multipleidentitynumberthree