Difference between revisions of "Solid Mechanics"
(→Connection between plane stress and plane strain) |
(→Two-dimensional stress distributions) |
||
| Line 49: | Line 49: | ||
1 & \nu & 0 \\ | 1 & \nu & 0 \\ | ||
\nu & 1 & 0 \\ | \nu & 1 & 0 \\ | ||
| − | 0 & 0 & \frac{1 | + | 0 & 0 & \frac{1}{2}(1-\nu) |
\end{bmatrix} | \end{bmatrix} | ||
\begin{bmatrix} | \begin{bmatrix} | ||
| Line 80: | Line 80: | ||
1-\nu & \nu & 0 \\ | 1-\nu & \nu & 0 \\ | ||
\nu & 1-\nu & 0 \\ | \nu & 1-\nu & 0 \\ | ||
| − | 0 & 0 & \frac{1-2\nu | + | 0 & 0 & \frac{1}{2}(1-2\nu) |
\end{bmatrix} | \end{bmatrix} | ||
\begin{bmatrix} | \begin{bmatrix} | ||
| Line 98: | Line 98: | ||
====Connection between plane stress and plane strain==== | ====Connection between plane stress and plane strain==== | ||
| − | For '''isotropic''' materials with elastic modulus $E$ and Poisson's ratio $\nu$ it is possible to go from plane stress to plane strain, or vice-versa, by replacing $E$ and $\nu$ in the stress-strain matrix with a fictitious modulus $E^*$ and fictitious Poisson ratio $\nu^*$. This allows us to "reuse" a plane stress program to solve plane strain or again vice-versa (as long as the material is isotropic). | + | For '''isotropic''' materials with elastic modulus $E$ and Poisson's ratio $\nu$ it is possible to go from plane stress to plane strain, or vice-versa, by replacing $E$ and $\nu$ in the stress-strain matrix with a fictitious modulus $E^*$ and fictitious Poisson ratio $\nu^*$. This allows us to "reuse" a plane stress program to solve plane strain or again vice-versa (as long as the material is isotropic). A few exercises on this topic are given [http://www.colorado.edu/engineering/cas/courses.d/IFEM.d/IFEM.Ch14.d/IFEM.Ch14.pdf at this link (page 13)]. |
To go from plane stres'''s''' ($s$) to plane strai'''n''' ($n$) insert the fictitious quantities | To go from plane stres'''s''' ($s$) to plane strai'''n''' ($n$) insert the fictitious quantities | ||
Revision as of 23:43, 30 October 2016
Contents
Basic equations of elasticity
To determine the distribution of static stresses and displacements in a solid body we must obtain a solution (either analytically or numerically) to the basic equations of the theory of elasticity, satisfying the boundary conditions on forces and/or displacements. For a general three dimensional solid object these equations are:
- strain-displacement equations ($6$)
- stress-strain equations ($6$)
- equations of equilibrium ($3$)
where the number in brackets indicates the number of equations. We see there are $15$ equations for $15$ unknown variables ($6$ strains, $6$ stresses and $3$ displacements). In two dimensions this simplifies to $8$ equations with $2$ displacements, $3$ stresses, and $3$ strains.
Stress formulation
Strain-displacement equations
Under the action of applied forces, a point in the solid originally at $(x,y,z)$ moves to position $(X,Y,Z)$. This movement can be described completely by the displacement vector \begin{equation} \boldsymbol{u}(\boldsymbol{x}) = \{u_x(x,y,z) \quad u_y(x,y,z) \quad u_z(x,y,z)\}, \end{equation} the brackets $\{\}$ indicating that this is a column vector.
Navier equation describes the dynamic of a solid through the displacement vector field $\b{u}$. The equation is as follows \begin{equation} \rho \frac{\partial^2 \b{u}}{\partial t^2} = \mu \nabla^2 \b{u} + (\lambda + \mu) \nabla(\nabla \cdot \b{u}) + \b{F}\ , \end{equation} where $\mu$ and $\lambda$ are Lamé constants, $\rho$ is the object density and $\b{F}$ are the external forces.
Stress-strain equations
Two-dimensional stress distributions
Many problems in elasticity can be simplified as two-dimensional problems described by plane theory of elasticity. In general there are two types of problems we may encounter in plane analysis: plane stress and plane strain. The first problem arises in analysis of thin plates loaded in the plane of the plate, while the second is used for elongated bodies of constant cross section subject to uniform loading.
Plane stress
Plane stress distributions build on the assumption that the normal stress and shear stresses directed perpendicular to the $x$-$y$ plane are assumed zero: \begin{equation}\label{eq:pstress_assump} \sigma_{zz} = \sigma_{zx} = \sigma_{zy} = 0. \end{equation} It is also assumed that the stress components do not vary through the thickness of the plate (the assumptions do violate some compatibility conditions, but are still sufficiently accurate for practical applications if the plate is thin).
Using (\ref{eq:pstress_assump}) the three-dimensional Hooke's law can be reduced to: \begin{equation}\label{eq:planestress} \begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{bmatrix} = \frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1}{2}(1-\nu) \end{bmatrix} \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \varepsilon_{xy} \end{bmatrix} \end{equation}
Plane strain
The plane strain problem arises in analysis of walls, dams, tunnels where one dimension of the structure is very large in comparison to the other two dimensions ($x$- and $y$- coordinates). It is also appropriate for small-scale problems such as bars and rollers compressed by forces normal to their cross section. In all such problems the body may be imagined as a prismatic cylinder with one dimension much larger that the other two. The applied forces act in the $x$-$y$ plane and do not vary in the $z$ direction, leading to the assumption \begin{equation} \frac{\partial}{\partial z} = u_z = 0. \end{equation}
With the above assumption it follows immediately that \begin{equation} \varepsilon_{zz} = \varepsilon_{zx} = \varepsilon_{zy} = 0. \end{equation} The three dimensional Hooke's law can now be reduced to \begin{equation}\label{eq:planestrain} \begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{bmatrix} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \frac{1}{2}(1-2\nu) \end{bmatrix} \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \varepsilon_{xy} \end{bmatrix} \end{equation} \begin{equation}\label{eq:sigmazz} \sigma_{zz} = \nu(\sigma_{xx}+\sigma_{yy}) \end{equation} \begin{equation} \sigma_{yz} = \sigma_{zx} = 0 \end{equation} The reason $\sigma_{zz}$ is not included in the matrix stress-strain equation (\ref{eq:planestrain}) is because it is linearly dependent on the normal stresses $\sigma_{xx}$ and $\sigma_{zz}$ as shown in (\ref{eq:sigmazz}).
Connection between plane stress and plane strain
For isotropic materials with elastic modulus $E$ and Poisson's ratio $\nu$ it is possible to go from plane stress to plane strain, or vice-versa, by replacing $E$ and $\nu$ in the stress-strain matrix with a fictitious modulus $E^*$ and fictitious Poisson ratio $\nu^*$. This allows us to "reuse" a plane stress program to solve plane strain or again vice-versa (as long as the material is isotropic). A few exercises on this topic are given at this link (page 13).
To go from plane stress ($s$) to plane strain ($n$) insert the fictitious quantities \begin{equation}\label{eq:ston1} E_n^* = \frac{E_s}{1-\nu_s^2}, \end{equation} \begin{equation}\label{eq:ston2} \nu_n^* = \frac{\nu_s}{1-\nu_s}. \end{equation}
Substitution from plane stress to plane strain
We start with the plane stress matrix with the inserted fictitious values $E^*$ and $\nu^*$ \[ \frac{E^*}{1-{\nu^*}^2} \begin{bmatrix} 1 & \nu^* & 0 \\ \nu^* & 1 & 0 \\ 0 & 0 & \frac{1}{2}(1-\nu^*) \end{bmatrix} \] and make the above substitutions (\ref{eq:ston1}) and (\ref{eq:ston2}) leading to: \[ \frac{E}{\left(1-\nu^2\right)\left(1-\left(\frac{\nu}{1-\nu}\right)^2\right)} \begin{bmatrix} 1 & \frac{\nu}{1-\nu} & 0 \\ \frac{\nu}{1-\nu} & 1 & 0 \\ 0 & 0 & \frac{1}{2}(1-\frac{\nu}{1-\nu}) \end{bmatrix}. \] We can then use the rule to convert sums of squares into products as well as bring the factor $1/(1-v)$ out of the matrix: \[ \frac{E}{\left(1-\nu\right)\left(1+\nu\right)\left(1-\frac{\nu}{1-\nu}\right)\left(1+\frac{\nu}{1-\nu}\right)\left(1-\nu\right)} \begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \frac{1}{2}(1-2\nu) \end{bmatrix}. \] By joining some of the factors to a common denominator and rearranging leads to \[ \frac{E\left(1-\nu\right)\left(1-\nu\right)}{\left(1-\nu\right)\left(1+\nu\right)\left(1-2\nu\right)1\left(1-\nu\right)} \begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \frac{1}{2}(1-2\nu) \end{bmatrix}. \] The final step is canceling the factors that occur in both the numerator and denominator \[ \frac{E}{\left(1+\nu\right)\left(1-2\nu\right)} \begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \frac{1}{2}(1-2\nu) \end{bmatrix} \] which leads exactly to the relationship for plane strain given in (\ref{eq:planestrain}).
In the opposite case we want to go from plane strain ($n$) to plane stress ($s$4) we can use:
\begin{equation}\label{eq:ntos1}
E_s^* = \frac{E_n(1+2\nu_n)}{(1+\nu_n)^2},
\end{equation}
\begin{equation}\label{eq:ntos2}
\nu_s^* = \frac{\nu_n}{1+\nu_n}.
\end{equation}
Substitution from plane strain to plane stress
Note that we have again omitted the indexes ($s$) and ($n$) since our wish is expressed by the bold title.
We start with the plane strain matrix with the inserted fictitious values $E^*$ and $\nu^*$ \[ \frac{E^*}{\left(1+\nu^*\right)\left(1-2\nu^*\right)} \begin{bmatrix} 1-\nu^* & \nu^* & 0 \\ \nu^* & 1-\nu^* & 0 \\ 0 & 0 & \frac{1}{2}(1-2\nu^*) \end{bmatrix} \] and make the above substitutions (\ref{eq:ntos1}) and (\ref{eq:ntos2}) leading to: \[ \frac{E(1+2\nu)}{\left(1+\nu\right)^2\left(1+\frac{\nu}{1+\nu}\right)\left(1-2\frac{\nu}{1+\nu}\right)} \begin{bmatrix} 1-\frac{\nu}{1+\nu} & \frac{\nu}{1+\nu} & 0 \\ \frac{\nu}{1+\nu} & 1-\frac{\nu}{1+\nu} & 0 \\ 0 & 0 & \frac{1}{2}(1-2\frac{\nu}{1+\nu}) \end{bmatrix}. \] Writing some of the sums with common denominators and rearranging leads to \[ \frac{E(1+2\nu)(1+\nu)(1+\nu)}{\left(1+\nu\right)^2\left(1+2\nu\right)\left(1-\nu\right)} \begin{bmatrix} \frac{1}{1+\nu} & \frac{\nu}{1+\nu} & 0 \\ \frac{\nu}{1+\nu} & \frac{1}{1+\nu} & 0 \\ 0 & 0 & \frac{1}{2}\left(\frac{1-\nu}{1+\nu}\right) \end{bmatrix}. \] We can also bring the factor $1/(1+\nu)$ out from the matrix components. After canceling all factors that occur in both denominator and numerator we are left with \[ \frac{E}{\left(1+\nu\right)\left(1-\nu\right)} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1}{2}(1-\nu) \end{bmatrix}. \] Rewriting the product of sums as a sum of squares gives us the matrix for plane stress in (\ref{eq:planestress}).
Equations of equilibrium
Tonti diagram
Point contact on a 2D half-plane
A starting point to solve problems in contact mechanics is to understand the effect of a point-load applied to a homogeneous, linear elastic, isotropic half-plane. This problem may be defined either as plane stress or plain strain (for solution with FreeFem++ we have choosen the latter). The traction boundary conditions for this problem are: \begin{equation}\label{eq:bc} \sigma_{xy}(x,0) = 0, \quad \sigma_{yy}(x,y) = -P\delta(x,y) \end{equation} where $\delta(x,y)$ is the Dirac delta function. Together these boundary conditions state that there is a singular normal force $P$ applied at $(x,y) = (0,0)$ and there are no shear stresses on the surface of the elastic half-plane.
The analytical relations for the stresses can be found from the Flamant solution (stress distributions in a linear elastic wedge loaded by point forces a the tip. When the "wedge" is flat we get a half-plane. The derivation uses polar coordinates.) and are given as: \begin{equation} \sigma_{xx} = -\frac{2P}{\pi} \frac{x^2y}{\left(x^2+y^2\right)^2}, \end{equation} \begin{equation} \sigma_{xy} = -\frac{2P}{\pi} \frac{y^3}{\left(x^2+y^2\right)^2}, \end{equation} \begin{equation} \sigma_{yy} = -\frac{2P}{\pi} \frac{xy^2}{\left(x^2+y^2\right)^2}, \end{equation} for some point $(x,y)$ in the half-plane. From this stress field the strain components and thus the displacements $(u,v)$ can be determined. The displacements are given by \begin{align} u_x &= -\frac{P}{4\pi\mu}\left((\kappa-1)\theta - \frac{2xy}{r^2}\right), \label{eq:dispx}\\ u_y &= -\frac{P}{4\pi\mu}\left((\kappa+1)\log r - \frac{2y^2}{r^2}\right), \label{eq:dispy} \end{align} where $$r = \sqrt{x^2+y^2}$$ and $$\tan \theta = \frac{x}{y}.$$ The symbol $\kappa$ is known as Dundars constant and is defined as \[ \kappa = \begin{cases} 3 - 4\nu & \qquad \text{plane strain}, \\ \cfrac{3 - \nu}{1+\nu} & \qquad \text{plane stress}. \end{cases} \] The last remaining symbol is $\mu$ which represents the shear modulus (sometimes also denoted with $G$).
Numerical solution with FreeFem++
Due to the known analytical solution the point-contact problem can be used for benchmarking numerical PDE solvers in terms of accuracy (as well as computational efficiency). The purpose of this section is to compare the numerical solution obtained by FreeFem++ with the analytical solution, as well as provide a reference numerical solution for the C++ library developed in our laboratory.
For purposes of simplicity we limit ourselves to the domain $(x,y) \in \Omega = [-1,1] \times[-1,-0.1]$ and prescribe Dirichlet displacement on the boundaries $\Gamma_D$ from the known analytical solution (\ref{eq:dispx}, \ref{eq:dispy}). This way we avoid having to deal with the Dirac delta traction boundary condition (\ref{eq:bc}). The problem can be describes as find $\boldsymbol{u(\boldsymbol{x})}$ that satisfies \begin{equation} \boldsymbol{\nabla}\cdot\boldsymbol{\sigma}= 0 \qquad \text{on }\Omega \end{equation} and \begin{equation} \boldsymbol{u} = \boldsymbol{u}_{\text{analytical}} \qquad \text{on }\Gamma_D \end{equation} where $\boldsymbol{u}_\text{analytical}$ is given in equations (\ref{eq:dispx}) and (\ref{eq:dispy}).
To solve the point-contact problem in FreeFem++ we must first provide the weak form of the balance equation: \begin{equation*} \boldsymbol{\nabla}\cdot\boldsymbol{\sigma} + \boldsymbol{b} = 0. \end{equation*} The corresponding weak formulation is \begin{equation}\label{eq:weak} \int_\Omega \boldsymbol{\sigma} : \boldsymbol{\varepsilon}(\boldsymbol{v}) \, d\Omega - \int_\Omega \boldsymbol{b}\cdot\boldsymbol{v}\,d\Omega = 0 \end{equation} where $:$ denotes the tensor scalar product (tensor contraction), i.e. $\boldsymbol{A}:\boldsymbol{B} =\sum_{i,j} A_{ij}B_{ij}$. The vector $\boldsymbol{v}$ is the test function or so-called "virtual displacement".
Equation (\ref{eq:weak}) can be handed to FreeFem++ with the help of Voigt or Mandel notation, that reduces the symmetric tensors $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}$ to vectors. The benefit of Mandel notation is that it allows the tensor scalar product to be performed as a scalar product of two vectors. For this reason we create the following macros:
macro u [ux,uy] // displacements macro v [vx,vy] // test function macro b [bx,by] // body forces macro e(u) [dx(u[0]),dy(u[1]),(dx(u[1])+dy(u[0]))/2] // strain (for post-processing) macro em(u) [dx(u[0]),dy(u[1]),sqrt(2)*(dx(u[1])+dy(u[0]))/2] // strain in Mandel notation macro A [[2*mu+lambda,mu,0],[mu,2*mu+lambda,0],[0,0,2*mu]] // stress-strain matrix
The weak form (\ref{eq:weak}) can then be expressed naturally in FreeFem++ syntax as
int2d(Th)((A*em(u))'*em(v)) - int2d(Th)(b'*v)
Stress and displacement fields
Convergence studies
For the convergence between analytical and numerical solution we vary the number of nodes by increasing the grid size in both $x$- and $y$- directions simultaneously by powers of two from $2^2$ (16 nodes all together) to $2^7$ (16384 nodes all together). The $L^2$ error norm is used to measure the "difference" between solutions. Since the displacements are the variable that obtain from FreeFem++, we use the displacement magnitude $|\boldsymbol{u}| = \sqrt{u_x^2+u_y^2}$ to define our $L^2$-error norm. The exact equation we have used is \begin{equation} L^2\text{-norm} = \sqrt{\frac{\int_\Omega (|\boldsymbol{u_{\text{numerical}}}|-|\boldsymbol{u_{\text{analytical}}}|)^2d\Omega}{\int_\Omega|\boldsymbol{u_{\text{analytical}}}|^2d\Omega}}. \end{equation} Results are shown in Figure 1.
Contact between parallel cylinders
References
- Theory of matrix structural analysis
