# 辅导案例-ME40064-Assignment 2

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.

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

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

Deadline: 4pm on Friday, 4th December 2020.  Email:51zuoyejun

@gmail.com