编程作业代写-TMA4220代写-Programming代写

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 heatequation
– the linear elasticityequation
ρ ∂t2 = ∇σ(u) – the free vibrationequation

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 编程作业代写

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

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 M as defined above.

b) 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

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₊₁ = y_n + hg(t_n, y_n).

More sophisticated include the improved Eulers methods

y_n₊₁ = y_n + (g(t_n, y_n) + g(t_n₊₁, y_n + hg(t_n, y_n)))

or the implicit trapezoid rule

y_n₊₁ = y_n + (g(t_n, y_n) + g(t_n₊₁, y_n₊₁))

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

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₌₀ = u₀(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))²+y²)

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

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 编程作业代写

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 = −νε_xx. This gives the following relations

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

(8)

and some appropriate boundary conditions

a) 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

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

Au = b

where

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

c) Test case

Show that

is a solution to the problem

where

and Ω = {(x, y) : max(|x|, |y|) ≤ 1} is the refereance square (−1, 1)².

d) 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

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⁰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 )

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

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:

mx¨ = −kx

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₁ is the spring k₁ dragging in negative direction and k₂ dragging in the positive direction.

m₁x¨₁ = −k₁x₁ + k₂(x₂ − x₁)

This is symmetric, and we have an analogue expression for m₂. 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)

Build the 2d mass matrix as given above.

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

ue^ωit (13)

which inserted into (12) yields

ω²M u = Au (14)

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)

Let x₀ be your initial geometric description (the nodal values). Plot an animation of the eigenmodes by

x = x₀ + α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)

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₀, where n is an integer and f₀ 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

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

-decreaseh

-increasep

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

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

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¹-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

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)

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

Computers operate with a limited number of decimal digits. For your typical computing environment, this is of the order O(10⁻¹⁵). 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

Let the domain Ω = [−1, 1]²/[0, 1]² 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?

合作平台：essay代写论文代写写手招聘英国留学生代写