User:Fem.tm7/HW3

Problem 3.1: Solving General 1-D Model with Weighted Residual Form and Polynomial Basis Functions

Problem Description

Consider $\left\{b_{j}(x);j=0,1,...,n\right\}$ when $\displaystyle b_{j}(x)=(x+k)^{j}$

Choose $\displaystyle k=1$, such that $\displaystyle b_{j}(x=0)\neq 0$

1. Let n=2 $\rightarrow$ ndof = n+1 = 2+1 = 3 with $\mathbf {d} =\left\{d_{j};j=0,...,n\right\}_{(n+1)\times 1}$.
2. Find 2 equations that enforce the boundary conditions for

$u^{h}(x)=\sum _{j=0}^{n}{d_{j}b_{j}(x)}$.

(3.1.1)

3. Find 1 more equation (i.e. j = 0,1,2) to solve for $\mathbf {d} =\left\{d_{j}\right\}_{3\times 1}$ by projecting the residue, $\displaystyle P(u^{h})$, on a basis function, $\displaystyle b_{k}(x)$, with k = 0,1,2 such that the additional equation is linearly independent from the equations in part 2.

• Notes regarding residue:

• The general form of the residue is defined as $\displaystyle P(u^{h})$. The exact solution is expressed as $\displaystyle P(u)=0$.
  
• $\displaystyle u^{h}$ approximates $\displaystyle u$. In other words $\displaystyle u^{h}\cong u$. Thus, $P(u^{h})\neq 0\quad \forall x\in \Omega$ (for all values of x in the domain).

• Notes regarding projection:

• $\displaystyle \int _{\Omega }b_{i}(x)P\left(u^{h}(x)\right)dx=0\quad \Rightarrow \quad \int _{0}^{1}b_{k}(x)P(u^{h})dx=0$
• $\displaystyle \mathbf {b} _{i}\cdot \mathbf {P} (\mathbf {v} )=0\quad \Rightarrow \quad \mathbf {b} _{i}\cdot \mathbf {P} (\mathbf {v} )\ =\ <\mathbf {b} _{i},\mathbf {P} (\mathbf {v} )>$
  
  
• Combining the previous two equations yields

$\displaystyle <b_{k},P(u^{h})>=\int _{0}^{1}b_{k}(x)P(u^{h})dx$.

(3.1.2)

4. Display 3 equations in matrix form, $\mathbf {K} \cdot \mathbf {d} =\mathbf {F}$.
5. Solve for $\mathbf {d}$.
6. Construct $\displaystyle u_{n}^{h}(x)$ and plot versus the exact solution, $\displaystyle u(x)$.
7. Repeat steps 1 through 6 for:

• n = 4
• n = 6

Given

The data set for the general 1-D model with "simple" boundary conditions (G1DM1.0/D2) is given in Vu-Quoc, lecture 12, page 1.

