User:Eml5526.s11.team4/HW4

= Problem 4.1 - Residue Projection for Modified Trigonometric Function Basis =

This problem refers to the previous Problem 2.9 and Problem 3.11

Problem Statement
Given the differential equation in strong form,

Consider the basis

Let, $$\displaystyle n=2,4$$ or $$\displaystyle 6$$ where $$\displaystyle ndof=n+1=3, 5$$ or $$\displaystyle 7$$ where ndof is the number of degrees of freedom and $$\displaystyle d=\left\{ {{d}_{j}};j=0,1,2.....n \right\}$$.


 * 1. Find two equations that enforce boundary conditions for


 * 2. Find $$\displaystyle 1,3$$ or $$\displaystyle 5$$ more equation(s) to solve for $$\displaystyle \mathbf{d}={{({{d}_{j}})}_{n\times 1}}

(j=0,1,2....n)$$by projecting the residue $$\displaystyle P({{u}^{h}}$$)on the basis function$$\displaystyle {{b}_{k}}(x)$$with $$\displaystyle k=1,2,3\cdots \cdots

$$ such that the additional equation is linearly independent from the above equation in part2.


 * 3. Display $$\displaystyle 3/5/7$$ equations in matrix form $$\displaystyle \mathbf{KD=F}$$ Observe the symmetry properities of$$\displaystyle \mathbf{K}$$


 * 4. Solve for $$\displaystyle \mathbf{d}$$


 * 5. Construct $$\displaystyle {{u}^{h}}(x)$$and plot $$\displaystyle {{u}^{h}}(x)$$ versus $$\displaystyle {{u}_{exact}}(x)$$

Exact Solution
With the general solution

After applying boundary conditions

For the case n=2
since $$\displaystyle {{b}_{j}}=\left\{ 1,\sin x,\sin 2x \right\}$$ with $$\displaystyle {{d}_{j}}=\left\{ {{d}_{0}},{{d}_{1}},{{d}_{2}} \right\}$$,

Boundary Conditions:

Additional equation:

For the third condition, we project the residue on the basis function.

Where:

hence,

and

thus the projection is

Putting equation 4.1.6, 4.1.7 and 4.1.11 into matrix form, we can have

Result:



For the case n=4
The approximating solution

Enforcing the boundary conditions:

There still three additional equations which been derived from enforcing the projection of $$\displaystyle P(u)$$ on basis to be zero:

putting equations 4.1.15 to 4.1.20 in matrix form yields:



For the case n=6
The approximating solution

Enforcing the boundary conditions:

There still three additional equations which been derived from enforcing the projection of $$\displaystyle P(u)$$on basis to be zero:

putting equations 4.1.25 to 4.1.31 in matrix form yields:

Thus



Comparison of error
The maximum error for each case is listed below

It is seen that for the case n=4 the error is numerically zero, however the error has risen for n=6 making the solution suspect.

Problem Statement
Obtain the weak form for the equations of heat conduction with the boundary conditions $$T(0)=100$$ and $$\bar{q}=hT$$ .The condition on the right is a convection condition.

Solution


w= weighted function

integrating by part yields

The weak form is

Find $$T(x)\epsilon \upsilon $$ such that

=Problem 4.3: Heat conduction problem in a circular plate=

Problem Statement
Given the strong form:

where


 * $$R$$ || is the total radius of the plate
 * $$s$$ || is the heat source per unit length along the plate radius
 * $$T$$ || is the temperature
 * $$k$$ || is the conductivity
 * }
 * $$k$$ || is the conductivity
 * }
 * }

Assume that $$k, s,$$ and $$R$$ are given.

Part A: Construct the weak form
Multiplying the governing equation (1) by a weight function $$ w(r) $$ and integrating over the domains on which they hold: We rewrite Equ (4) in the equivalent form: Taking the derivative of a product for the first term, Since $$w(R)=0$$ and because of (2), the first term in the above vanishes. So the solution is obtained as follows according to Equ (5) and (6):

Part B: Consider quadratic trial solutions and weight functions
From the given information the equations for $$T(r)$$ and $$w(r)$$ are as follows: The weight function $$w(r)$$ must vanish at $$r=R$$. Also, the trial solution $$T(r)$$ must satisfy the essential boundary condition (3). Differentiating $$w(r)$$ and $$T(r)$$ to obtain the solution of the weak form,

Part C: Obtain a solution
Show that the temperature distribution along the radius is given by

