**Programming project in TMA4220 – part 2**

**Programming project in TMA4220 – part 2**

编程作业代写 In the second part of the project you are going to make use of your finite element library which you have now built.

**Do real-life experimentation using your FEM code**

In the second part of the project you are going to make use of your finite element library which you have now built. We are going to apply this to one of several real-life applications. The three first tasks introduce the main equations you can choose from

- – the linear elasticityequation
*ρ**∂t*2 =*∇*(*σ*) – the free vibrationequation*u*

and include some very rough ”getting-started” theory on these.

You should choose **one **of the following exercises (i.e. either 1, 2, 3 og 4).

**Exercise 1: The heat**** ****equation 编程作业代写**

**Exercise 1: The heat**

**equation 编程作业代写**

** **

** **Figure 1: A candle

In this exercise we will model heat transfer in a thin metal sheet when a candle is placed under it.

The heat equation reads

(1)

where *α *is a positive constant defined by

* *

with *κ** *being the thermal conductivity, *ρ** *the mass density, *c*_{p}* *the specific heat capacity of the material and *f** *a heat source term.

We will let the computational domain Ω be a two dimensional metal sheet, and we want to find the temperature distribution in the metal sheet when a candle fire is placed under it.

**a) ****Semidiscretization**

**Semidiscretization**

We are going to *semidiscretize *the system by projecting the spatial variables to a finite element subspace *X** _{h}*. First, assume the source term

*f*is 0, and multiply (1) by a test function

*v*and integrate over the domain Ω to get

Note that we have only semidiscretized the system, so our unknown *u *is given as a linear combination of the *spatial *basis functions, and is continuous in time, i.e.

The variational form of the problem then reads: Find *u*_{h}* **∈ **X*^{D}* *such that

which in turn can be written as the linear system

(2)

which is an ordinary differential equation (ODE) with the matrices defined as

Construct the matrix ** A **and

**as defined above.**

*M***b) ****Include the source**** ****term**

**Include the source**

**term**

Assume the heat from the candle light is given by a function

for some parameter *β > *0.

What changes in the variational formulation from (a) do you need in order to include this?

**c) ****Time**** ****integration**

**Time**

**integration**

The system (2) is an ODE, which should be familiar from previous courses. Very briefly an ODE is an equation on the form

where *y *may be a vector. The simplest ODE solver available is Eulers method

*y*_{n}_{+1} = *y*_{n}* *+ *hg*(*t*_{n}*, y** _{n}*)

*.*

More sophisticated include the improved Eulers methods

*y*_{n}_{+1} = *y*_{n}* *+ (*g*(*t*_{n}*, y** _{n}*) +

*g*(

*t*

_{n}_{+1}

*, y*

_{n}*+*

*hg*(

*t*

_{n}*, y*

*)))*

_{n}or the implicit trapezoid rule

*y*_{n}_{+1} = *y*_{n}* *+ (*g*(*t*_{n}*, y** _{n}*) +

*g*(

*t*

_{n}_{+1}

*, y*

_{n}_{+1}))

and the famous Runge-Kutta methods.

Choose an ODE scheme (based on your previous experience and expertise) and implement your time integration. Why did you choose the solver you did?

**d) ****Experimentation**

**Experimentation**

The parameters *κ*, *ρ *and *c*_{p}* *are problem specific. Try to find (by for example googling it) typical values for these for different types of metal.

The free parameters that you can change in this problem are the boundary conditions

*u*(*t, x, y*)*|*_{∂}_{Ω} = *u** ^{D}*(

*t, x, y*)

and the initial condition

*u*(*t,** **x,** **y*)*|*_{t}_{=0} = *u*_{0}(*x,** **y*)*.*

What do you think are appropriate boundary and initial conditions? Why?

The parameter *β *in the source term decides how large the candle light is. What happens if you let *β *be large (for example *β *= 1000)? 编程作业代写

Use a source function

*f** *(*t, x, **y*) = *e**−**β*((*x**−**a *sin(*t*))^{2}+*y*^{2})

for different values of *a*. What happens to the solution when *a *is large? Do you have to change the time integration somehow?

**e)**

**e)**

When the analytical solution is not known, a way to compute the convergence rate is to compare the numerical solution with the numerical solution computed using very small elements.

Compute the solution for element sizes *h **∈ {*__ ____1__ *, i *= 0*, *1*, *2*, *3*} *and compare these with the solution using element size *h *= . Compute the relative error and the convergence rate.

**Exercise 2: Linear elasticity equation ****编程作业代写**

**Exercise 2: Linear elasticity equation**

**编程作业代写**

We are in this problem going to consider the linear elasticity equation. The equations describe deformation and motion in a continuum. While the entire theory of continuum mechanics is an entire course by itself, it will here be sufficient to only study a small part of this: the linear elasticity. This is governed by three main variables *u**, *** ε **and

**(see Table 1).**

*σ*Table 1: Linear elasticity variables in two dimensions

Note that the subscript denotes vector component and *not *derivative, i.e. *u*_{x}* **ƒ*= ^{∂u}* *.These three variables can be expressed in terms of each other in the following way:

u=u(x) (3)