$\displaystyle \Omega =]{0,1}[$

$\displaystyle a_{2}=2$

$\displaystyle f=3$

Static case since $\displaystyle {\frac {\partial ^{s}u}{\partial t^{s}}}(x,t)=0$.

The natural (eq 3.1.3) and essential (ed 3.1.4) boundary conditions are defined as

$-{\frac {d{u}^{h}}{dx}}(0)=4$

(3.1.3)

$\displaystyle u(1)=0$

(3.1.4)

Solution

For the case when $\displaystyle n=2$

Parts 1 & 2

$\left\{b_{j}(x);j=0,1,...,n\right\}$ where $\displaystyle b_{j}(x)=(x+k)^{j}$

$\Rightarrow \quad b_{0}(x)=(x+k)^{0}=1,b_{1}(x)=(x+k)^{1}=(x+k),\ and\ b_{2}(x)=(x+k)^{2}$.

Choosing $k=1$ to avoid $b_{j}'(x)=0$ yields

$\quad b_{0}(x)=1,b_{1}(x)=(x+1),\ and\ b_{2}(x)=(x+1)^{2}$.

The natural boundary condition can be implemented by differentiating $\displaystyle u^{h}(x)$ with respect to $\displaystyle x$ and taking its value at x=0.

${\frac {d{u}^{h}}{dx}}(x)=\sum _{j=0}^{2}{d_{j}b_{j}'(x)}=d_{0}b_{0}'(x)+d_{1}b_{1}'(x)+d_{2}b_{2}'(x)=d_{1}+2d_{2}(x+1)$

${\frac {d{u}^{h}}{dx}}(0)=\sum _{j=0}^{2}{d_{j}b_{j}'(0)}=d_{0}b_{0}'(0)+d_{1}b_{1}'(0)+d_{2}b_{2}'(0)=d_{1}+2d_{2}=-4$

So the relevant equation is

$\displaystyle d_{1}+2d_{2}=-4$

(3.1.5)

The essential boundary condition, $\displaystyle u^{h}(1)=0$, is implemented as follows:

$u^{h}(x)=\sum _{j=0}^{2}{d_{j}b_{j}(x)}=d_{0}b_{0}(x)+d_{1}b_{1}(x)+d_{2}b_{2}(x)=d_{0}+d_{1}(x+1)+d_{2}(x+1)^{2}$

$u^{h}(1)=\sum _{j=0}^{2}{d_{j}b_{j}(1)}=d_{0}b_{0}(1)+d_{1}b_{1}(1)+d_{2}b_{2}(1)=d_{0}+2d_{1}+4d_{2}=0$

So the relevant equation is

$\displaystyle d_{0}+2d_{1}+4d_{2}=0$

(3.1.6)

Part 3

Equation 3.1.2 can be used to project the residue onto the basis function. First, the partial differential equation $\displaystyle P(u^{h})$ is defined.

$\displaystyle P(u^{h})={\frac {\partial }{\partial x}}\left[a_{2}(x){\frac {\partial u}{\partial x}}\right]+f(x,t)-{\overline {m}}(x){\frac {\partial ^{s}u}{\partial t^{s}}}$

(3.1.7)

Substituting the given conditions into this equation and recalling equation 3.1.1,

$\displaystyle P(u^{h})={\frac {\partial }{\partial x}}\left\{2\left[d_{1}+2d_{2}(x+1)\right]\right\}+3$

(3.1.8)

Differentiating with respect to $\displaystyle x$ and substituting known values,

$\displaystyle P(u^{h})=\left[4d_{2}\right]+3$

(3.1.9)

Substituting into 3.1.2

$\displaystyle <b_{k},P(u^{h})>=\int _{0}^{1}b_{k}(x)P(u^{h})dx=\int _{0}^{1}b_{k}(x)\left[4d_{2}+3\right]dx=0$

(3.1.10)

Since the equation is valid for all $\displaystyle b_{k}(x),\ b_{0}(x)=(x+1)$ is utilized.

$\displaystyle <b_{k},P(u^{h})>=\int _{0}^{1}(x+1)\left[4d_{2}+3\right]dx=\int _{0}^{1}[(3+4d_{2})x+4d_{2}+3]dx=(3+4d_{2}){\frac {1^{2}}{2}}+(4d_{2}+3)1=6d_{2}+{\frac {9}{2}}=0$

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Expected width > 0."): {\displaystyle \displaystyle }

Simplifying, the relevant equation becomes

$\displaystyle 4d_{2}=-3$

(3.1.11)

Part 4

The coefficient ($K$) and constant ($f$) matrices are constructed from the equations determined aboved.

Notes:

• The first row of the matrix $\mathbf {K}$ is based on equation 3.1.1 and the essential boundary condition, 3.1.4, when $\displaystyle x=1$.
  
• The second row of the matrix $\mathbf {K}$ is from equations 3.1.5 and 3.1.6 (for the general case).
  
• The third row of the matrix $\mathbf {K}$ is from equation 3.1.11.
  

${\begin{bmatrix}0&1&2\\1&2&4\\0&0&4\end{bmatrix}}{\begin{bmatrix}d_{0}\\d_{1}\\d_{2}\end{bmatrix}}={\begin{bmatrix}-4\\0\\-3\end{bmatrix}}$ Note $\displaystyle \mathbf {K}$ is not Symmetric.

(3.1.12)

Part 5

To solve for the matrix $\mathbf {d}$, first recognize that $\mathbf {K} \cdot \mathbf {d} =\mathbf {F}$ can be rewritten as $\mathbf {d} =\mathbf {K} ^{-1}\cdot \mathbf {F}$.

${\begin{bmatrix}d_{0}\\d_{1}\\d_{2}\end{bmatrix}}={\begin{bmatrix}0&1&2\\1&2&4\\0&0&4\end{bmatrix}}^{-1}{\begin{bmatrix}-4\\0\\-3\end{bmatrix}}$ ${\begin{bmatrix}d_{0}\\d_{1}\\d_{2}\end{bmatrix}}={\begin{bmatrix}-2&1&0\\1&0&-0.5\\0&0&0.25\end{bmatrix}}{\begin{bmatrix}-4\\0\\-3\end{bmatrix}}$ $\displaystyle {\underline {d}}_{n=2}={\begin{bmatrix}8\\-2.5\\-0.75\end{bmatrix}}$ (3.1.13)

The inverse matrix and d coefficients were found using the following matlab script:

Matlab script

function []=fem3N2(); % enter boundary conditions f matrix fmat(1,1)=0; fmat(2,1)=-4; N=2; % create rows of k matrix for boundary conditions matrix(1,1)=1; matrix(2,1)=0; n=1; for kk=2:2:N matrix(1,kk)=cos(n*1); matrix(1,kk+1)=sin(n*1); matrix(2,kk)=-1*n*sin(n*0); matrix(2,kk+1)=n*cos(n*0); n=n+1; end %Enter in bo innerproducts for first column of matrix n=-1; % create additional rows of K matrix for given basis for jj=3:N+1 matrix(jj,1)=0; l=1; m=1; n1=1; n=n+1; for kk=2:2:N % loops over number of terms desired (N) if jj==3 % for b0 or when b is equal to 1 matrix(jj,kk)=2*((-1*n1*sin(n1*1))+(n1*sin(n1*0))); matrix(jj,kk+1)=2*((n1*cos(n1*1))-(n1*cos(n1*0))); n1=n1+1; elseif mod(jj,2)==0 % for when first b of inner product is even if n==n1 % for when incrementing terms of the inner % product are equal matrix(jj,kk)=(-2*(n1^2))*((((2*n1*1)+sin((2*n1*1)))/(4*n1))-(((2*n1*0)+(sin(2*n1*0)))/(4*n1))); matrix(jj,kk+1)=(-2*(n1^2))*(((-1*(cos(n1*1)^2)/(2*n1)))-((-1*(cos(n1*0)^2)/(2*n1)))); n1=n1+1; else % for when incrementing terms of the inner % product are not equal matrix(jj,kk)=(-2*(n1^2))*((((n*sin(n*1)*cos(n1*1))-(n1*cos(n*1)*sin(n1*1)))/((n^2)-(n1^2)))-(((n*sin(n*0)*cos(n1*0))-(n1*cos(n*0)*sin(n1*0)))/((n^2)-(n1^2)))); matrix(jj,kk+1)=(-2*(n1^2))*((((n*sin(n*1)*sin(n1*1))+(n1*cos(n*1)* ... cos(n1*1)))/ ... ((n^2)-(n1^2)))-(((n*sin(n*0)*sin(n1*0))+(n1*cos(n*0)*cos(n1*0)))/((n^2)-(n1^2)))); n1=n1+1; end else % when first b of inner product is odd if n==n1 % for when incrementing terms of the inner % product are equal a=1; matrix(jj,kk)=(-2*(n1^2))*(((-1*(cos(n1*1)^2)/(2*n1)))-((-1*(cos(n1*0)^2)/(2*n1)))); matrix(jj,kk+1)=(-2*(n1^2))*(((1/2)-((sin(2*n1*1))/(4*n1)))-((0/2)-((sin(2*n1*0))/(4*n1)))); n1=n1+1; else % for when incrementing terms of the inner % product are not equal matrix(jj,kk)=(-2*(n1^2))*((-1*((n1*sin(n*1)*sin(n1*1))+(n*cos(n*1)*cos(n1*1)))/((n^2)-(n1^2)))+(((n1*sin(n*0)*sin(n1*0))+(n*cos(n*0)*cos(n1*0)))/((n^2)-(n1^2)))); matrix(jj,kk+1)=(-2*(n1^2))*((((n1*sin(n*1)*cos(n1*1))-(n*cos(n*1)*sin(n1*1)))/((n^2)-(n1^2)))-(((n1*sin(n*0)*cos(n1*0))-(n*cos(n*0)*sin(n1*0)))/((n^2)-(n1^2)))); n1=n1+1; end end end end n=1; % <b0,f> fmat(3,1)=-3; % <bi,f> if N > 3 for m=4:2:N fmat(m,1)=-3*(((1/n)*sin(n*1))-((1/n)*sin(n*0))); fmat(m+1,1)=-3*(((-1/n)*cos(n*1))+((1/n)*cos(n*0))); n=n+1; end end % Calculate d coefficients dmat=(matrix^-1)*fmat;

Parts 6 & 7

The approximate function, $\displaystyle u_{n}^{h}(x)$, can now be written using the determined coefficients in equation 3.1.13. The exact solution, was obtained and is displayed here:

$\displaystyle u(x)={\frac {-3}{4}}x^{2}-4x+{\frac {19}{4}}$

The same procedure is repeated for $n=4$ and $n=6$ in equation 3.1.1.

As the procedure is the same shown above, a Matlab script was written to enforce the boundary conditions, perform the inner products and integration, solve the system of equations, and plot the results.

Matlab script

% enter boundary conditions f matrix fmat(1,1)=0; fmat(2,1)=-4; N=6; k=1; % create rows of k matrix for voundary conditions for kk=1:N+1 %essential matrix(1,kk)=(1+k)^(kk-1); %natural matrix(2,kk)=(kk-1)*(k)^(kk-2); end matrix(2,1)=0; x=sym('x');d0=sym('d0');d1=sym('d1');d2=sym('d2'); d3=sym('d3');d4=sym('d4');d5=sym('d5');d6=sym('d6'); % create additional rows of K matrix for given basis for kk=1:N+1 % loops over number of terms desired (N) %matrix(jj+2,kk)=2*(kk-1)*(kk-2)*(x+k)^(kk-3); Puh(1,kk)=2*(kk-1)*(kk-2)*(x+k)^(kk-3); end for kk=1:N-1 bkpu=(Puh(1)*d0+Puh(2)*d1+Puh(3)*d2+Puh(4)*d3+Puh(5)*d4+Puh(6)*d5+Puh(7)*d6+3)*(x+k)^(kk-1); eqn=int(bkpu,x,0,1) end %**evaluate equations and append them to matrix and fmat matrix(3,:)=[0 0 4 6 8 10 12];fmat(3)=-3; matrix(4,:)=[0 0 2 4 6 8 10];fmat(4)=-3/2; matrix(5,:)=[0 0 4/3 3 24/5 20/3 60/7];fmat(5)=-1; matrix(6,:)=[0 0 1 12/5 4 40/7 15/2];fmat(6)=-3/4; matrix(7,:)=[0 0 4/5 2 24/7 5 20/3];fmat(7)=-3/5; % Calculate d coefficients dmat=matrix^-1*fmat % Y values for N=6 xp1=zeros(21,1); n2error=0; n4error=0; n6error=0; n=1; for ii=0:.05:1 for kk=1:N+1 xp1(n,1)=dmat(kk)*(ii+k)^(kk-1)+xp1(n,1); if ii==0 n2error=dmat(kk)*(.55+k)^(kk-1)+n2error; end end n=n+1; end % Y values for N=4 xp2=zeros(21,1); n=1; N=4; for ii=0:.05:1 for kk=1:N+1 xp2(n,1)=dmat(kk)*(ii+k)^(kk-1)+xp2(n,1); if ii==0 n4error=dmat(kk)*(.55+k)^(kk-1)+n2error; end end n=n+1; end % Y values for N=2 xp3=zeros(21,1); n=1; N=2; for ii=0:.05:1 for kk=1:N+1 xp3(n,1)=dmat(kk)*(ii+k)^(kk-1)+xp3(n,1); u(n)=((-3/4)*ii^2)-(4*ii)+19/4; if ii==0 n6error=dmat(kk)*(.55+k)^(kk-1)+n6error; end end n=n+1; end x=0:.05:1; uexact=((-3/4)*.5^2)-(4*.5)+19/4; error(1)=uexact-n2error; n(1)=2; error(2)=uexact-n4error; n(2)=4; error(3)=uexact-n6error; n(3)=6; figure (1) plot(x,xp3,'-b',x,xp2,'--r',x,xp1,':k',x,u,'+') xlabel('x') ylabel('u') title('(x+1)^j basis functions') legend('N=2','N=4','N=6','uact') figure (2) plot(n,error) xlabel('n') ylabel('error')

Figure 3.1.1 shows the approximate solutions for $n=2,4,6$ plotted against the exact solution in the domain $\Omega$. It is clear that $n=2$ provides a good approximation. This is likely due to the fact that the approximate and exact solutions are both second order polynomial functions.

Plot of u and u^h versus x

Part 8

The previous Matlab code also calculates the error in the estimate for u or the following equation:

$\displaystyle e_{n}(x)=u(x)-u^{h}(x)$

(3.1.14)

This equation was evaluated at .5 and plotted vs n giving the following figure:

Plot of u-u^h versus n

Author

Johnathan Whittaker Bullard

Reviewers / Team Members

Reviewers Table

Name	Reviewed (y/n)
Johnathan Whittaker Bullard	Y
Philip Flater	Y
Brandon Hua	Y
Jiang Jin	Y
Braden Snook	Y
Srilalithkumar Swaminathan	Y
Zongyi Yang

Problem 3.2: Fish and Belytschko 2.1

Problem Description

Consider figure 2.16 from "A First Course in Finite Elements" by Fish and Belytschko p. 37 pb.2.1

Part a: Number the elements and nodes.

Part b: Assemble the global stiffness and force matrix.

Part c: Partition the system and solve for the nodal displacements.

Part d: Compute the reaction forces.

Solution

Part a

The elements are numbered in $\displaystyle {\color {red}red}$ and the nodes are labled in $\displaystyle {\color {blue}blue}$

Part b

Assemble the global stiffness and force matrix.

$\displaystyle {\color {red}Element1}$, $\displaystyle {\color {blue}I=1}$, $\displaystyle {\color {blue}J=4}$

$\displaystyle \mathbf {K} ^{(1)}=3k{\begin{vmatrix}1&-1\\-1&1\end{vmatrix}}\Rightarrow \mathbf {\tilde {K}} ^{(1)}={\begin{vmatrix}3k&0&0&-3k\\0&0&0&0\\0&0&0&0\\-3k&0&0&3k\end{vmatrix}}$

$\displaystyle {\color {red}Element2}$, $\displaystyle {\color {blue}I=1}$, $\displaystyle {\color {blue}J=3}$

$\displaystyle \mathbf {K} ^{(2)}=k{\begin{vmatrix}1&-1\\-1&1\end{vmatrix}}\Rightarrow \mathbf {\tilde {K}} ^{(2)}={\begin{vmatrix}k&0&-k&0\\0&0&0&0\\-k&0&k&0\\0&0&0&0\end{vmatrix}}$

$\displaystyle {\color {red}Element3}$, $\displaystyle {\color {blue}I=3}$, $\displaystyle {\color {blue}J=4}$

$\displaystyle \mathbf {K} ^{(3)}=2k{\begin{vmatrix}1&-1\\-1&1\end{vmatrix}}\Rightarrow \mathbf {\tilde {K}} ^{(3)}={\begin{vmatrix}0&0&0&0\\0&0&0&0\\0&0&2k&-2k\\0&0&-2k&2k\end{vmatrix}}$

$\displaystyle {\color {red}Element4}$, $\displaystyle {\color {blue}I=4}$, $\displaystyle {\color {blue}J=2}$

$\displaystyle \mathbf {K} ^{(4)}=k{\begin{vmatrix}1&-1\\-1&1\end{vmatrix}}\Rightarrow \mathbf {\tilde {K}} ^{(4)}={\begin{vmatrix}0&0&0&0\\0&k&0&-k\\0&0&0&0\\0&-k&0&k\end{vmatrix}}$

Assembled system matrix:

$\displaystyle \mathbf {K} =\sum _{e=1}^{4}\mathbf {\tilde {K}} ^{e}={\begin{vmatrix}4k&0&-k&-3k\\0&k&0&-k\\-k&0&3k&-2k\\-3k&-k&-2k&6k\end{vmatrix}}$

The displacement and force matrices for the system are:

$\displaystyle \mathbf {d} ={\begin{vmatrix}0\\0\\u_{3}\\u_{4}\\\end{vmatrix}},\;\mathbf {f} ={\begin{vmatrix}0\\0\\0\\50\\\end{vmatrix}}N,\;\mathbf {r} ={\begin{vmatrix}r_{1}\\r_{2}\\0\\0\\\end{vmatrix}}\;$

The global system of equations is given by:

$\displaystyle \mathbf {K} =\sum _{e=1}^{4}\mathbf {\tilde {K}} ^{e}={\begin{vmatrix}4k&0&-k&-3k\\0&k&0&-k\\-k&0&3k&-2k\\-3k&-k&-2k&6k\end{vmatrix}}\displaystyle {\begin{vmatrix}0\\0\\u_{3}\\u_{4}\\\end{vmatrix}}={\begin{vmatrix}r_{1}\\r_{2}\\0\\50\\\end{vmatrix}}$

Part c

Partition the system and solve for the nodal displacements.

As the first two displacements are prescribed, we partition after two rows and columns:

$\displaystyle {\begin{vmatrix}\mathbf {K} _{E}&\mathbf {K} _{EF}\\\mathbf {K} _{EF}^{T}&\mathbf {K} _{F}\end{vmatrix}}{\begin{vmatrix}\mathbf {\overline {d}} _{E}\\\mathbf {d} _{F}\end{vmatrix}}={\begin{vmatrix}\mathbf {r} _{E}\\\mathbf {f} _{F}\end{vmatrix}}$

Where

$\displaystyle \mathbf {K} _{E}={\begin{vmatrix}4k&0\\0&k\end{vmatrix}}$, $\displaystyle \mathbf {K} _{F}={\begin{vmatrix}3k&2k\\-2k&6k\end{vmatrix}}$, $\displaystyle \mathbf {K} _{EF}={\begin{vmatrix}-k&-3k\\0&-k\end{vmatrix}}$, $\displaystyle \mathbf {\overline {d}} _{E}={\begin{vmatrix}0\\0\end{vmatrix}}$, $\displaystyle \mathbf {d} _{f}={\begin{vmatrix}u_{3}\\u_{4}\end{vmatrix}}$, $\displaystyle \mathbf {f} _{F}={\begin{vmatrix}0\\50\end{vmatrix}}$, $\displaystyle \mathbf {e} _{E}={\begin{vmatrix}r_{1}\\r_{2}\end{vmatrix}}$

Solving for the nodal displacements:

$\displaystyle {\begin{vmatrix}\mathbf {K} _{F}\\\end{vmatrix}}{\begin{vmatrix}\mathbf {d} _{F}\\\end{vmatrix}}={\begin{vmatrix}\mathbf {f} _{F}\\\end{vmatrix}}$

$\displaystyle {\begin{vmatrix}3k&-2k\\-2k&6k\end{vmatrix}}{\begin{vmatrix}u_{3}\\u_{4}\end{vmatrix}}={\begin{vmatrix}0\\50\end{vmatrix}}$

The systems of equations are:

$\displaystyle 3k\cdot u_{3}-2k\cdot u_{4}=0$

Muliply by -3

$\displaystyle -3(3k\cdot u_{3}-2k\cdot u_{4})=0$

$\displaystyle -9k\cdot u_{3}-6k\cdot u_{4}=0$

and

$\displaystyle -2k\cdot u_{3}-6k\cdot u_{4}=50$

Summing the two equations lets us solve for the displacements as:

$\displaystyle u_{3}={\frac {-50}{7k}},\;u_{4}={\frac {75}{7k}}$

Part d

Compute the reaction forces:

The global system now becomes:

$\displaystyle \mathbf {K} =\sum _{e=1}^{4}\mathbf {\tilde {K}} ^{e}={\begin{vmatrix}4k&0&-k&-3k\\0&k&0&-k\\-k&0&3k&-2k\\-3k&-k&-2k&6k\end{vmatrix}}\displaystyle {\begin{vmatrix}0\\0\\{\frac {-50}{11k}}\\{\frac {75}{11k}}\\\end{vmatrix}}={\begin{vmatrix}r_{1}\\r_{2}\\0\\50\\\end{vmatrix}}$

Solving for the reation forces:

$\displaystyle {\begin{vmatrix}\mathbf {K} _{EF}\\\end{vmatrix}}{\begin{vmatrix}\mathbf {d} _{F}\\\end{vmatrix}}={\begin{vmatrix}\mathbf {r} _{E}\\\end{vmatrix}}$

The equations for the reaction forces are:

$\displaystyle -k\left({\frac {-50}{11k}}\right)-3k\left({\frac {75}{11k}}\right)=r_{1}$

$\displaystyle -k\left({\frac {75}{11k}}\right)=r_{2}$

The reaction forces equal:

$\displaystyle r_{1}={\frac {-175}{11}}N$ and $r_{2}={\frac {75}{11}}N$

Author

Brandonhua 04:48, 15 February 2011 (UTC)

Reviewers / Team Members

Reviewers Table

Name	Reviewed (y/n)
Johnathan Whittaker Bullard
Philip Flater	Y
Brandon Hua	Y
Jiang Jin
Braden Snook	Y
Srilalithkumar Swaminathan	Y
Zongyi Yang

Problem 3.3: Fish and Belytschko 2.3

Given

Looking at the truss structure below where nodes A and B are fixed. A 10 N force in the positive x-direction acts at node C. The joint locations are in meters. The corresponding Young‘s modulus is $E={10^{11}}{\text{Pa}}$ and cross-sectional area for all bars are $A=2\cdot {10^{-2}}{{\text{m}}^{2}}$.

Find

1. Number the elements and nodes.
2. Assemble the global stiffness and force matrix.
3. Partition the system and solve for the nodal displacements.
4. Compute the stresses and reactions.

Solution

1. The elements and nodes numbers are shown in the figure below:

2. Dividing the structure into 4 elements and 4 nodes, and deals with each element starting with element 1: Element 1 is numbered with global nodes 1 and 4. It is positioned at an angle ${\phi ^{\left(1\right)}}={90^{\circ }}$with respect to positive x-axis. Then,

$\displaystyle \cos {90^{\circ }}=0$

(Eq 3.1)

$\displaystyle \sin {90^{\circ }}=1$

(Eq 3.2)

$\displaystyle {l^{\left(1\right)}}=1{\text{m}}$

(Eq 3.3)

$\displaystyle {k^{\left(1\right)}}={\frac {{A^{\left(1\right)}}{E^{\left(1\right)}}}{l^{\left(1\right)}}}=2\times {10^{9}}{\text{N/m}}$

(Eq 3.4)

$\displaystyle {{\mathbf {K} }^{\left(1\right)}}=2\times {10^{9}}\left[{\begin{array}{ccccccccccccccc}0&0&0&0\\0&1&0&{-1}\\0&0&0&0\\0&{-1}&0&1\end{array}}\right]$

(Eq 3.5)

Element 2 is numbered with global nodes 2 and 4. It is positioned at an angle ${\phi ^{\left(2\right)}}={135^{\circ }}$with respect to positive x-axis.

$\displaystyle \cos {135^{\circ }}=-{\frac {1}{\sqrt {2}}}$

(Eq 3.6)

$\displaystyle \sin {135^{\circ }}={\frac {1}{\sqrt {2}}}$

(Eq 3.7)

$\displaystyle {l^{\left(2\right)}}={\sqrt {2}}{\text{m}}$

(Eq 3.8)

$\displaystyle {k^{\left(2\right)}}={\frac {{A^{\left(2\right)}}{E^{\left(2\right)}}}{l^{\left(2\right)}}}={\frac {2}{\sqrt {2}}}\times {10^{9}}{\text{N/m}}$

(Eq 3.9)

$\displaystyle {{\mathbf {K} }^{\left(2\right)}}={\frac {2}{\sqrt {2}}}\times {10^{9}}\left[{\begin{array}{ccccccccccccccc}{\frac {1}{2}}&{-{\frac {1}{2}}}&{-{\frac {1}{2}}}&{\frac {1}{2}}\\{-{\frac {1}{2}}}&{\frac {1}{2}}&{\frac {1}{2}}&{-{\frac {1}{2}}}\\{-{\frac {1}{2}}}&{\frac {1}{2}}&{\frac {1}{2}}&{-{\frac {1}{2}}}\\{\frac {1}{2}}&{-{\frac {1}{2}}}&{-{\frac {1}{2}}}&{\frac {1}{2}}\end{array}}\right]$

(Eq 3.10)

Element 3 is numbered with global nodes 3 and 4. It is positioned at an angle ${\phi ^{\left(3\right)}}={180^{\circ }}$with respect to positive x-axis.

$\displaystyle \cos {180^{\circ }}=-1$

(Eq 3.11)

$\displaystyle \sin {180^{\circ }}=0$

(Eq 3.12)

$\displaystyle {l^{\left(3\right)}}=1{\text{m}}$

(Eq 3.13)

$\displaystyle {k^{\left(3\right)}}={\frac {{A^{\left(3\right)}}{E^{\left(3\right)}}}{l^{\left(3\right)}}}=2\times {10^{9}}{\text{N/m}}$

(Eq 3.14)

$\displaystyle {{\mathbf {K} }^{\left(2\right)}}=2\times {10^{9}}\left[{\begin{array}{ccccccccccccccc}1&0&{-1}&0\\0&0&0&0\\{-1}&0&1&0\\0&0&0&0\end{array}}\right]$

(Eq 3.15)

Element 4 is numbered with global nodes 2 and 3. It is positioned at an angle ${\phi ^{\left(4\right)}}={90^{\circ }}$with respect to positive x-axis.

$\displaystyle \cos {90^{\circ }}=0$

(Eq 3.16)

$\displaystyle \sin {90^{\circ }}=1$

(Eq 3.17)

$\displaystyle {l^{\left(4\right)}}=1{\text{m}}$

(Eq 3.18)

$\displaystyle {k^{\left(4\right)}}={\frac {{A^{\left(4\right)}}{E^{\left(4\right)}}}{l^{\left(4\right)}}}=2\times {10^{9}}{\text{N/m}}$

(Eq 3.19)

$\displaystyle {{\mathbf {K} }^{\left(4\right)}}=2\times {10^{9}}\left[{\begin{array}{ccccccccccccccc}0&0&0&0\\0&1&0&{-1}\\0&0&0&0\\0&{-1}&0&1\end{array}}\right]$

(Eq 3.20)

Assemble the global stiffness matrix

$\displaystyle {\mathbf {K} }=2\times {10^{9}}\left[{\begin{array}{ccccccccccccccc}0&0&0&0&0&0&0&0\\0&1&0&0&0&0&0&{-1}\\0&0&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}\\0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{1+{\frac {1}{2{\sqrt {2}}}}}&0&{-1}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&0&0&0&1&0&{-1}&0\\0&0&0&{-1}&0&1&0&0\\0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}&{-1}&0&{1+{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&{-1}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{1+{\frac {1}{2{\sqrt {2}}}}}\end{array}}\right]$

(Eq 3.21)

The external force and reaction matrix

$\displaystyle {\mathbf {f} }=\left[{\begin{array}{ccccccccccccccc}0\\0\\0\\0\\{10}\\0\\0\\0\end{array}}\right]$

(Eq 3.22)

$\displaystyle {\mathbf {r} }=\left[{\begin{array}{ccccccccccccccc}{r_{1x}}\\{r_{1y}}\\{r_{2x}}\\{r_{2y}}\\0\\0\\0\\0\end{array}}\right]$

(Eq 3.23)

And the displacement matrix

$\displaystyle {\mathbf {d} }=\left[{\begin{array}{ccccccccccccccc}0\\0\\0\\0\\{u_{3x}}\\{u_{3y}}\\{u_{4x}}\\{u_{4y}}\end{array}}\right]$

(Eq 3.24)

3. Global system of equation

$\displaystyle 2\times {10^{9}}\left[{\begin{array}{ccccccccccccccc}0&0&0&0&0&0&0&0\\0&1&0&0&0&0&0&{-1}\\0&0&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}\\0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{1+{\frac {1}{2{\sqrt {2}}}}}&0&{-1}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&0&0&0&1&0&{-1}&0\\0&0&0&{-1}&0&1&0&0\\0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}&{-1}&0&{1+{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&{-1}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{1+{\frac {1}{2{\sqrt {2}}}}}\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}0\\0\\0\\0\\{u_{3x}}\\{u_{3y}}\\{u_{4x}}\\{u_{4y}}\end{array}}\right]={\mathbf {r} }=\left[{\begin{array}{ccccccccccccccc}{r_{1x}}\\{r_{1y}}\\{r_{2x}}\\{r_{2y}}\\{10}\\0\\0\\0\end{array}}\right]$

(Eq 3.25)

Partition the system

$\displaystyle {{\overline {\mathbf {d} }}_{E}}=\left[{\begin{array}{ccccccccccccccc}{{\overline {u}}_{1x}}\\{{\overline {u}}_{1y}}\\{{\overline {u}}_{2x}}\\{{\overline {u}}_{2y}}\end{array}}\right]=\left[{\begin{array}{ccccccccccccccc}0\\0\\0\\0\end{array}}\right],{{\mathbf {d} }_{F}}=\left[{\begin{array}{ccccccccccccccc}{u_{3x}}\\{u_{3y}}\\{u_{4x}}\\{u_{4y}}\end{array}}\right],{{\mathbf {f} }_{F}}=\left[{\begin{array}{ccccccccccccccc}{10}\\0\\0\\0\end{array}}\right],{{\mathbf {K} }_{F}}=\left[{\begin{array}{ccccccccccccccc}1&0&{-1}&0\\0&1&0&0\\{-1}&0&{1+{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{1+{\frac {1}{2{\sqrt {2}}}}}\end{array}}\right]$

$\displaystyle {{\mathbf {r} }_{E}}=\left[{\begin{array}{ccccccccccccccc}{r_{1x}}\\{r_{1y}}\\{r_{2x}}\\{r_{2y}}\end{array}}\right],{{\mathbf {K} }_{EF}}=\left[{\begin{array}{ccccccccccccccc}0&0&0&0\\0&0&0&{-1}\\0&0&{-{\frac {2}{2{\sqrt {2}}}}}&{\frac {2}{2{\sqrt {2}}}}\\0&{-1}&{\frac {2}{2{\sqrt {2}}}}&{-{\frac {2}{2{\sqrt {2}}}}}\end{array}}\right]$

Solve for the nodal displacements

$\displaystyle 2\times {10^{9}}\left[{\begin{array}{ccccccccccccccc}1&0&{-1}&0\\0&1&0&0\\{-1}&0&{1+{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{1+{\frac {1}{2{\sqrt {2}}}}}\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}{u_{3x}}\\{u_{3y}}\\{u_{4x}}\\{u_{4y}}\end{array}}\right]=\left[{\begin{array}{ccccccccccccccc}{10}\\0\\0\\0\end{array}}\right]$

(Eq 3.26)

$\displaystyle \left[{\begin{array}{ccccccccccccccc}{u_{3x}}\\{u_{3y}}\\{u_{4x}}\\{u_{4y}}\end{array}}\right]={\frac {1}{2}}\times {10^{-9}}\left[{\begin{array}{ccccccccccccccc}{48.2843}\\0\\{38.2843}\\{10}\end{array}}\right]$

(Eq 3.27)

4. The reaction matrix is

$\displaystyle {{\mathbf {r} }_{E}}=\left[{\begin{array}{ccccccccccccccc}{r_{1x}}\\{r_{1y}}\\{r_{2x}}\\{r_{2y}}\end{array}}\right]={{\mathbf {K} }_{E}}{{\mathbf {\bar {d}} }_{E}}+{{\mathbf {K} }_{EF}}{{\mathbf {d} }_{F}}=\left[{\begin{array}{ccccccccccccccc}0&0&0&0\\0&0&0&{-1}\\0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}\\0&{-1}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}{48.2843}\\0\\{38.2843}\\{10}\end{array}}\right]=\left[{\begin{array}{ccccccccccccccc}0\\{-10}\\{-10}\\{10}\end{array}}\right]$

(Eq 3.28)

The stresses in the two elements are

$\displaystyle {\sigma ^{e}}={\frac {E^{e}}{l^{e}}}\left[{\begin{array}{ccccccccccccccc}{-\cos {\phi ^{e}}}&{-\sin {\phi ^{e}}}&{\cos {\phi ^{e}}}&{\sin {\phi ^{e}}}\end{array}}\right]{{\mathbf {d} }^{e}}$

(Eq 3.29)

For element 1:

$\displaystyle {{\mathbf {d} }^{\left(1\right)}}=\left[{\begin{array}{ccccccccccccccc}{u_{1x}}\\{u_{1y}}\\{u_{4x}}\\{u_{4y}}\end{array}}\right]={\frac {1}{2}}\times {10^{-9}}\left[{\begin{array}{ccccccccccccccc}0\\0\\{38.2843}\\{10}\end{array}}\right]$

(Eq 3.30)

$\displaystyle {\sigma ^{\left(1\right)}}=\left[{\begin{array}{ccccccccccccccc}0&{-1}&0&1\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}0\\0\\{38.2843}\\{10}\end{array}}\right]{\frac {1}{2\times {{10}^{-2}}}}=500{\text{Pa}}$

(Eq 3.31)

For element 2:

$\displaystyle {{\mathbf {d} }^{\left(2\right)}}=\left[{\begin{array}{ccccccccccccccc}{u_{2x}}\\{u_{2y}}\\{u_{4x}}\\{u_{4y}}\end{array}}\right]={\frac {1}{2}}\times {10^{-9}}\left[{\begin{array}{ccccccccccccccc}0\\0\\{38.2843}\\{10}\end{array}}\right]$

(Eq 3.32)

$\displaystyle {\sigma ^{\left(2\right)}}=\left[{\begin{array}{ccccccccccccccc}{\frac {1}{\sqrt {2}}}&{-{\frac {1}{\sqrt {2}}}}&{-{\frac {1}{\sqrt {2}}}}&{\frac {1}{\sqrt {2}}}\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}0\\0\\{38.2843}\\{10}\end{array}}\right]{\frac {1}{2\times {{10}^{-2}}}}=-1000{\text{Pa}}$

(Eq 3.33)

For element 3:

$\displaystyle {{\mathbf {d} }^{\left(3\right)}}=\left[{\begin{array}{ccccccccccccccc}{u_{3x}}\\{u_{3y}}\\{u_{4x}}\\{u_{4y}}\end{array}}\right]={\frac {1}{2}}\times {10^{-9}}\left[{\begin{array}{ccccccccccccccc}{48.2843}\\0\\{38.2843}\\{10}\end{array}}\right]$

(Eq 3.34)

$\displaystyle {\sigma ^{\left(3\right)}}=\left[{\begin{array}{ccccccccccccccc}1&0&{-1}&0\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}{48.2843}\\0\\{38.2843}\\{10}\end{array}}\right]{\frac {1}{2\times {{10}^{-2}}}}=500{\text{Pa}}$

