Problem 3.1: Solving General 1-D Model with Weighted Residual Form and Polynomial Basis Functions[edit]
Problem Description[edit]
- Consider
when 
- Choose
, such that 
1. Let n=2
ndof = n+1 = 2+1 = 3 with
.
2. Find 2 equations that enforce the boundary conditions for
.
|
(3.1.1)
|
3. Find 1 more equation (i.e. j = 0,1,2) to solve for
by projecting the residue,
, on a basis function,
, with k = 0,1,2 such that the additional equation is linearly independent from the equations in part 2.
-
- The general form of the residue is defined as
. The exact solution is expressed as
. [1]
approximates
. In other words
. Thus,
(for all values of x in the domain).
- Notes regarding projection:
[2]
[3]
- Combining the previous two equations yields
.
|
(3.1.2)
|
4. Display 3 equations in matrix form,
.
5. Solve for
.
6. Construct
and plot versus the exact solution,
.
7. Repeat steps 1 through 6 for:
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 \displaystyle \Omega =]{0,1}[}](https://wikimedia.org/api/rest_v1/media/math/render/svg/86d648b6143ea55c2f1f57ee00b7f98872f68fa8)


Static case since
.
The natural (eq 3.1.3) and essential (ed 3.1.4) boundary conditions are defined as
|
(3.1.3)
|
|
(3.1.4)
|
Solution[edit]
For the case when 
Parts 1 & 2[edit]
where 
.
Choosing
to avoid
yields
.
The natural boundary condition can be implemented by differentiating
with respect to
and taking its value at x=0.


So the relevant equation is
|
(3.1.5)
|
The essential boundary condition,
, is implemented as follows:


So the relevant equation is
|
(3.1.6)
|
Equation 3.1.2 can be used to project the residue onto the basis function. First, the partial differential equation
is defined. [4]
|
(3.1.7)
|
Substituting the given conditions into this equation and recalling equation 3.1.1,
|
(3.1.8)
|
Differentiating with respect to
and substituting known values,
|
(3.1.9)
|
Substituting into 3.1.2
|
(3.1.10)
|
Since the equation is valid for all
is utilized.
![{\displaystyle \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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0701ec4f3097caf8ec956bb6126288ca208793c7)
- 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
|
(3.1.11)
|
The coefficient (
) and constant (
) matrices are constructed from the equations determined aboved.
- Notes:
- The first row of the matrix
is based on equation 3.1.1 and the essential boundary condition, 3.1.4, when
.
- The second row of the matrix
is from equations 3.1.5 and 3.1.6 (for the general case).
- The third row of the matrix
is from equation 3.1.11.
Note is not Symmetric.
|
(3.1.12)
|
To solve for the matrix
, first recognize that
can be rewritten as
.

|
(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[edit]
The approximate function,
, can now be written using the determined coefficients in equation 3.1.13. The exact solution, was obtained and is displayed here:

The same procedure is repeated for
and
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
plotted against the exact solution in the domain
. It is clear that
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
The previous Matlab code also calculates the error in the estimate for u or the following equation:
|
(3.1.14)
|
This equation was evaluated at .5 and plotted vs n giving the following figure:
Johnathan Whittaker Bullard
Reviewers / Team Members[edit]
Problem 3.2: Fish and Belytschko 2.1[edit]
Problem Description[edit]
Consider figure 2.16 from "A First Course in Finite Elements" by Fish and Belytschko p. 37 pb.2.1[5]

- 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[edit]

The elements are numbered in and the nodes are labled in
|
|
Assemble the global stiffness and force matrix.
,
,
,
,
,
,
,
,
Assembled system matrix:
|
|
The displacement and force matrices for the system are:
The global system of equations is given by:
|
|
Partition the system and solve for the nodal displacements.
As the first two displacements are prescribed, we partition after two rows and columns:

Where
,
,
,
,
,
,

Solving for the nodal displacements:

The systems of equations are:
Muliply by -3
and
Summing the two equations lets us solve for the displacements as:
|
|
Compute the reaction forces:
The global system now becomes:
Solving for the reation forces:
The equations for the reaction forces are:
The reaction forces equal:
and
|
|
Brandonhua 04:48, 15 February 2011 (UTC)
Reviewers / Team Members[edit]
Problem 3.3: Fish and Belytschko 2.3[edit]
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
and cross-sectional area for all bars are
.
- Number the elements and nodes.
- Assemble the global stiffness and force matrix.
- Partition the system and solve for the nodal displacements.
- Compute the stresses and reactions.
Solution[edit]
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
with respect to positive x-axis.
Then,
|
(Eq 3.1)
|
|
(Eq 3.2)
|
|
(Eq 3.3)
|
|
(Eq 3.4)
|
|
(Eq 3.5)
|
Element 2 is numbered with global nodes 2 and 4. It is positioned at an angle
with respect to positive x-axis.
|
(Eq 3.6)
|
|
(Eq 3.7)
|
|
(Eq 3.8)
|
|
(Eq 3.9)
|
|
(Eq 3.10)
|
Element 3 is numbered with global nodes 3 and 4. It is positioned at an angle
with respect to positive x-axis.
|
(Eq 3.11)
|
|
(Eq 3.12)
|
|
(Eq 3.13)
|
|
(Eq 3.14)
|
|
(Eq 3.15)
|
Element 4 is numbered with global nodes 2 and 3. It is positioned at an angle
with respect to positive x-axis.
|
(Eq 3.16)
|
|
(Eq 3.17)
|
|
(Eq 3.18)
|
|
(Eq 3.19)
|
|
(Eq 3.20)
|
Assemble the global stiffness matrix
|
(Eq 3.21)
|
The external force and reaction matrix
|
(Eq 3.22)
|
|
(Eq 3.23)
|
And the displacement matrix
|
(Eq 3.24)
|
3. Global system of equation
|
(Eq 3.25)
|
Partition the system
|
|
Solve for the nodal displacements
|
(Eq 3.26)
|
|
(Eq 3.27)
|
4. The reaction matrix is
|
(Eq 3.28)
|
The stresses in the two elements are
|
(Eq 3.29)
|
For element 1:
|
(Eq 3.30)
|
|
(Eq 3.31)
|
For element 2:
|
(Eq 3.32)
|
|
(Eq 3.33)
|
For element 3:
|
(Eq 3.34)
|
|
(Eq 3.35)
|
For element 4:
|
(Eq 3.36)
|
|
(Eq 3.37)
|
Jiang Jin
Reviewers / Team Members[edit]
Problem 3.4: Show Asymmetry of Weighted Residual Stiffness Matrix[edit]
Problem Description[edit]
- Prove the asymmetry of the weighted residual stiffness matrix by showing:
|
(3.4.1)
|
The stiffness matrix is defined as:
|
(3.4.2)
|
Where < > is for the inner product. Here uh is an approximation of u from a summation of basis functions given as:
|
(3.4.3)
|
Also, P(u) is given by:
|
(3.4.4)
|
Plugging Eq. 3.4.3 into Eq. 3.4.4 and looking only at one term of the summation gives:
|
(3.4.5)
|
Which for a different numbered basis function than bj(x) is written as:
|
(3.4.6)
|
Plugging 3.4.5 and 3.4.6 into 3.4.2 gives:
|
(3.4.7)
|
|
(3.4.8)
|
Which our are equations for Eq.3.4.1, so we now must prove:
|
(3.4.9)
|
Solution[edit]
First, expanding out the integrands of Eq. 3.4.1:
|
(3.4.10)
|
|
(3.4.11)
|
So now we're simply proving:
|
(3.4.12)
|
This will be proven using the following basis functions:
|
(3.4.13)
|
|
(3.4.14)
|
Where
.
Plugging Eq. 3.4.13 and 3.4.14 into Eq. 3.4.10 gives:
|
(3.4.15)
|
Repeating for Eq. 3.4.11:
|
(3.4.16)
|
Expanding out Eq. 3.4.15 and 3.4.16:
|
(3.4.17)
|
|
(3.4.18)
|
Inspecting Eq. 3.4.17 and Eq. 3.4.18 shows:
|
(3.4.19)
|
Braden Snook
Reviewers / Team Members[edit]
Problem 3.5: Fish and Belytschko 2.2[edit]
Problem Statement[edit]
Jacob Fish and Ted Belytschko, " A first Course in Finite Elements",John Wiley & Sons, Ltd, Chapter 2, P37 Problem 2.2.[6]
" 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
|
(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[edit]
The bar can be sub-divided into 3 elements.
The node numbers 1,2,3,4 are assgined.
The element have stiffness .
We know that stiffness 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
|
(3.5.2)
|
Where
|
(3.5.3)
|
For element 2, the stiffness matrix is
|
(3.5.4)
|
Where
|
(3.5.5)
|
For element 3, the stiffness matrix is
|
(3.5.6)
|
Where
|
(3.5.7)
|
Now assembling the stiffness matrices we get,
|
(3.5.8)
|
|
(3.5.9)
|
We Know that
Where d is displacement and F is force. Taking node(1) to be fixed end.
|
(3.5.10)
|
|
(3.5.12)
|
Solving the above we get,
|
(3.5.13)
|
|
|
(3.5.13)
|
Srilalithkumar 15:41, 15 February 2011 (UTC)
Reviewers / Team Members[edit]
Problem 3.6: Fish and Belytschko 2.4[edit]
Problem Description[edit]
"A First Course in Finite Elements" by Fish and Belytschko p. 38 Problem 2.4[7]
Given the three-bar structure subjected to the prescribed load ad point C equal to as shown. The Young’s modulus is Pa. The cross-sectional area of the bar BC is and that of BD and BF is
. 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[edit]
a)
First subdivide the structure into elements. There are three elements and two nodes. Shown in Fig.
Element 1:
|
|
|
(3.6.1)
|
Element 2:
|
|
|
(3.6.2)
|
Element 3:
|
|
|
(3.6.3)
|
The global stiffness matrix is:
|
(3.6.4)
|
And
|
(3.6.5)
|
|
(3.6.6)
|
|
(3.6.7)
|
b)
Global system of equations:
|
(3.6.8)
|
Then reduce the global system of equations:
The global system is partitioned four rows and four columns:
And:
|
(3.6.9)
|
So:
|
(3.6.10)
|
|
(3.6.11)
|
To sum up, the displacement in X-direction of B is mm, in Y-direction is mm. The displacement in X-direction of C is mm.
c)
Element 1:
|
(3.6.12)
|
|
(3.6.13)
|
Element 2:
|
(3.6.14)
|
|
(3.6.15)
|
Element 3:
|
(3.6.16)
|
|
(3.6.17)
|
d)
Considering the Global system of equations, then:
|
(3.6.18)
|
So
|
(3.6.19)
|
|
(3.6.20)
|
|
(3.6.21)
|
|
(3.6.22)
|
Zongyi Yang
Reviewers / Team Members[edit]
Problem 3.7: Solving General 1-D Model with Weighted Residual Form and Polynomial Basis Functions[edit]
Problem Description[edit]
Repeat problem 3.1 with
Solution[edit]
For the case when 
Parts 1 & 2[edit]
where 
.
The natural boundary condition can be implemented by differentiating with respect to and taking its value at x=0.