ε=ε(u)(4)

σ=σ(ε)(5)

The primary unknown ** u **(the displacement) is the one we are going to find in our finite element implementation. From (3) we will have two displacement values for each finite element ”node”, one in each of the spatial directions.

#### The relation (4) is a purely geometric one. Consider an infinitesimal small square of size *dx *and *dy*, and its deformed geometry as depicted in figure 2.

The strain is defined as the stretching of the element, i.e. *ε*_{xx}* *= ^{length}^{(}^{ab}^{)}^{−}^{length}^{(}^{AB}__ ^{)}__ . The complete derivations of these quantities is described well in the Wikipedia article on strain, and the result is the following relations

(6)

Note that these relations are the *linearized *quantities, which will only be true for small deformations.

For the final relation, which connects the deformation to the forces acting upon it, we turn to the material properties. Again, there is a rich literature on the subject, and different

Figure 2: An infinitesimal small deformed rectangle

relations or physical laws to describe different materials. In our case, we will study small deformations on solid materials like metal, wood or concrete. It is observed that such materials behave elastically when under stress of a certain limit, i.e. a deformed geometry will return to its initial state if all external forces are removed. Experiment has shown that the Generalized Hooks Law is proving remarkable accurate under such conditions. It states the following.

Consider a body being dragged to each side by some stress *σ*_{xx}* *as depicted in Figure 3. Hooks law states that the forces on a spring is linearly dependant on the amount of stretching multiplied by some stiffness constant, i.e. *σ*_{xx}* *= *Eε** _{xx}*. The constant

*E*is called Young’s modulus. Generalizing upon this law, we see that materials typically contract in the

*y*-direction, while being dragged in the

*x*-direction. The ratio of compression vs expansion is called Poisson’s ratio

*ν*

*and is expressed as*

*ε*

_{yy}*=*

*−*

*νε*

*. This gives the following relations*

_{xx}Due to symmetry conditions, we see that when applying a stress *σ*_{yy}* *in addition to *σ*_{xx}

we get

Finally, it can be shown (but we will not) that the relation between the shear strain and shear stress is . Collecting the components of ** ε **and

**in a vector, gives us**

*σ*Figure 3: Deformed geometry under axial stresses

#### the compact notation

or conversely

For a body at static equilibrium, we have the governing equations

and some appropriate boundary conditions

**a) ****Weak**** ****form**

**Weak**

**form**

Show that (8) can be written as the scalar equation

(where we have exchanged the subscripts (*x, **y*) with (1*, *2)) by multiplying with a test function and integrating over the domain Ω Moreover, show that this can be written in compact vector form as

**b) ****Galerkin**** ****projection**

**Galerkin**

**projection**

As in 2b) let ** v **be a test function in the space

*X*

_{h}*of piecewise linear functions on some triangulation*

*T*. Note that unlike before, we now have

*vector*test functions. This means that for each node ˆ

*i*, we will have two test functions

Let these functions be numbered by a single running index *i** *= 2ˆ*i** *+ *d*, where *i** *is the node number in the triangulation and *d *is the vector component of the function.

Show that by inserting * *into (11) you get the system of linear equations

*A*** u **=

*b*where

(**Hint: ***ε*¯(*·*) is a linear operator)

**c) ****Test**** ****case**

**Test**

**case**

Show that

is a solution to the problem

where

and Ω = *{*(*x,** **y*) : max(*|**x**|**,** **|**y**|*) *≤** *1*}** *is the refereance square (*−*1*,** *1)^{2}.

**d) ****Implementation**

**Implementation**

Modify your Poisson solver to solve the problem (11). Verify that you are getting the correct result by comparing with the exact solution. The mesh may be obtained through the Grid function getPlate( ).

**e) ****Stress**** ****analysis**

**Stress**

**analysis**

Solving (8) with a finite element method gives you the primary unknown: the displacement

*u*. If you are interested in derived quantities such as the stresses, these can be calculated from (7). Note that*σ*is in essence the derivative of*u*which means that since*u*is*C*^{0}across element boundaries, then*σ*will be discontinuous. To get stresses at the nodal values, we propose to average the stresses over all neighbouring elements.

Loop over all elements and evaluate (the constant) stresses on that element. For each node, assign the stresses to be the average stress over all neighbouring elements. This method is called ”Stress Recovery”.

**f )**

**f )**

When the analytical solution is not known, a way to compute the convergence rate is to compare the numerical solution with the numerical solution computed using very small elements. 编程作业代写

Compute the solution for element sizes *h **∈ {*__ __*, i *= 0*, *1*, *2*, *3*} *and compare these with the solution using element size *h *= . Compute the relative error and the convergence rate.

**Exercise 3: Vibration analysis**

**Exercise 3: Vibration analysis**

Do problem **2a) – 2d) **and read the theory on linear elasticity.

Figure 4: Mass-spring-model

The forces acting on a point mass *m *by a spring is given by the well known Hooks law:

*m**x*¨ = *−**k**x*

#### This can be extended to multiple springs and multiple bodies as in figure 5

Figure 5: 2 degree-of-freedom mass spring model