(Eq 3.35)

For element 4:

$\displaystyle {{\mathbf {d} }^{\left(4\right)}}=\left[{\begin{array}{ccccccccccccccc}{u_{2x}}\\{u_{2y}}\\{u_{3x}}\\{u_{3y}}\end{array}}\right]={\frac {1}{2}}\times {10^{-9}}\left[{\begin{array}{ccccccccccccccc}0\\0\\{48.2843}\\0\end{array}}\right]$

(Eq 3.36)

$\displaystyle {\sigma ^{\left(4\right)}}=\left[{\begin{array}{ccccccccccccccc}0&{-1}&0&1\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}0\\0\\{48.2843}\\0\end{array}}\right]{\frac {1}{2\times {{10}^{-2}}}}=0$

(Eq 3.37)

Author

Jiang Jin

Reviewers / Team Members

Reviewers Table

Name	Reviewed (y/n)
Johnathan Whittaker Bullard
Philip Flater	Y
Brandon Hua	Y
Jiang Jin	Y
Braden Snook	Y
Srilalithkumar Swaminathan	Y
Zongyi Yang

Problem 3.4: Show Asymmetry of Weighted Residual Stiffness Matrix

Problem Description

Prove the asymmetry of the weighted residual stiffness matrix by showing:

$\displaystyle \int \mathbf {b} _{j}(x){\big [}{\frac {d}{dx}}a_{2}(x){\frac {d}{dx}}\mathbf {b} _{i}(x)]dx\neq \int \mathbf {b} _{i}(x){\big [}{\frac {d}{dx}}a_{2}(x){\frac {d}{dx}}\mathbf {b} _{j}(x)]dx$

