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作业君