So the relevant equation is
|
|
The essential boundary condition, , is implemented as follows:


So the relevant equation is
|
|
Equation 3.1.2 can be used to project the residue onto the basis function. First, the partial differential equation is defined. [8]
|
|
Substituting the given conditions into this equation and recalling equation 3.1.1,
|
|
Differentiating with respect to and substituting known values,
|
|
Substituting into 3.1.2
|
|
Since the equation is valid for all is utilized.
![{\displaystyle \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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0701ec4f3097caf8ec956bb6126288ca208793c7)
- 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
|
|
The coefficient ( ) and constant ( ) matrices are constructed from the equations determined aboved.
- Notes:
- The first row of the matrix
is based on equation 3.1.1 and the essential boundary condition, 3.1.4, when .
- The second row of the matrix
is from equations 3.1.5 and 3.1.6 (for the general case).
- The third row of the matrix
is from equation 3.1.11.
Note is not symmetric
|
|
To solve for the matrix , first recognize that can be rewritten as .

|
|
|
Parts 6 & 7[edit]
The approximate function, , can now be written using the determined coefficients. The exact solution, was previously obtained and is displayed here:

The same procedure is repeated for and 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 plotted against the exact solution in the domain . It is clear that provides a good approximation. This is likely due to the fact that the approximate and exact solutions are both second order polynomial functions.
The previous Matlab code also calculates the error in the estimate for u or the following equation:
|
|
This equation was evaluated at .5 and plotted vs n giving the following figure:
Johnathan Whittaker Bullard
Reviewers / Team Members[edit]
Whitaker Bullard/Braden Snook
Problem 3.8: Fish and Belytschko 3.1[edit]
Problem Description[edit]
- [From Fish & Belytschko, Chapter 3, Problem 3.1] [9]
Show that the weak form of



is given by

Solution[edit]
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] [10] First, the governing equation and traction boundary condition are multiplied by an arbitrary (or weight) function, , and integrated from 1 to 3.
|
(3.8.1)
|
|
(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
|
(3.8.3)
|
Using integration by parts we obtain
|
(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 , respectively. is defined as the boundary where the displacements are prescribed (x=3 in this case); 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
|
(3.8.5)
|
The original traction boundary can be applied to the second term. Recalling that removes the first term. Thus we are left with the weak form in the case of one-dimensional stress analysis.
|
(3.8.6)
|
Philip Flater
Reviewers / Team Members[edit]
Problem 3.9: Fish and Belytschko 3.3[edit]
Problem Description[edit]
"A First Course in Finite Elements" by Fish and Belytschko p. 72, pb.3.3[11]
Consider a trial (candidate) solution of the form 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[edit]
The candidate solution:

|
(3.9.1)
|
The weight solution:

|
(3.9.2)
|
From pb.3.1 from Fish and Belytschko (eqn 3.8.6 from this hw submission):

|
(3.9.3)
|
The weight function 


|
(3.9.4)
|
From pb.3.1 from Fish and Belytschko:

|
(3.9.5)
|

|
(3.9.6)
|

|
(3.9.6)
|
The candidate solution and the weight function is now:

|
(3.9.7)
|

|
(3.9.8)
|
Recall that the weak form is

|
(3.9.9)
|
From (3.9.7) and (3.9.8):

|
(3.9.10)
|

|
(3.9.11)
|
Plug (3.9.7), (3.9.8), (3.9.10) and (3.9.11) into (3.9.9):

|
(3.9.12)
|
Integrage (3.9.12)

|
(3.9.13)
|
Bring everything to one side and factor out 

|
(3.9.14)
|

|
(3.9.15)
|

|
(3.9.16)
|

|
(3.9.17)
|
Plug (3.9.17) into (3.9.7) and the candidate solution becomes:

|
(3.9.18)
|
|
|
From pb.3.1 from Fish and Belytschko, the equilibrium equation is:

|
(3.9.19)
|
Substitute (3.9.18) into (3.9.19)

|
(3.9.20)
|

|
(3.9.21)
|
The equilibrium equation is not satisfied.
|
|
From pb.3.1 from Fish and Belytschko, the boundary condition is:

|
(3.9.22)
|
Substitute (3.9.18) into (3.9.22):
|
|