(3.4.1)

Given

The stiffness matrix is defined as:

$\displaystyle K_{i,j}=<\mathbf {b} _{i},\mathbf {P} (u^{h})>$

(3.4.2)

Where < > is for the inner product. Here u^h is an approximation of u from a summation of basis functions given as:

$\displaystyle u(x)\approx u^{h}(x)=\sum _{j}\mathbf {b} _{j}(x)$

(3.4.3)

Also, P(u) is given by:

$\displaystyle \mathbf {P} (\mathbf {u} (x))={\frac {d}{dx}}a_{2}(x){\frac {d}{dx}}\mathbf {u} (x)$

(3.4.4)

Plugging Eq. 3.4.3 into Eq. 3.4.4 and looking only at one term of the summation gives:

$\displaystyle {\frac {d}{dx}}a_{2}(x){\frac {d}{dx}}\mathbf {b} _{j}(x)=(a_{2}*\mathbf {b} _{j}')'$

(3.4.5)

Which for a different numbered basis function than b_j(x) is written as:

$\displaystyle {\frac {d}{dx}}a_{2}(x){\frac {d}{dx}}\mathbf {b} _{i}(x)=(a_{2}*\mathbf {b} _{i}')'$

(3.4.6)

Plugging 3.4.5 and 3.4.6 into 3.4.2 gives:

$\displaystyle \int \mathbf {b} _{j}(x){\big [}{\frac {d}{dx}}a_{2}(x){\frac {d}{dx}}\mathbf {b} _{i}(x)]dx$

(3.4.7)

$\displaystyle \int \mathbf {b} _{i}(x){\big [}{\frac {d}{dx}}a_{2}(x){\frac {d}{dx}}\mathbf {b} _{j}(x)]dx$

(3.4.8)

Which our are equations for Eq.3.4.1, so we now must prove:

$\displaystyle K_{i,j}\neq K_{j,i}$

(3.4.9)

Solution

First, expanding out the integrands of Eq. 3.4.1:

$\displaystyle \alpha =\mathbf {b} _{i}{\big [}a_{2}'*\mathbf {b} _{j}'+a_{2}*\mathbf {b} _{j}''{\big ]}$

(3.4.10)

$\displaystyle \beta =\mathbf {b} _{j}{\big [}a_{2}'*\mathbf {b} _{i}'+a_{2}*\mathbf {b} _{i}''{\big ]}$

(3.4.11)

So now we're simply proving:

$\displaystyle \alpha \neq \beta$

(3.4.12)

This will be proven using the following basis functions:

$\displaystyle \mathbf {b} _{i}(x)=cos(i*x)$

(3.4.13)

$\displaystyle \mathbf {b} _{j}(x)=cos(j*x)$

(3.4.14)

Where $\displaystyle i\neq j$.

Plugging Eq. 3.4.13 and 3.4.14 into Eq. 3.4.10 gives:

$\displaystyle \alpha =cos(i*x){\big [}a_{2}'(-jsin(j*x))+a_{2}(-j^{2}cos(j*x)){\big ]}$

(3.4.15)

Repeating for Eq. 3.4.11:

$\displaystyle \beta =cos(j*x){\big [}a_{2}'(-isin(i*x))+a_{2}(-i^{2}cos(i*x)){\big ]}$

(3.4.16)

Expanding out Eq. 3.4.15 and 3.4.16:

$\displaystyle \alpha =-ja_{2}'cos(i*x)sin(j*x)-j^{2}a_{2}cos(i*x)cos(j*x)$

(3.4.17)

$\displaystyle \beta =-ia_{2}'cos(j*x)sin(i*x)-i^{2}a_{2}cos(j*x)cos(i*x)$

(3.4.18)

Inspecting Eq. 3.4.17 and Eq. 3.4.18 shows:

$\displaystyle \alpha \neq \beta {\text{ when }}i\neq j$

(3.4.19)

Author

Braden Snook

Reviewers / Team Members

Reviewers Table

Name	Reviewed (y/n)
Johnathan Whittaker Bullard
Philip Flater	Y
Brandon Hua	Y
Jiang Jin	Y
Braden Snook	Y
Srilalithkumar Swaminathan	Y
Zongyi Yang

Problem 3.5: Fish and Belytschko 2.2

Problem Statement

Jacob Fish and Ted Belytschko, " A first Course in Finite Elements",John Wiley & Sons, Ltd, Chapter 2, P37 Problem 2.2.

" Show that the Equivalent stiffness of a spring aligned in the x-direction for the bar of thickness t with a centered square hole in figure is

$\displaystyle \mathbf {k={\frac {5Etab}{\left(a+b\right)l}}}$ (3.5.1) Where E is the Young's modulus and t is the width of the bar." Reference figure in textbook: 2.17 Bar of thickness t with a centered square hole Solution The bar can be sub-divided into 3 elements. The node numbers 1,2,3,4 are assgined. The element have stiffness ${\color {red}k^{\left(1\right)}},{\color {green}k^{\left(2\right)}},{\color {blue}k^{\left(3\right)}}$ . We know that stiffness $\mathbf {k={\frac {AE}{L}}}$ Where A is cross sectional area, E is the Elastic modulus and L is the length of the element For element 1, the stiffness matrix is $\displaystyle K^{(1)}={\begin{bmatrix}{\color {red}k^{(1)}}&-{\color {red}k^{(1)}}\\-{\color {red}k^{(1)}}&{\color {red}k^{(1)}}\end{bmatrix}}$ (3.5.2) Where ${\color {red}k^{(1)}}={\frac {A^{1}E^{1}}{L^{1}}}={\frac {(a\times t)\times E}{\frac {l}{10}}}={\frac {10atE}{l}}$ (3.5.3) For element 2, the stiffness matrix is $\displaystyle K^{(2)}={\begin{bmatrix}{\color {green}k^{(2)}}&-{\color {green}k^{(2)}}\\-{\color {green}k^{(2)}}&{\color {green}k^{(2)}}\end{bmatrix}}$ (3.5.4) Where ${\color {green}k^{(2)}}={\frac {A^{2}E^{2}}{L^{2}}}={\frac {(4b\times t)\times E}{\frac {l}{10}}}={\frac {5btE}{l}}$ (3.5.5) For element 3, the stiffness matrix is $\displaystyle K^{(3)}={\begin{bmatrix}{\color {blue}k^{(3)}}&-{\color {blue}k^{(3)}}\\-{\color {blue}k^{(3)}}&{\color {blue}k^{(3)}}\end{bmatrix}}$ (3.5.6) Where ${\color {blue}k^{(3)}}={\frac {A^{3}E^{3}}{L^{3}}}={\frac {(a\times t)\times E}{\frac {l}{10}}}={\frac {10atE}{l}}$ (3.5.7) Now assembling the stiffness matrices we get, $\displaystyle \mathbf {K_{eq}} ={\begin{bmatrix}{\color {red}k^{(1)}}&0&-{\color {red}k^{(1)}}&0\\0&{\color {blue}k^{(3)}}&0&{\color {blue}k^{(3)}}\\-{\color {red}k^{(1)}}&0&{\color {red}k^{(1)}}+{\color {green}k^{(2)}}&-{\color {green}k^{(2)}}\\0&-{\color {blue}k^{(3)}}&-{\color {green}k^{(2)}}&{\color {green}k^{(2)}}+{\color {blue}k^{(3)}}\end{bmatrix}}$ (3.5.8) $\displaystyle \mathbf {K_{eq}} ={\frac {5Et}{l}}\times {\begin{bmatrix}2a&0&-2a&0\\0&2a&0&-2a\\-2a&0&(2a+b)&-b\\0&-2a&-b&(2a+b)\end{bmatrix}}$ (3.5.9) We Know that $\displaystyle \mathbf {\mathit {K.d=F}}$ Where d is displacement and F is force. Taking node(1) to be fixed end. (3.5.10) $\displaystyle {\frac {5Et}{l}}\times {\begin{bmatrix}2a&0&-2a&0\\0&2a&0&-2a\\-2a&0&(2a+b)&-b\\0&-2a&-b&(2a+b)\end{bmatrix}}\cdot {\begin{bmatrix}0\\d_{2}\\d_{3}\\d_{4}\end{bmatrix}}={\begin{bmatrix}-F\\F\\0\\0\end{bmatrix}}$ (3.5.12) Solving the above we get, $\displaystyle {\frac {5Etab}{(a+b)l}}\cdot d_{2}=F$ (3.5.13) $\displaystyle \Rightarrow k={\frac {5Etab}{\left(a+b\right)l}}$ (3.5.13) Author Srilalithkumar 15:41, 15 February 2011 (UTC) Reviewers / Team Members Name Reviewed (y/n) Johnathan Whittaker Bullard Philip Flater Y Brandon Hua Y Jiang Jin Y Braden Snook Y Srilalithkumar Swaminathan Y Zongyi Yang Problem 3.6: Fish and Belytschko 2.4 Problem Description "A First Course in Finite Elements" by Fish and Belytschko p. 38 Problem 2.4 Given the three-bar structure subjected to the prescribed load ad point C equal to ${10^{3}}N$ as shown. The Young’s modulus is $E={{10}^{11}}$ Pa. The cross-sectional area of the bar BC is ${{10}^{-2}}{{m}^{2}}$ and that of BD and BF is ${{10}^{-2}}{{m}^{2}}$. Note that point D is free to move in the x-direction. Coordinates of joins are given in meters. Construct the global stiffness matrix and load matrix. Partition the matrices and solve for the unknown displacement at Point B and displacement in the x-direction at point D. Find the stressed in the three bars. Find the reaction at the nodes C,D and F. Solution a) First subdivide the structure into elements. There are three elements and two nodes. Shown in Fig. Element 1: $\displaystyle \varphi =135^{\circ },\cos \varphi =-{\frac {1}{\sqrt {2}}},\sin \varphi ={\frac {1}{\sqrt {2}}}$ $\displaystyle {K^{(1)}}={\frac {0.01E}{\sqrt {2}}}\left[{\begin{array}{ccccccccccccccc}{\frac {1}{2}}&{-{\frac {1}{2}}}&{-{\frac {1}{2}}}&{\frac {1}{2}}\\{-{\frac {1}{2}}}&{\frac {1}{2}}&{\frac {1}{2}}&{-{\frac {1}{2}}}\\{-{\frac {1}{2}}}&{\frac {1}{2}}&{\frac {1}{2}}&{-{\frac {1}{2}}}\\{\frac {1}{2}}&{-{\frac {1}{2}}}&{-{\frac {1}{2}}}&{\frac {1}{2}}\end{array}}\right]$ (3.6.1) Element 2: $\displaystyle \varphi =90^{\circ },\cos \varphi =0,\sin \varphi =1$ $\displaystyle {K^{(2)}}={\frac {0.02E}{1}}\left[{\begin{array}{ccccccccccccccc}0&0&0&0\\0&1&0&{-1}\\0&0&0&0\\0&{-1}&0&1\end{array}}\right]$ (3.6.2) Element 3: $\displaystyle \varphi =45^{\circ },\cos \varphi ={\frac {1}{\sqrt {2}}},\sin \varphi ={\frac {1}{\sqrt {2}}}$ $\displaystyle {K^{(3)}}={\frac {0.01E}{\sqrt {2}}}\left[{\begin{array}{ccccccccccccccc}{\frac {1}{2}}&{\frac {1}{2}}&{-{\frac {1}{2}}}&{-{\frac {1}{2}}}\\{\frac {1}{2}}&{\frac {1}{2}}&{-{\frac {1}{2}}}&{-{\frac {1}{2}}}\\{-{\frac {1}{2}}}&{-{\frac {1}{2}}}&{\frac {1}{2}}&{\frac {1}{2}}\\{-{\frac {1}{2}}}&{-{\frac {1}{2}}}&{\frac {1}{2}}&{\frac {1}{2}}\end{array}}\right]$ (3.6.3) The global stiffness matrix is: $\displaystyle K={\frac {0.01E}{1}}\left[{\begin{array}{ccccccccccccccc}{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&0&0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}\\{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}&0&0&0&0&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&0&0&0&0&0&0&0\\0&0&0&2&0&0&0&{-2}\\0&0&0&0&{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&0&0&0&{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}&0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{\sqrt {2}}}&0\\{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&{-2}&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&{2+{\frac {1}{\sqrt {2}}}}\end{array}}\right]$ (3.6.4) And $\displaystyle d={\left[{\begin{array}{ccccccccccccccc}0&0&0&0&{u_{3x}}&0&{u_{4x}}&{u_{4y}}\end{array}}\right]^{T}}$ (3.6.5) $\displaystyle F={\left[{\begin{array}{ccccccccccccccc}0&0&0&0&0&0&{{10}^{3}}&0\end{array}}\right]^{T}}$ (3.6.6) $\displaystyle r={\left[{\begin{array}{ccccccccccccccc}{r_{1x}}&{r_{1y}}&{r_{2x}}&{r_{1y}}&0&{r_{3x}}&0&0\end{array}}\right]^{T}}$ (3.6.7) b) Global system of equations: $\displaystyle {\frac {0.01E}{1}}\left[{\begin{array}{ccccccccccccccc}{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&0&0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}\\{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}&0&0&0&0&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&0&0&0&0&0&0&0\\0&0&0&2&0&0&0&{-2}\\0&0&0&0&{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&0&0&0&{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}&0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{\sqrt {2}}}&0\\{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&{-2}&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&{2+{\frac {1}{\sqrt {2}}}}\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}0\\0\\0\\0\\{u_{3x}}\\0\\{u_{4x}}\\{u_{4y}}\end{array}}\right]=\left[{\begin{array}{ccccccccccccccc}{r_{1x}}\\{r_{1y}}\\{r_{2x}}\\{r_{2y}}\\0\\{r_{3y}}\\{{10}^{3}}\\0\end{array}}\right]$ (3.6.8) Then reduce the global system of equations: The global system is partitioned four rows and four columns: And: $\displaystyle {\frac {0.01E}{1}}\left[{\begin{array}{ccccccccccccccc}{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{\sqrt {2}}}&0\\{-{\frac {1}{2{\sqrt {2}}}}}&{-{\frac {1}{2{\sqrt {2}}}}}&0&{2+{\frac {1}{\sqrt {2}}}}\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}{u_{3x}}\\0\\{u_{4x}}\\{u_{4y}}\end{array}}\right]=\left[{\begin{array}{ccccccccccccccc}0\\{r_{3y}}\\{{10}^{3}}\\0\end{array}}\right]$ (3.6.9) So: $\displaystyle \left[{\begin{array}{ccccccccccccccc}{u_{3x}}\\0\\{u_{4x}}\\{u_{4y}}\end{array}}\right]=\left[{\begin{array}{ccccccccccccccc}{(1+2{\sqrt {2}})\times {{10}^{-6}}}\\0\\{({\frac {1}{2}}+2{\sqrt {2}})\times {{10}^{-6}}}\\{0.5\times {{10}^{-6}}}\end{array}}\right]=\left[{\begin{array}{ccccccccccccccc}{3.828\times {{10}^{-6}}}\\0\\{3.328\times {{10}^{-6}}}\\{0.5\times {{10}^{-6}}}\end{array}}\right]$ (3.6.10) $\displaystyle {r_{3y}}=0$ (3.6.11) To sum up, the displacement in X-direction of B is ${{10}^{-6}}$ mm, in Y-direction is ${{10}^{-6}}$ mm. The displacement in X-direction of C is ${{10}^{-6}}$ mm. c) Element 1: $\displaystyle \varphi =135^{\circ },\cos \varphi =-{\frac {1}{\sqrt {2}}},\sin \varphi ={\frac {1}{\sqrt {2}}}$ (3.6.12) $\displaystyle {\sigma _{1}}={\frac {E}{\sqrt {2}}}\left[{\begin{array}{ccccccccccccccc}{\frac {1}{\sqrt {2}}}&{-{\frac {1}{\sqrt {2}}}}&{-{\frac {1}{\sqrt {2}}}}&{\frac {1}{\sqrt {2}}}\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}{3.328\times {{10}^{-6}}}\\{0.5\times {{10}^{-6}}}\\0\\0\end{array}}\right]=1.414\times {10^{5}}Pa$ (3.6.13) Element 2: $\displaystyle \varphi =90^{\circ },\cos \varphi =0,\sin \varphi =1$ (3.6.14) $\displaystyle {\sigma _{2}}={\frac {E}{1}}\left[{\begin{array}{ccccccccccccccc}0&{-1}&0&1\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}{3.328\times {{10}^{-6}}}\\{0.5\times {{10}^{-6}}}\\0\\0\end{array}}\right]=-0.5\times {10^{5}}Pa$ (3.6.15) Element 3: $\displaystyle \varphi =45^{\circ },\cos \varphi ={\frac {1}{\sqrt {2}}},\sin \varphi ={\frac {1}{\sqrt {2}}}$ (3.6.16) $\displaystyle {\sigma _{3}}={\frac {E}{\sqrt {2}}}\left[{\begin{array}{ccccccccccccccc}{-{\frac {1}{\sqrt {2}}}}&{-{\frac {1}{\sqrt {2}}}}&{\frac {1}{\sqrt {2}}}&{\frac {1}{\sqrt {2}}}\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}{3.328\times {{10}^{-6}}}\\{0.5\times {{10}^{-6}}}\\{3.828\times {{10}^{-6}}}\\0\end{array}}\right]=0Pa$ (3.6.17) d) Considering the Global system of equations, then: $\displaystyle {\frac {0.01E}{1}}\left[{\begin{array}{ccccccccccccccc}0&0&{-{\frac {1}{2{\sqrt {2}}}}}&{\frac {1}{2{\sqrt {2}}}}\\0&0&{\frac {1}{2{\sqrt {2}}}}&{-{\frac {1}{2{\sqrt {2}}}}}\\0&0&0&0\\0&0&0&{-2}\end{array}}\right]\left[{\begin{array}{ccccccccccccccc}{u_{3x}}\\0\\{u_{4x}}\\{u_{4y}}\end{array}}\right]=\left[{\begin{array}{ccccccccccccccc}{r_{1x}}\\{r_{1y}}\\{r_{2x}}\\{r_{2y}}\end{array}}\right]$ (3.6.18) So $\displaystyle \left[{\begin{array}{ccccccccccccccc}{r_{1x}}\\{r_{1y}}\\{r_{2x}}\\{r_{2y}}\end{array}}\right]=\left[{\begin{array}{ccccccccccccccc}{-1\times {{10}^{3}}}\\{1\times {{10}^{3}}}\\0\\{-1\times {{10}^{3}}}\end{array}}\right]$ (3.6.19) $\displaystyle {F_{Dx}}={F_{Dy}}=0$ (3.6.20) $\displaystyle {F_{Cx}}={r_{2x}}=0,{F_{Cy}}={r_{2y}}=-1\times {10^{3}}N$ (3.6.21) $\displaystyle {F_{Fx}}={r_{1x}}-1\times {10^{3}},{F_{Fy}}={r_{1y}}=1\times {10^{3}}N$ (3.6.22) Author Zongyi Yang Reviewers / Team Members Name Reviewed (y/n) Johnathan Whittaker Bullard Philip Flater Y Brandon Hua Y Jiang Jin Y Braden Snook Y Srilalithkumar Swaminathan Y Zongyi Yang Problem 3.7: Solving General 1-D Model with Weighted Residual Form and Polynomial Basis Functions Problem Description Repeat problem 3.1 with $k=0$ Solution For the case when $\displaystyle n=2$ Parts 1 & 2 $\left\{b_{j}(x);j=0,1,...,n\right\}$ where $\displaystyle b_{j}(x)=(x+k)^{j}$ $\Rightarrow \quad b_{0}(x)=(x+k)^{0}=1,b_{1}(x)=(x+k)^{1}=(x+k),\ and\ b_{2}(x)=(x+k)^{2}$. The natural boundary condition can be implemented by differentiating $\displaystyle u^{h}(x)$ with respect to $\displaystyle x$ and taking its value at x=0. ${\frac {d{u}^{h}}{dx}}(x)=\sum _{j=0}^{2}{d_{j}b_{j}'(x)}=d_{0}b_{0}'(x)+d_{1}b_{1}'(x)+d_{2}b_{2}'(x)=d_{1}+2d_{2}(x)$ ${\frac {d{u}^{h}}{dx}}(0)=\sum _{j=0}^{2}{d_{j}b_{j}'(0)}=d_{0}b_{0}'(0)+d_{1}b_{1}'(0)+d_{2}b_{2}'(0)=d_{1}=-4$ So the relevant equation is $\displaystyle d_{1}=-4$ The essential boundary condition, $\displaystyle u^{h}(1)=0$, is implemented as follows: $u^{h}(x)=\sum _{j=0}^{2}{d_{j}b_{j}(x)}=d_{0}b_{0}(x)+d_{1}b_{1}(x)+d_{2}b_{2}(x)=d_{0}+d_{1}(x)+d_{2}(x)^{2}$ $u^{h}(1)=\sum _{j=0}^{2}{d_{j}b_{j}(1)}=d_{0}b_{0}(1)+d_{1}b_{1}(1)+d_{2}b_{2}(1)=d_{0}+d_{1}+d_{2}=0$ So the relevant equation is $\displaystyle d_{0}+d_{1}+d_{2}=0$ Part 3 Equation 3.1.2 can be used to project the residue onto the basis function. First, the partial differential equation $\displaystyle P(u^{h})$ is defined. $\displaystyle P(u^{h})={\frac {\partial }{\partial x}}\left[a_{2}(x){\frac {\partial u}{\partial x}}\right]+f(x,t)-{\overline {m}}(x){\frac {\partial ^{s}u}{\partial t^{s}}}$ Substituting the given conditions into this equation and recalling equation 3.1.1, $\displaystyle P(u^{h})={\frac {\partial }{\partial x}}\left\{2\left[d_{1}+2d_{2}(x)\right]\right\}+3$ Differentiating with respect to $\displaystyle x$ and substituting known values, $\displaystyle P(u^{h})=\left[4d_{2}\right]+3$ Substituting into 3.1.2 $\displaystyle <b_{k},P(u^{h})>=\int _{0}^{1}b_{k}(x)P(u^{h})dx=\int _{0}^{1}b_{k}(x)\left[4d_{2}+3\right]dx=0$ Since the equation is valid for all $\displaystyle b_{k}(x),\ b_{1}(x)=(x+1)$ is utilized. $\displaystyle <b_{k},P(u^{h})>=\int _{0}^{1}(x+1)\left[4d_{2}+3\right]dx=\int _{0}^{1}[(3+4d_{2})x+4d_{2}+3]dx=(3+4d_{2}){\frac {1^{2}}{2}}+(4d_{2}+3)1=6d_{2}+{\frac {9}{2}}=0$ Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Expected width > 0."): {\displaystyle \displaystyle } Simplifying, the relevant equation becomes $\displaystyle 4d_{2}=-3$ Part 4 The coefficient ($K$) and constant ($f$) matrices are constructed from the equations determined aboved. Notes: The first row of the matrix $\mathbf {K}$ is based on equation 3.1.1 and the essential boundary condition, 3.1.4, when $\displaystyle x=1$. The second row of the matrix $\mathbf {K}$ is from equations 3.1.5 and 3.1.6 (for the general case). The third row of the matrix $\mathbf {K}$ is from equation 3.1.11. ${\begin{bmatrix}0&1&0\\1&1&1\\0&0&4\end{bmatrix}}{\begin{bmatrix}d_{0}\\d_{1}\\d_{2}\end{bmatrix}}={\begin{bmatrix}-4\\0\\-3\end{bmatrix}}$ Note $\displaystyle \mathbf {K}$is not symmetric Part 5 To solve for the matrix $\mathbf {d}$, first recognize that $\mathbf {K} \cdot \mathbf {d} =\mathbf {F}$ can be rewritten as $\mathbf {d} =\mathbf {K} ^{-1}\cdot \mathbf {F}$. ${\begin{bmatrix}d_{0}\\d_{1}\\d_{2}\end{bmatrix}}={\begin{bmatrix}0&1&0\\1&1&1\\0&0&4\end{bmatrix}}^{-1}{\begin{bmatrix}-4\\0\\-3\end{bmatrix}}$ ${\begin{bmatrix}d_{0}\\d_{1}\\d_{2}\end{bmatrix}}={\begin{bmatrix}-1&1&-0.25\\1&0&0\\0&0&0.25\end{bmatrix}}{\begin{bmatrix}-4\\0\\-3\end{bmatrix}}$ $\displaystyle {\underline {d}}_{n=2}={\begin{bmatrix}4.75\\-4\\-0.75\end{bmatrix}}$ Parts 6 & 7 The approximate function, $\displaystyle u_{n}^{h}(x)$, can now be written using the determined coefficients. The exact solution, was previously obtained and is displayed here: $\displaystyle u(x)={\frac {-3}{4}}x^{2}-4x+{\frac {19}{4}}$ The same procedure is repeated for $n=4$ and $n=6$ in equation 3.1.1. As the procedure is the same shown above, a Matlab script was written to enforce the boundary conditions, perform the inner products and integration, solve the system of equations, and plot the results. Matlab script % enter boundary conditions f matrix fmat(1,1)=0; fmat(2,1)=-4; N=6; k=0; % create rows of k matrix for voundary conditions for kk=1:N+1 %essential matrix(1,kk)=(1+k)^(kk-1); %natural matrix(2,kk)=(kk-1)*(k)^(kk-2); end matrix(2,1)=0; x=sym('x');d0=sym('d0');d1=sym('d1');d2=sym('d2'); d3=sym('d3');d4=sym('d4');d5=sym('d5');d6=sym('d6'); % create additional rows of K matrix for given basis for kk=1:N+1 % loops over number of terms desired (N) %matrix(jj+2,kk)=2*(kk-1)*(kk-2)*(x+k)^(kk-3); Puh(1,kk)=2*(kk-1)*(kk-2)*(x+k)^(kk-3); end for kk=1:N-1 bkpu=(Puh(1)*d0+Puh(2)*d1+Puh(3)*d2+Puh(4)*d3+Puh(5)*d4+Puh(6)*d5+Puh(7)*d6+3)*(x+k)^(kk-1); eqn=int(bkpu,x,0,1) end %**evaluate equations and append them to matrix and fmat matrix(3,:)=[0 0 4 6 8 10 12];fmat(3)=-3; matrix(4,:)=[0 0 2 4 6 8 10];fmat(4)=-3/2; matrix(5,:)=[0 0 4/3 3 24/5 20/3 60/7];fmat(5)=-1; matrix(6,:)=[0 0 1 12/5 4 40/7 15/2];fmat(6)=-3/4; matrix(7,:)=[0 0 4/5 2 24/7 5 20/3];fmat(7)=-3/5; % Calculate d coefficients dmat=matrix^-1*fmat % Y values for N=6 xp1=zeros(21,1); n2error=0; n4error=0; n6error=0; n=1; for ii=0:.05:1 for kk=1:N+1 xp1(n,1)=dmat(kk)*(ii+k)^(kk-1)+xp1(n,1); if ii==0 n2error=dmat(kk)*(.55+k)^(kk-1)+n2error; end end n=n+1; end % Y values for N=4 xp2=zeros(21,1); n=1; N=4; for ii=0:.05:1 for kk=1:N+1 xp2(n,1)=dmat(kk)*(ii+k)^(kk-1)+xp2(n,1); if ii==0 n4error=dmat(kk)*(.55+k)^(kk-1)+n2error; end end n=n+1; end % Y values for N=2 xp3=zeros(21,1); n=1; N=2; for ii=0:.05:1 for kk=1:N+1 xp3(n,1)=dmat(kk)*(ii+k)^(kk-1)+xp3(n,1); u(n)=((-3/4)*ii^2)-(4*ii)+19/4; if ii==0 n6error=dmat(kk)*(.55+k)^(kk-1)+n6error; end end n=n+1; end x=0:.05:1; uexact=((-3/4)*.5^2)-(4*.5)+19/4; error(1)=uexact-n2error; n(1)=2; error(2)=uexact-n4error; n(2)=4; error(3)=uexact-n6error; n(3)=6; figure (1) plot(x,xp3,'-b',x,xp2,'--r',x,xp1,':k',x,u,'+') xlabel('x') ylabel('u') title('(x+1)^j basis functions') legend('N=2','N=4','N=6','uact') figure (2) plot(n,error) xlabel('n') ylabel('error') Figure 3.1.1 shows the approximate solutions for $n=2,4,6$ plotted against the exact solution in the domain $\Omega$. It is clear that $n=2$ provides a good approximation. This is likely due to the fact that the approximate and exact solutions are both second order polynomial functions. Plot of u and uh vs x Part 8 The previous Matlab code also calculates the error in the estimate for u or the following equation: $\displaystyle e_{n}(x)=u(x)-u^{h}(x)$ This equation was evaluated at .5 and plotted vs n giving the following figure: Plot of u - uh vs n Author Johnathan Whittaker Bullard Reviewers / Team Members Reviewers Table Name Reviewed (y/n) Johnathan Whittaker Bullard Y Philip Flater Y Brandon Hua Y Jiang Jin Y Braden Snook Srilalithkumar Swaminathan Y Zongyi Yang Whitaker Bullard/Braden Snook Problem 3.8: Fish and Belytschko 3.1 Problem Description [From Fish & Belytschko, Chapter 3, Problem 3.1] Show that the weak form of ${\frac {d}{dx}}\left(AE{\frac {du}{dx}}\right)+2x=0\qquad {\text{on}}\qquad 1<x<3,$ $\displaystyle \sigma (1)=\left(E{\frac {du}{dx}}\right)_{x=1}=0.1,$ $\displaystyle u(3)=0.001$ is given by $\int _{1}^{3}{\frac {dw}{dx}}AE{\frac {du}{dx}}dx=-0.1(wA)_{x=1}+\int _{1}^{3}2xw\ dx\qquad \forall w\ {\text{ with }}\ w(3)=0.$ Solution The strong form of this problem can be directly compared to the series of equations for the case of one-dimensional stress analysis, and the derivation is therefore adapted from the text [F&B, Section 3.5] First, the governing equation and traction boundary condition are multiplied by an arbitrary (or weight) function, $w(x)$, and integrated from 1 to 3. $\int _{1}^{3}w\left[{\frac {d}{dx}}\left(AE{\frac {du}{dx}}\right)+2x\right]dx=0\qquad \forall w,$ (3.8.1) $\left(wA\left(E{\frac {du}{dx}}-0.1\right)\right)_{x=1}=0\qquad \forall w.$ (3.8.2) In its current form, equation 3.8.1 is twice differentiable and the stiffness matrix is not symmetric. Therefore it will need to be reduced to first derivatives only and hence a symmetric stiffness matrix. It can be rewritten as $\int _{1}^{3}w{\frac {d}{dx}}\left(AE{\frac {du}{dx}}\right)dx+\int _{1}^{3}w(2x)dx=0\qquad \forall w.$ (3.8.3) Using integration by parts we obtain $\left(wAE{\frac {du}{dx}}\right)|_{1}^{3}-\int _{1}^{3}{\frac {dw}{dx}}AE{\frac {du}{dx}}dx+\int _{1}^{3}w(2x)dx=0\qquad \forall w\ {\text{with}}\ w(3)=0.$ (3.8.4) It is important to note that in the derivation of the weak form, the weight functions and trial solutions are constructed so that $w(3)=0{\text{ and }}u={\bar {u}}{\text{ on }}\Gamma _{u}$, respectively. $\displaystyle \Gamma _{u}$ is defined as the boundary where the displacements are prescribed (x=3 in this case); $\displaystyle \Gamma _{t}$ is the boundary where the traction is prescribed (x=1). Recalling Hooke's Law and the definition of strain equation 3.8.4 is rewritten as $(wA\sigma )_{x=3}-(wA\sigma )_{x=1}-\int _{1}^{3}{\frac {dw}{dx}}AE{\frac {du}{dx}}dx+\int _{1}^{3}w(2x)dx=0\qquad \forall w\ {\text{with}}\ w(3)=0.$ (3.8.5) The original traction boundary can be applied to the second term. Recalling that $\displaystyle w(3)=0$ removes the first term. Thus we are left with the weak form in the case of one-dimensional stress analysis. $\int _{1}^{3}{\frac {dw}{dx}}AE{\frac {du}{dx}}dx=-0.1(wA)_{x=1}+\int _{1}^{3}w(2x)dx=0\qquad \forall w\ {\text{with}}\ w(3)=0.$ (3.8.6) Author Philip Flater Reviewers / Team Members Name Reviewed (y/n) Johnathan Whittaker Bullard Philip Flater Y Brandon Hua Y Jiang Jin Y Braden Snook Y Srilalithkumar Swaminathan Y Zongyi Yang Problem 3.9: Fish and Belytschko 3.3 Problem Description "A First Course in Finite Elements" by Fish and Belytschko p. 72, pb.3.3 Consider a trial (candidate) solution of the form $\displaystyle u(x)=\alpha _{0}+\alpha _{1}(x-3)$ and a weight function of the same form. Obtain a solution to the weak form in pb.3.1 from Fish and Belytschko. Check the equilibrium equation in the strong from in pb.3.1 from Fish and Belytschko; is it satisfied? Check the natural boundary condition; is it satisfied? Solution The candidate solution: $\displaystyle u(x)=\alpha _{0}+\alpha _{1}(x-3)$ (3.9.1) The weight solution: $\displaystyle w(x)=\beta _{0}+\beta _{1}(x-3)$ (3.9.2) From pb.3.1 from Fish and Belytschko (eqn 3.8.6 from this hw submission): $\displaystyle \int _{1}^{3}{\frac {dw}{dx}}AE{\frac {du}{dx}}dx=-0.1(wA)_{x=1}+\int _{1}^{3}w(2x)dx=0\qquad \forall w\ {\text{with}}\ w(3)=0$ (3.9.3) The weight function $\displaystyle w(3)=0$ $\displaystyle w(3)=0=\beta _{0}+\beta _{1}(3-3)$ $\displaystyle \beta _{0}=0$ (3.9.4) From pb.3.1 from Fish and Belytschko: $\displaystyle \alpha _{0}=.001$ (3.9.5) $\displaystyle u(3)=\alpha _{0}+\alpha _{1}(3-3)=.001$ (3.9.6) $\displaystyle \alpha _{0}=.001$ (3.9.6) The candidate solution and the weight function is now: $\displaystyle u(x)=.001+\alpha _{1}(x-3)$ (3.9.7) $\displaystyle w(x)=\beta _{1}(x-3)$ (3.9.8) Recall that the weak form is $\displaystyle \int _{1}^{3}{\frac {dw}{dx}}AE{\frac {du}{dx}}dx=-0.1(wA)_{x=1}+\int _{1}^{3}2xw\ dx\qquad \forall w\ {\text{ with }}\ w(3)=0$ (3.9.9) From (3.9.7) and (3.9.8): $\displaystyle {\frac {du(x)}{dx}}=\alpha _{1}$ (3.9.10) $\displaystyle {\frac {dw(x)}{dx}}=\beta _{1}$ (3.9.11) Plug (3.9.7), (3.9.8), (3.9.10) and (3.9.11) into (3.9.9): $\displaystyle \int _{1}^{3}\beta _{1}AE\alpha _{1}dx=-0.1(\beta _{1}(x-3)A)_{x=1}+\int _{1}^{3}2x\beta _{1}(x-3)\ dx\qquad \forall w\ {\text{ with }}\ w(3)=0$ (3.9.12) Integrage (3.9.12) $\displaystyle 2\beta _{1}AE\alpha _{1}=-0.1(\beta _{1}(1-3)A)-2{\frac {10}{3}}\beta _{1}$ (3.9.13) Bring everything to one side and factor out $\displaystyle \beta _{1}$ $\displaystyle \beta _{1}(2AE\alpha _{1}-0.2A+{\frac {20}{3}})=0$ (3.9.14) $\displaystyle \beta _{1}=0$ (3.9.15) $\displaystyle 2AE\alpha _{1}-0.2A+{\frac {20}{3}}=0$ (3.9.16) $\displaystyle \alpha _{1}=\left(0.2A-{\frac {20}{3}}\right)\left({\frac {1}{2AE}}\right)$ (3.9.17) Plug (3.9.17) into (3.9.7) and the candidate solution becomes: $\displaystyle u(x)=.001+\left({\frac {0.1}{E}}-{\frac {10}{3AE}}\right)(x-3)$ (3.9.18) From pb.3.1 from Fish and Belytschko, the equilibrium equation is: $\displaystyle {\frac {d}{dx}}\left(AE{\frac {du}{dx}}\right)+2x=0$ (3.9.19) Substitute (3.9.18) into (3.9.19) $\displaystyle {\frac {d}{dx}}\left(AE\left({\frac {0.1}{E}}-{\frac {10}{3AE}}\right)\right)+2x=0$ (3.9.20) $\displaystyle 0+2x\neq 0$ (3.9.21) The equilibrium equation is not satisfied. From pb.3.1 from Fish and Belytschko, the boundary condition is: $\displaystyle \sigma (1)=\left(E{\frac {du}{dx}}\right)_{x=1}=0.1$ (3.9.22) Substitute (3.9.18) into (3.9.22): $\mathrm{\sigma <!--\; \sigma \; -->}(1)={(E\left(\frac{0.1}{E}-<!--\; -\; -->\frac{10}{3AE}\right)}_{}$