Dr A. Cookson Department of Mechanical Engineering

2020/2021 University of Bath

Page 1 of 4

Assignment 2: Transient MATLAB-Based FEM Modelling

ME40064: System Modelling & Simulation

ME50344: Engineering Systems Simulation

You must extend your finite element code to be able to solve the transient form of the

diffusion-reaction equation, Eq. 1, using either the backward Euler or Crank-Nicolson time

integration scheme:

=

2

2

+ +

Eq. 1

You will then use this code to model the effects of heat damage to human skin, thereby being

able to estimate the level of thermal shielding that a piece of protective clothing must provide.

Part 1: Software Verification [35%]

Demonstrate that your code is working correctly by solving the transient diffusion equation,

Eq. 2, for the domain x = 0 to 1, subject to the following initial and Dirichlet boundary

conditions, and compare to the analytical solution in Eq. 3:

=

2

2

Eq. 2

= [0,1], (, 0) = 0, (0, ) = 0, (1, ) = 1, > 0

In order to ensure an accurate solution, it is suggested that you use an element size of 0.1 (i.e.

a 10-element mesh) and a time step, Δt = 0.01.

As shown in Lecture 13, this problem has the analytical solution in Eq. 3, a MATLAB function

for which, TransientAnalyticSoln.m, is provided on Moodle:

(, ) = +

2

∑

(−1)

∞

=1

−

22 sin()

Eq. 3

Create the following two figures to demonstrate that your code is correct:

a. Plot your solution c(x) vs. x, showing the solutions at t = 0.05, 0.1, 0.3, 1.0, in the format

shown in Lecture 13.

b. Plot both the analytical solution and your numerical solution at x=0.8, for

t = 0 to 1.0.

Dr A. Cookson Department of Mechanical Engineering

2020/2021 University of Bath

Page 2 of 4

Part 2: Modelling & Simulation Results [35%]

Now use your code to solve the tissue burn model introduced in Lecture 14. The temperature

distribution in tissue is modelled by Eq. 4:

= (

)

2

2

− (

) + (

)

Eq. 4

Your finite element mesh and material parameters must represent the three-layered structure

of skin, as shown in Figure 1.

Figure 1: Schematic of finite element mesh required to represent skin tissue layers

Parameter values for each layer are provided in Table 1 – with the exception that in Question

1 you should assume zero blood flow i.e. G=0 everywhere. The layers of the tissue are

defined at the following x coordinates: E=0.00166667, D=0.005, B=0.01.

Throughout you should discuss the meaning and accuracy of your results in terms of the

physics, initial conditions, boundary conditions, numerical methods, and modelling

assumptions. You should consider whether your results are converged in both space and time,

and therefore the validity of your conclusions. You may wish to explore differing mesh

densities for each layer, and you should explain how you represent the layers in your code.

1. Run your code for the following initial conditions and Dirichlet boundary conditions, for a

maximum of 50 seconds:

(, 0) = 310.15, ( = , ) = 310.15, ( = 0, ) = 393.15

Plot the temperature distribution through the tissue for a range of time-points between 0-

50 seconds and discuss the results.

2. Use your solution to determine the level of tissue damage that will occur by numerically

evaluating the integral in Eq. 5 at x=E and x=D:

Γ = ∫ 2 × 1098

=50

exp (−

12017

( − 273.15)

)

Eq. 5

Integrate between the time points, tburn, (at which the temperature T becomes greater than

317.15 Kelvin), and t=50. Note that the MATLAB function, trapz(x), assumes an interval

of 1, therefore multiply the output by timestep, Δt, to obtain the final value of the integral.

Dr A. Cookson Department of Mechanical Engineering

2020/2021 University of Bath

Page 3 of 4

If Γ > 1 at x=E, this produces a second-degree burn. State the value of Γ you have

calculated at this position (don’t be alarmed if it is very large or small), and hence whether

you expect a second-degree burn to occur.

If Γ > 1 at x=D, this produces a third-degree burn. State the value of Γ you have calculated

at this position (don’t be alarmed if it is very large or small), and hence whether you

expect a third-degree burn to occur.

3. Determine the minimum temperature reduction (to a reasonable level of precision) that

must be achieved by the protective clothing at the outer boundary, x=0, in order that

second-degree and/or third-degree burns do not occur within 50 seconds.

State the final Dirichlet boundary condition at x=0 that achieves this, and explain how you

used your code, and any other calculations, to estimate this value. Include any data or

figures that you feel are relevant.

4. Re-run your results for question 2, but now including the blood-flow related reaction &

source terms, for the values given in Table 1. How much does this change the temperature

distribution in space and time? Discuss whether this is an important effect to consider in

future models of this problem.

Parameter Epidermis Dermis Sub-cutaneous

k 25 40 20

G 0 0.0375 0.0375

ρ 1200 1200 1200

c 3300 3300 3300

ρb - 1060 1060

cb - 3770 3770

Tb - 310.15 310.15

Table 1: Parameter values for Part 2. Note that some are realistic and others are not,

having been altered to allow you to solve your model in a reasonable time frame.

Additional Marking Criteria

Presentation [5%]

• Clear graphical presentation of equations, text, diagrams, and plots

• Standard of written English and explanations

• Proper use of references, figure numbering & captions, and table headings

• Clear layout of report

Code Quality & Verification [5%]

• Appropriate use of unit tests

• Elegance & generality of code for representing a layered material

• Readability – use of meaningful variable names and code commenting

Dr A. Cookson Department of Mechanical Engineering

2020/2021 University of Bath

Page 4 of 4

Extra Code Features [20%]

None of these extra features are required in order to solve Parts A and B. In order of difficulty

of implementation, they are:

• Testing forward & backward Euler and Crank-Nicolson methods

• Gaussian quadrature for evaluating integrals

• Implementing quadratic basis functions

• Investigating errors & convergence using the L2 norm

SUBMISSION GUIDELINES

You may structure your report as a set of answers to these questions – there is no

requirement to write this in a lab report format. Please note, you do not need to re-explain the

entirety of the finite element method. However, your report must be self-contained and

therefore must not assume that the reader knows the content in this document.

• You must include all your MATLAB source code as text in the Appendices – failure to

do so will cause you to lose marks i.e.

- Do not paste your code into the document as an OLE item or as an image.

- Do not upload archived/zipped/compressed folders of these source files.

• Do not use MATLAB’s symbolic algebra toolbox – you will lose marks.

• Word limit of 2000 words (not including source code).

• You are reminded of the university’s policies regarding academic plagiarism: your code

and your report must be your own work.

Submit your work via the online submission link on the unit Moodle page.

Deadline: 4pm on Friday, 4th December 2020.

欢迎咨询51作业君