Substituting Equ (10,12) into the weak form (7) to solove the differential equation with the boundary conditions, Evaluating the integrals and factoring out $$\beta_1$$ and $$\beta_2$$ gives As the above must hold for all $$\beta_1$$ and $$\beta_2$$, it follows that the terms in the parentheses must vanish. Rearranging it,

We can get $$\alpha_1$$ and $$\alpha_2$$ from the above 2 equations. Substituting (17) into (11) gives the weak solution.

Problem Statement:
Given the strong form for the circular bar in torsion shown above: $$ \begin{align} \frac{\mathrm{d}}{\mathrm{d}x}(JG\frac{\mathrm{d\phi}}{\mathrm{d}x})+m=0 \end{align} $$ $$ \begin{align} 0\leq x\leq 1 \end{align} $$ The natural boundary condition is given by: $$ \begin{align} M(x=l)=(JG\frac{\mathrm{d\phi}}{\mathrm{d}x})_l=\bar{M} \end{align} $$ The essential boundary condition is given by: $$ \begin{align} \phi(x=0)=\bar{\phi} \end{align} $$ Where m(x) is a distributed moment per unit length, M is the torsion moment, $$\phi$$ is the angle of rotation, G is the shear modulus, and J is the polar moment of inertia given by $$J=\frac{\pi C^{4}}{2}$$. In this equation, C is the radius of the circular shaft. 1.) Construct the weak form for the circular bar in torsion. 2.) Assume that m(x)=0 and integrate the differential equation given above. Find the integration constants using boundary conditions.

Solution:
1.) To construct the weak form, first, the strong form is multiplied by the weighting function $$w(x)$$ producing equation 4.4.1. Next, equation 4.4.1 is integrated over the domain $$\Omega$$ as shown in equation 4.4.2. This integral can be solved using integration by parts as derived in Fish and Belytschko Section 3.5.2 (shown below in general form). For the more specific case: $$ \Gamma_u=\left\{x=0\right\}$$ $$ \Gamma_t=\left\{x=L\right\}$$ $$ n\left\{x=0\right\}=-1$$ $$ n\left\{x=L\right\}=1$$ $$ f=JG\frac{\partial\phi}{\partial x} $$ Combining this with equations 4.4.2 and 4.4.3 produces: Given that $$w(x=0)=0$$ and by applying the natural boundary condition produces the weak form shown in equation 4.4.5. 2.) Assuming $$m(x)=0$$ produces the new strong form equation 4.4.6. Integrating with respect to $$x$$ results in equation 4.4.7. By applying the previously stated natural boundary condition, the coefficient of integration $$C_{1}$$ can be obtained as shown below.

References:
Solved by Dustan Simko

=Problem 4.5: Constraint Breaking for Weak Form Basis Functions=

Problem Statement:
Consider the following potential basis functions: Let $$\Omega=\left[\alpha,\beta\right]=\left[-2,4\right]$$. For $$\Gamma_{g}=\left\{\beta\right\}$$ and each of the above listed families of basis functions, find the corresponding family $$\tilde{F_{i}}=\left\{\mathbf{\bar{b}_{j}}\right\}$$ satisfying the constraint breaking solution (CBS) requirements. Recall that the CBS is a family of basis functions $$\left\{\bar{b}_j(x),j=0,1,...,n\right \}$$ such that $$\bar{b}_{0}(\beta)\neq 0$$ and $$ \bar{b}_{j}(\beta)=0$$ for $$j=1,...,n$$ 1.) Plot $$\mathbf{\left\{b_{j}\right\}}$$ and $$\left\{\mathbf{\bar{b}_{j}}\right\}$$ for $$j=1,2,3$$ 2.)Show $$\mathfrak{F}_{e}=\left\{e^{jx},j=0,1,2,3\right\}$$ is linearly independent on $$\Gamma$$
 * $$\mathbf{Polynomial}$$
 * $$\mathbf{Cosine}$$
 * $$\mathbf{Sine}$$
 * $$\mathbf{Fourier}$$
 * $$\mathbf{Exponential}$$

Polynomial
The polynomial basis function is given by $$\mathfrak{F}_{p}=\left\{x^{j},j=0,1,...,n\right\}$$. For $$\left\{j=0,1,2,3\right\}$$ the basis functions $$\left\{b_{j}\right\}$$ are as follows: The constraint breaking solution$$\left\{\bar{b}_j\right\}$$ for the polynomial basis function is shown below: In the specific case where $$\left\{\beta=4\right\}$$ the CBS becomes: The basis functions$$\left\{\bar{b}_j\right\}$$ satisfy the CBS requirements since$$\left\{\bar{b}_0(\beta)=constant\neq0\right\}$$ and $$\left\{\bar{b}_j(\beta)=0, for j=1,2,...,n\right\}$$. The following figure shows the original polynomial basis functions and the CBS basis functions for $$\left\{j=1,2,3\right\}$$:

Cosine
The cosine basis function is given by $$\mathfrak{F}_{c}=\left\{cos(jx),j=0,1,...,n\right\}$$. For $$\left\{j=0,1,2,3\right\}$$ the basis functions $$\left\{b_{j}\right\}$$ are as follows: The constraint breaking solution$$\left\{\bar{b}_j\right\}$$ for the cosine basis function is shown below: In the specific case where $$\left\{\beta=4\right\}$$ the CBS becomes: The basis functions$$\left\{\bar{b}_j\right\}$$ satisfy the CBS requirements since$$\left\{\bar{b}_0(\beta)=constant\neq0\right\}$$ and $$\left\{\bar{b}_j(\beta)=0, for j=1,2,...,n\right\}$$. The following figure shows the original cosine basis functions and the CBS basis functions for $$\left\{j=1,2,3\right\}$$:

Sine
The sine basis function is given by $$\mathfrak{F}_{c}=\left\{1,sin(jx),j=1,...,n\right\}$$. For $$\left\{j=0,1,2,3\right\}$$ the basis functions $$\left\{b_{j}\right\}$$ are as follows: The constraint breaking solution$$\left\{\bar{b}_j\right\}$$ for the sine basis function is shown below: In the specific case where $$\left\{\beta=4\right\}$$ the CBS becomes: The basis functions$$\left\{\bar{b}_j\right\}$$ satisfy the CBS requirements since$$\left\{\bar{b}_0(\beta)=constant\neq0\right\}$$ and $$\left\{\bar{b}_j(\beta)=0, for j=1,2,...,n\right\}$$. The following figure shows the original sine basis functions and the CBS basis functions for $$\left\{j=1,2,3\right\}$$:

Fourier
The Fourier basis function is given by So, the basis functions $$\left\{b_{j,k}\right\}$$ are as follows: The corresponding basis functions $$\left\{\bar{b}_{j,k}\right\}$$ satisfying CBS are shown below: The basis functions$$\left\{\bar{b}_{j,k}\right\}$$ satisfy the CBS requirements


 * $$\left\{\bar{b}_0(\beta)=1\neq0\right\}$$


 * $$\left\{\bar{b}_{j,k}(\beta)=0, for j,k=1,2,...\right\}$$.

The following figure shows the original fourier basis functions and the CBS basis functions for $$\left\{j=1,2,3;k=1,2\right\}$$:

Exponential
The exponential basis function is given by So, the basis functions $$\left\{b_j\right\}$$ are as follows: The corresponding basis functions $$\left\{\bar{b}_j\right\}$$ satisfying CBS are shown below: The basis functions$$\left\{\bar{b}_j\right\}$$ satisfy the CBS requirements


 * $$\left\{\bar{b}_0(\beta)=1\neq0\right\}$$


 * $$\left\{\bar{b}_{j}(\beta)=0, for j=1,2,...\right\}$$.

The following figure shows the original exponential basis functions and the CBS basis functions for $$\left\{j=1,2,3\right\}$$:

References:
Solved by Dustan Simko and Woo Rhee

= Problem 4.6 - Strong and weak form of elastic bar problem with variable distributed spring =

Problem Statement
Consider an elastic bar with a variable distributed spring p(x) along its length .The distributed spring imposes an axial force on the bar in proportion to the displacement. Consider a bar of length l,cross-sectional area A(x),young's modulus E(x)with body force b(x) and boundary conditions as follows

a. Construct the strong form

b. Construct the weak form

Solution
First consider an element of this elastic bar.



From the force equilibrium i.e D'Alembert's principle,

taking the limit as $$\Delta x\rightarrow 0$$

for the governing ODE

Boundary conditions are

So the strong form is

The weak form derivation:

integrating by parts for the first term yields

which can be expanded as follows

Let w (0)=0 be the essential boundary condition.

The the weak form is as follows

Find u(x) smooth enough and satisfy u(0)=0,such that

=Problem 4.7: CGX Project and Tutorial=

From the lecture slide Mtg23