The physical laws will now become a system of equations instead of the scalar one above. The forces acting on *m*_{1} is the spring *k*_{1} dragging in negative direction and *k*_{2} dragging in the positive direction.

*m*_{1}*x*¨_{1} = *−**k*_{1}*x*_{1} + *k*_{2}(*x*_{2} *− **x*_{1})

This is symmetric, and we have an analogue expression for *m*_{2}. The system can be written in matrix form as

When doing continuum mechanics, it is the exact same idea, but the actual equations differ some. Instead of discrete equations, we have continuous functions in space and the governing equations are

semi-discretization yields the following system of equations

with the usual stiffness and mass matrix

**e)**

**e)**

Build the 2d mass matrix as given above.

We are now going to search for solutions of the type:

*u**e** ^{ωit}* (13)

which inserted into (12) yields

*ω*^{2}*M** *** u **=

*A*

**(14)**

*u***f )**

**f )**

Equation (14) is called a generalized eigenvalue problem (the traditional being with *M *=*I*). Find the 20 first eigenvalues *ω*_{i}* *and eigenvectors *u*_{i}* *corresponding to this problem.

**g)**

**g)**

Let *x*_{0} be your initial geometric description (the nodal values). Plot an animation of the eigenmodes by

** x **=

*x*_{0}+

*α*

*u*

_{i}*sin(*

*t*)

You may want to scale the vibration amplitude by some visually pleasing scalar *α*, and choose the time steps appropriately. Note that for visualization purposes, you will not use the eigenfrequency *ω*_{i}* *since you are interested in viewing (say) 1-5 complete periods of the vibration, but for engineering purposes this is a very important quantity.**编程作业代写**

**h)**

**h)**

An interesting question is with regards to ”harmonic” frequencies. The 1st nonzero fre- quency is often called the fundamental frequency, and the rest is called overtones. If the overtones are multiplies of the natural frequency (i.e. *f*_{i}* *= *nf*_{0}, where *n *is an integer and *f*_{0} is the fundamental frequency) the sound is said to be harmonic. This is an important part in all musical instruments. More information can be found in the Wikipedia article on pitch.

**Exercise 4: The best**** ****answer**

**Exercise 4: The best**

**answer**

In this problem you are tasked to give the most accurate finite element code you can make. It is going to be a technical problem, and will focus on improving your finite element library. You are going to compute a finite element solution to problems with known exact solutions

*u *and the goal is to get the error *“**u **− **u*_{h}*” *as small as possible. From the theory we know that *“**u **− **u*_{h}*” ≈ **Ch** ^{p}*, so there are two ways of improving the error

-decrease*h*

-increase*p*

Both strategies will require major revisions to your existing library from **part 1 **of the project; either from implementing new basis functions (*p*-refinement), or optimizing exist- ing algorithms to handle large number of unknowns (*h*-refinement).

**a) ****Profile your**** ****code**

**Profile your**

**code**

Read about *profiling *at the Python documentation page. Use some profile library to find out which part of your program is the slowest. 编程作业代写

**b) ****Compute the**** ****error**

**Compute the**

**error**

To compute the error of you finite element solution, you will need to integrate this in some norm. One usually measures the error of finite element problems in one of the three norms,

The *H*^{1}-seminorm is usually denoted as the energy norm, and satisfies *a*(*u, u*) = *|**u**|*2 1 . From the Galerkin orthogonality, we know that the finite element solution is the best possible solution (in our solution space) as measured in the energy norm. It is the most natural, and predictive way of computing the error.

Compute your error in energy norm on the problem from task 2 and task 3 in **part 1 **of the project. To compute the continuous integral of the norms, you will need to split it into a sum of integration over single elements, and evaluate the error functions here. You will find your assembly procedure from task 2a to be helpful.

**c) ****Get convergence**** ****plots**

**Get convergence**

**plots**

Again, we will use the analytical solutions from part I. Compute your solution on a series of meshes of different sizes. Plot the error you are getting vs the element size *h** *in a log-log plot. What do you see? Explain your findings.

**d)**

**d)**

When the analytical solution is not known, a way to compute the convergence rate is to compare the numerical solution with the numerical solution computed using very small elements.

Compute the solution for element sizes *h **∈ {**, i *= 0*, *1*, *2*, *3*} *and compare these with the solution using element size *h *= . Compute the relative error and the convergence rate.

**e) ****Machine**** ****precision**

**Machine**

**precision**

Computers operate with a limited number of decimal digits. For your typical computing environment, this is of the order *O*(10^{−}^{15}). Can you create a 2D finite element solution which reaches machine precision? How large system did you need, and how long time did the computations take? You can use the time-module in Python to measure this

Figure 6: The L shape

**f) ****) The**** ****L-shape**

**) The**

**L-shape**

Let the domain Ω = [*−*1*,** *1]^{2}*/*[0*,** *1]^{2} be the L-shaped domain around the origin (see Figure 6). Solve the Poisson problem

With the known exact solution

compute *g *from (15) and run your finite element program on this problem. How fine solution are you able to obtain? How long did your simulation take?

更多代写：monash代码查重 澳洲pte代考 Accounting代考价格 cover letter怎么写 literature review點寫 论文conclusion