Problem Statement

 * 1) install cgx.
 * 2) Read manual; sign up with user group to ask questions if any, also access archive.
 * 3) Reproduce the basic examples:disk, cylinder, sphere, sphere-volume, airfoil.
 * 4) Write a report for "dummies": explain how to install and run cgx.

Final Report
Installation and run:
 * 1) Download the windows-version of Calculix from its website, and then install it on computer.
 * 2) Run the CalculiX Command: [Start] [Programs] [bConverged][CalculiX] [CalculiX Command]
 * 3) Reproduce the basic examples

A. The “disc” problem
1. First switch to the directory of working ‘E:\Documents\ccx’.

2. Then input the command ‘cgx –b disc.fbd’, there will appear another command “CalculiX Graphix”.

3. Input the command” pnt [  ]” and then plot all the points with the command ‘plot p all b’. For example, ”PNT py  0.0  1.0 0.0” is used to define a point which is located at (0,1,0).

PNT py     -0.00000        1.00000        0.00000

The figure is showed below.



4. Input the command” line    ” and the “gsur '+|-'BLEND| ' '+|-'  '+|-' -> .. (3-5 times)and then plot all the lines and surfaces with command ‘plus’. For example, “LINE L001 P00I P001 p0 4 ” is used to define a arc line which pass through points P00I and P001 and its center point is p0 and radius is 4.

LINE L001 P00I P001 p0 4

GSUR A001 + BLEND - L003 - L002 - L001 - L004

The figure is showed below.



5. You can use mouse on the Graphix command to rotate and zoom in or out the figure.

6. Input the command ‘ELTY all QU4’ to assign the element type and then ‘mesh all’ to mesh the geometry, and then left-click on the marginal area on the Graphix, then on ’Viewing’ and continue to ‘show all elements with light’ and ‘Toggle Element edges’, as showed in following figures





B. The “cylinder” problem
1. First switch to the directory of working ‘E:\Documents\ccx’.

2. Then input the command ‘cgx –b cylinder.fbd’, there will appear another command “CalculiX Graphix”.

3. Input the command” pnt [  ]” and then plot all the points. Input the command” line    ” and the “gsur '+|-'BLEND| ' '+|-'  '+|-' -> .. (3-5 times)and then plot all the lines and surfaces. The figure is showed below



4. Use the command “seta |'n'|'e'|'p'|'l'|'c'|'s'|'b'|'S'|'L'|'se' <name ..> ['n'|'e' <name” to define some sets to be used in further loading operations. For example, “SETA p1 p p1” is to define a set of points which contains point p1.

5. You can use mouse on the Graphix command to rotate and zoom in or out the figure.

6. Input the command ‘ELTY all QU4’ and ‘mesh all’, and then left-click on the marginal area on the Graphix, then on ’Viewing’ and continue to ‘show all elements with light’ and ‘Toggle Element edges’, as showed in following figure.



C. The “sphere” problem
1. First switch to the directory of working ‘E:\Documents\ccx’

2. Then input the command ‘cgx –b sphere.fbd’, there will appear another command “CalculiX Graphix”.

3. Input the command” pnt [<x> <y> <z>]” and then plot all the points with the command ‘plot p all b’. For example, ”PNT py  0.0  1.0 0.0” is used to define a point which is located at (0,1,0). The figure is showed below.



4. Input the command” line <p1> <p2> <cp|seq> ” and the “gsur '+|-'BLEND| ' '+|-' <line|lcmb> '+|-' -><line|lcmb> .. (3-5 times)and then plot all the lines and surfaces with command ‘plus’. For example, “LINE L001 P00I P001 p0 104 ” is used to define a arc line which pass through points P00I and P001 and its center point is p0 and radius is 104. The figure is showed below.



5. Input the command ‘ELTY all QU4’ and ‘mesh all’, and then left-click on the marginal area on the Graphix, then on ’Viewing’ and continue to ‘show all elements with light’ and ‘Toggle Element edges’, as showed in following figure.



D. The “sphere-volume” problem
1. First switch to the directory of working ‘E:\Documents\ccx’

2. Then input the command ‘cgx –b sphere-volu.fbd’, there will appear another command “CalculiX Graphix”.

3. Input the command” pnt [<x> <y> <z>]” and then plot all the points. Input the command” line <p1> <p2> <cp|seq> ” and the “gsur '+|-'BLEND| ' '+|-' <line|lcmb> '+|-' -><line|lcmb> .. (3-5 times)and then plot all the lines and surfaces. The figure is showed below



4. Input the command” gbod 'NORM' '+|-' '+|-' ->( 5-7 surfaces )” to define a volume. For example, “GBOD B001 NORM + A006 - A003 - A004 + A002 + A001” is to define a volume which consists of these above areas. GBOD B001 NORM + A006 - A003 - A004 + A002 + A001

5. Input the command ‘ELTY all he20’ and ‘mesh all’, and then left-click on the marginal area on the Graphix, then on ’Viewing’ and continue to ‘show all elements with light’ and ‘Toggle Element edges’, as showed in following figure.



E. The “airfoil” problem
1. First switch to the directory of working ‘E:\Documents\ccx’

2. Then input the command ‘cgx –b airfoil.fbd’, there will appear another command “CalculiX Graphix”.

3. Input the command” pnt [<x> <y> <z>]” and then plot all the points with the command ‘plot p all b’. For example, ” PNT P05D     0.00118       -5.37045        0.00000” is used to define a point which is located at (0.00118,-5.37045,0.00000). The figure is showed below.



4. Input the command ”seqa ['pnt' .. <=>]|['afte'|'befo' .. <=>]|['end' .. <=>]” to define a sequential set. For example, “SEQA S005 pnt  P010 P00C P00A P01C” is used to define a sequential point set which pass through points P010 P00C P00A P01C sequentially.

SEQA S006    pnt  P01A P00A P00C P00E P00G P00I P00J P00K P00L P019

5. Input the command ”lcmb ['+|-' '+|-' '+|-' ->  ..(up to 14 lines)]| ['ADD' '+|-' '+|-' -> '+|-' ..(up to 14 lines)]” to combine lines. For example, “ LCMB C001 + L015 + L01X + L01Y + L016” is used to combine the lines of C001, L015 + L01X + L01Y, L016.

LCMB C001 + L001 + L002 + L00A + L00L + L009

6. Input the command” line <p1> <p2> <cp|seq>” and the “gsur '+|-'BLEND| ' '+|-' <line|lcmb> '+|-' -><line|lcmb> .. (3-5 times)and then plot all the lines and surfaces with command ‘plus’. For example, “LINE L001 P00I P001 p0 104 ” is used to define a arc line which pass through points P00I and P001 and its center point is p0 and radius is 104. The figure is showed below.



7. Input the command ‘ELTY all QU4’ and ‘mesh all’, and then left-click on the marginal area on the Graphix, then on ’Viewing’ and continue to ‘show all elements with light’ and ‘Toggle Element edges’, as showed in following figure.





References:
=Problem 4.8 Finding the mass and stiffness matrix for a dynamic system =

The following problem was assigned in Lecture 23 and is a modification to Fish and Belytschko's problem 3.4 on page 72.

Given data
The following given partial differential equation (PDE) is provided in strong form with boundary and initial conditions. The PDE is a second-order constant coefficient wave equation with a first order term.

The basis and weighting functions are provided respectively as

With

Problem statement
Find $$\textbf{F}(t)$$ and the mass matrix, $$\displaystyle{ \tilde{ \mathbf{ M}}}$$

Developing the Weak Form
Starting with the strong form we derive the weak form by multiplying an arbitrary weight function and then integrate by parts.

Evaluating with the basis and weighting functions
Evaluating term by term.

Here we note that $$u^h$$ is a function of time as well as space. For this to be true, the coefficients $$a_0,a_1,a_2, a_3$$ must be functions of time. Note that in a final Finite Element formulation, such as a Galerkin Formulation, $$u^h$$ is assumed to be separable into spatially dependent and time dependent terms. This would allow us to explicitly integrate in time (the exception being for nonlinear equations).

Further, note

This reduces the unknowns to three, and provides three governing equations. Evaluating the temporal term gives the last components.

Solving the system of equations
By the generality in $$\beta_i$$ we now have a system of equations, namely

We now have three equations and three unknowns $$(\alpha_1(t), \alpha_2(t), \alpha_3(t))$$. We note that $$\alpha_0''=-4\sin(2t)$$

Linearizing the system of equations
We now linearize the system by introducing three new equations

We now have six equations and six unknowns. The task now remains to solve six simultaneous ODE's.

So that we now have the linear system of ODE's, which may be solved as follows

Calculate an integrating factor

Now multiplying by the integrating factor and recognize the left hand side is the total derivative

In theory the solution can be calculated with a symbolic toolbox. However, the solve times were excessively large are intractable. This suggests another approach or a complicated temporal behavior. We note that numerical integration would overcome this challenge.

The mass and stiffness matrices
As requested we list the mass and stiffness matrices from equation (8.12)