- 23,708
- 5,924
There are only two unknowns: the two displacements u and v.pyroknife said:Yes, don't we have 4 unknowns in $$\frac{du}{dx} \\ \frac{du}{dy} \\ \frac{dv}{dx} \\ \frac{dv}{dy} $$? And not enough knowns?
There are only two unknowns: the two displacements u and v.pyroknife said:Yes, don't we have 4 unknowns in $$\frac{du}{dx} \\ \frac{du}{dy} \\ \frac{dv}{dx} \\ \frac{dv}{dy} $$? And not enough knowns?
Oops I deleted that post after reading your post prior to the quoted one.Chestermiller said:There are only two unknowns: the two displacements u and v.
No. You are going to be solving the two differential equations for u and v in post #22, subject to the boundary conditions on u and v.pyroknife said:Oops I deleted that post after reading your post prior to the quoted one.
I didn't realize that you can't specify ##\sigma_{yy}## on that face. So yes, I agree that u&v are the unknowns, but to obtain u&v, wouldn't I need to obtain the displacement gradients, and then apply the finite difference formula to compute u&v?
Yes, but the issue is when I discretize the domain into cells, there are issues when I get to cells that lie on the boundary. I plan to use a central difference for the second derivatives, so I would need to know the u&v values in ghost cells that lie outside the real domain. So before I can compute anything, I would need to solve for the u&v values in all the ghost cells given the BCs. I am applying the stress BCs on the face between the boundary cell and the ghost cell, and the displacement BCs are applied directly in the ghost cells. So, for the stress BCs, I do not know the u&v values in the ghost cell. Wouldn't I need to compute the u&v values in the ghost cells first using the stress BC?Chestermiller said:No. You are going to be solving the two differential equations for u and v in post #22, subject to the boundary conditions on u and v.
No. Now that we have everything on the boundaries expressed in terms of the displacements and their derivatives, we no longer have to be working in terms of the stresses. This is where having some experience with solving PDEs numerically comes in. I can help with this. The only problem is properly handling the BCs in conjunction with the finite difference approximations to the differential equations at x = L and y = L.pyroknife said:Yes, but the issue is when I discretize the domain into cells, there are issues when I get to cells that lie on the boundary. I plan to use a central difference for the second derivatives, so I would need to know the u&v values in ghost cells that lie outside the real domain. So before I can compute anything, I would need to solve for the u&v values in all the ghost cells given the BCs. I am applying the stress BCs on the face between the boundary cell and the ghost cell, and the displacement BCs are applied directly in the ghost cells. So, for the stress BCs, I do not know the u&v values in the ghost cell. Wouldn't I need to compute the u&v values in the ghost cells first using the stress BC?
Thanks for all the help. I will type it out in the next post here in a little bit.Chestermiller said:No. Now that we have everything on the boundaries expressed in terms of the displacements and their derivatives, we no longer have to be working in terms of the stresses. This is where having some experience with solving PDEs numerically comes in. I can help with this. The only problem is properly handling the BCs in conjunction with the finite difference approximations to the differential equations at x = L and y = L.
It's late now, so I'll be back tomorrow. Meanwhile, please write out finite difference approximations to the differential equations at the interior grid points.
Nice job.pyroknife said:The central difference formulas are:
\begin{equation}
\frac{\partial ^2 u}{\partial x^2} \approx \frac{u_{i+1,j}-2u_{i,j} + u_{i-1,j}}{\Delta x^2}
\end{equation}
\begin{equation}
\frac{\partial ^2 u}{\partial x \partial y} \approx \frac{u_{i+1,j+1}-u_{i+1,j-1} - u_{i-1,j+1} + u_{i-1,j-1}}{4\Delta x \Delta y}
\end{equation}
Plugging these into the PDEs yield:
\begin{equation}
(2\mu + \lambda) \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{\Delta x^2} + (\mu+\lambda)\frac{v_{i+1,j+1} - v_{i+1,j-1} - v_{i-1,j+1} + v_{i-1,j-1}}{4\Delta x\Delta y} + \mu \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{\Delta y^2} = 0
\end{equation}
\begin{equation}
(2\mu + \lambda) \frac{v_{i,j+1} - 2v_{i,j} + u_{i,j-1}}{\Delta y^2} + (\mu+\lambda)\frac{u_{i+1,j+1} - u_{i+1,j-1} - u_{i-1,j+1} + u_{i-1,j-1}}{4\Delta x\Delta y} + \mu \frac{v_{i+1,j} - 2v_{i,j} + v_{i-1,j}}{\Delta x^2} = 0
\end{equation}
these equations are valid from i,j=1,N. ##\Delta x## and ##\Delta y## are the spatial step size between grid points and are constant.
Assuming that if i=1 or j=1 or i=N, j=N is the boundary node. e.g., i = 0=>x=0, j=0=>y=0,i=N=>x=L,j=N=>y=L. Assuming that the displacement boundary conditions are applied at i = 0 and j = 0, which lie off the real domain, i.e., ##x=-\Delta x##, ##y=-\Delta y##. The stress BCs are applied half way in between i,j = N+1 and i,j=N. i.e., ##x = L + \Delta x/2##, ##y=L + \Delta y/2##
In the above finite difference implementation issues would arise at i or j = N, since we do not know the u or v values at i or j = N+1.
I agree! However, n this scenario, aren't the the stress BCs being applied at x,y=N? If so, then in Eqns (9) and (10), aren't the u&v values unknown at N?Chestermiller said:Nice job.
HANDLING THE BOUNDARY CONDITION AT x= L
We set the grid up so that there are grid points along the boundaries x = L and y = L. We need to derive appropriate finite difference relationships that apply at these boundaries so that we can determine the displacements on these boundaries. At x = L, we showed that the boundary conditions are as follows:
$$\frac{\partial v}{\partial x}=-\frac{\partial u}{\partial y}\tag{1}$$ and
$$\frac{\partial u}{\partial x}=-\frac{\sigma_1}{(2\mu+\lambda)}-\frac{\lambda}{(2\mu+\lambda)}\frac{\partial v}{\partial y}\tag{2}$$
With no loss of generality, we can take the partial derivatives of these equations in the direction tangent to the boundary (i.e., the y direction) to obtain:
$$\frac{\partial^2 v}{\partial x\partial y}=-\frac{\partial^2 u}{\partial y^2}\tag{3}$$ and
$$\frac{\partial^2 u}{\partial x\partial y}=-\frac{\lambda}{(2\mu+\lambda)}\frac{\partial^2 v}{\partial y^2}\tag{4}$$
If we substitute these relationships into the stress equilibrium equations, we obtain (specifically for the boundary at x = L):
$$ \frac{\partial^2 u}{\partial x^2} + \frac{\lambda}{(2\mu +\lambda)} \frac{\partial^2 u}{\partial y^2} = 0\tag{5}$$
$$ \frac{\partial^2 v}{\partial x^2}+\left(2+\frac{\lambda}{(2\mu +\lambda)}\right)\frac{\partial^2 v}{\partial y^2} = 0\tag{6}$$
If we follow the usual procedure of introducing a fictitious grid point at ##x=L+\Delta x##, then we can use the boundary conditions to obtain finite difference approximations to the dependent variable values at the fictitious grid point:
$$u_{N+1,j}\approx u_{N-1,j}-2\Delta x\frac{\sigma_1}{(2\mu +\lambda)}-\frac{\lambda}{(2\mu +\lambda)} \frac{\Delta x}{\Delta y}(v_{N,j+1}-v_{N,j-1}) \tag{7}$$
$$v_{N+1,j}\approx v_{N-1,j}-\frac{\lambda}{(2\mu +\lambda)} \frac{\Delta x}{\Delta y}(u_{N,j+1}-u_{N,j-1})\tag{8} $$
If we substitute these into our finite difference approximations for ##\partial^2 u/\partial x^2## and ##\partial^2 v/\partial x^2##, we obtain:
$$\frac{\partial^2u}{\partial x^2}\approx \frac{2(u_{N-1,j}-u_{N,j})}{(\Delta x)^2}-\frac{2}{\Delta x}\frac{\sigma_1}{(2\mu +\lambda)}-\frac{\lambda}{(2\mu +\lambda)}\frac{(v_{N,j+1}-v_{N,j-1})}{\Delta y} \tag{9}$$
$$\frac{\partial^2v}{\partial x^2}\approx \frac{2(v_{N-1,j}-v_{N,j})}{(\Delta x)^2}-\frac{\lambda}{(2\mu +\lambda)}\frac{(u_{N,j+1}-u_{N,j-1})}{\Delta y} \tag{10}$$
Comments?
Sure. We are solving for them along with the rest of ythe u's and v's at all the other grid points.pyroknife said:I agree! However, n this scenario, aren't the the stress BCs being applied at x,y=N? If so, then in Eqns (9) and (10), aren't the u&v values unknown at N?
I should have mentioned this early, but I am solving the transient equations, which includes the time derivative term. So instead of the 2 PDEs in post #22 being set equal to 0 on the RHS, they are instead set equal to ##\rho \frac{\partial^2 u}{\partial t^2}##. In this case, I think I would need to know the u&v values at x&y=N?Chestermiller said:Sure. We are solving for them along with the rest of ythe u's and v's at all the other grid points.
I don't think so. What is the nature of the transient loading?pyroknife said:I should have mentioned this early, but I am solving the transient equations, which includes the time derivative term. So instead of the 2 PDEs in post #22 being set equal to 0 on the RHS, they are instead set equal to ##\rho \frac{\partial^2 u}{\partial t^2}##. In this case, I think I would need to know the u&v values at x&y=N?
I am working on a research project for modeling thermal protection systems that involves 2 major components: (1) Thermochemical physics modeling and to a lesser extent (2) Structural mechanics modeling (linear elastic should be sufficient). The thermochemical modeling involves solving transient solid decomposition (arrhenius rate equation), continuity, and mixture-energy equations. I am working with a class of materials where the thermo properties such as thermal conductivity, porosity, permeability are expected to be a function of stress/strain.Chestermiller said:I don't think so. What is the nature of the transient loading?
You would be solving the steady state equation if things aren't changing with time. Otherwise, you need to include the time-dependence.pyroknife said:I am working on a research project for modeling thermal protection systems that involves 2 major components: (1) Thermochemical physics modeling and to a lesser extent (2) Structural mechanics modeling (linear elastic should be sufficient). The thermochemical modeling involves solving transient solid decomposition (arrhenius rate equation), continuity, and mixture-energy equations. I am working with a class of materials where the thermo properties such as thermal conductivity, porosity, permeability are expected to be a function of stress/strain.
I'm not too familiar with structural mechanics, but when would I be solving the steady state as opposed to the transient equation?
Aren't the displacement fields always changing with time?Chestermiller said:You would be solving the steady state equation if things aren't changing with time. Otherwise, you need to include the time-dependence.
No way. If you exert a tensile force on a rod, is the displacement field always changing? If you hang a mass from a spring (at a point that it does not oscillate), is the displacement field of the points along the spring always changing?pyroknife said:Aren't the displacement fields always changing with time?
No. Hmmm. So if the boundary conditions are not a function of time, does that essentially mean the displacement field will not change?Chestermiller said:No way. If you exert a tensile force on a rod, is the displacement field always changing? If you hang a mass from a spring (at a point that it does not oscillate), is the displacement field of the points along the spring always changing?
You are familiar with the term "static equilibrium," correct? We covered this in freshman physics. You also probably had a course in Statics and Dynamics, which again covered static equilibrium. Do you see how this relates to what you are asking?pyroknife said:No. Hmmm. So if the boundary conditions are not a function of time, does that essentially mean the displacement field will not change?
Where is the mass going?Also, technically the density of the material I am modeling would change with time because as I mentioned the thermophysics part of the code models solid decomposition. The level of decomposition (level of density change) would depend on the material being modeled. So this could be a mechanism that could make the structural mechanics transient.
So you are including transient heat conduction in your code? If so, then, yes, thermal strain makes the stresses and elastic strains transient. If the temperature is at steady state, however, they you are back to static equilibrium.Furthermore, I am also thinking about incorporating thermal strain into the structural computations, and thermal strain is a function of temperature, which is time-independent in the thermophysics part of the code.
Yes. If you are going to be including this, then you need to familiarize yourself with the subject of poroelasticity. You need to learn about the Terzaghi effective stress, and how to describe the coupling between the fluid flow and the solid mechanics.In addition, the internal pressure of the material may also be changing due to internal gas flow (the materials I am working with are usually highly porous), and this would change the stress I presume.
Yes I know what static equilibrium means. Essentially some body at rest in which net forces sum to 0. I took a statics and dynamics class long ago. I'm still kind of confused on applied loading. With your tensile force applied to a rod example, if you apply this constant tensile force for 10minutes, wouldn't the displacement field within the rod change for the duration the force is applied?Chestermiller said:OMGYou are familiar with the term "static equilibrium," correct? We covered this in freshman physics. You also probably had a course in Statics and Dynamics, which again covered static equilibrium. Do you see how this relates to what you are asking?
Where is the mass going?
So you are including transient heat conduction in your code? If so, then, yes, thermal strain makes the stresses and elastic strains transient. If the temperature is at steady state, however, they you are back to static equilibrium.
Yes. If you are going to be including this, then you need to familiarize yourself with the subject of poroelasticity. You need to learn about the Terzaghi effective stress, and how to describe the coupling between the fluid flow and the solid mechanics.
From what you describe above, I am beginning to sense a geophysical application, involving subsurface flow of liquids, and the effect of increased pore pressure on the solid mechanics (and the feedback of the solid mechanics on the storativity of the rock).
No. The solid would respond virtually instantaneously.pyroknife said:Yes I know what static equilibrium means. Essentially some body at rest in which net forces sum to 0. I took a statics and dynamics class long ago. I'm still kind of confused on applied loading. With your tensile force applied to a rod example, if you apply this constant tensile force for 10minutes, wouldn't the displacement field within the rod change for the duration the force is applied?
This all sounds very interesting. I'm a ChE, so this is the kind of thing I respond to.The material is made up of a solid matrix where its pores are resin infused with some resin. At high temperatures, the resin decomposes and diffuses through the pores into the boundary layer. In addition, the solid matrix also ablates at the surface, resulting in mass law at the surface and a change in the computational domain.
Transient heat conduction is included. At some point the temperature does become steady state obviously.
OK. This helps. I don't think the pore pressures would be high enough to feed back into the solid mechanics.I don't know if I want to go too in-depth with the structural modeling as my focus is in thermophysics&fluids and I don't know very much about structural mechanics at all (and frankly don't have nearly as much interest in it compared to thermophysics and fluids). Poroelasticity sounds like it would primarily involve microscale modeling, which could be interesting and microscale modeling is definitely something I think I will need to do in the future. I just briefly google'd poroelasticity.
"Modeling poroelasticity requires the coupling of two laws. The first of these is Darcy's law...The second law is the structural displacement of the porous matrix. Biot poroelasticity describes this coupled physics." I'm currently using Darcy's law to compute the velocity of the porous gas flow in the continuity equation. And the linear elastic equation solves for displacement.
Cool!This is actually an atmospheric entry application for space exploration. A lot of the thermal protection systems (e.g., PICA from NASA) found on modern atmospheric entry vehicles/probes are resin-infused.
I'm still confused about the tensile loading. You mentioned that the solid would respond virtually instantaneously; are you implying that after this instantaneous response, the displacement would no longer be changing even though the constant load is still being applied?Chestermiller said:No. The solid would respond virtually instantaneously.
This all sounds very interesting. I'm a ChE, so this is the kind of thing I respond to.
OK. This helps. I don't think the pore pressures would be high enough to feed back into the solid mechanics.
Cool!
Sure.pyroknife said:I'm still confused about the tensile loading. You mentioned that the solid would respond virtually instantaneously; are you implying that after this instantaneous response, the displacement would no longer be changing even though the constant load is still being applied?
You've been holding out on me. What department? Do you follow the football team? I've become obsessed with Michigan football. Do you ever go down and hang out at the athletic campus (just to catch a glimpse of the team and Jim)?By the way, I am a UM student.
Hmm maybe that's where I'm confused. In thermophysics, if I apply a constant heat flux, that heat flux would continue to cause changes for the duration it is applied. I don't see why a constant compression load would not continue to cause displacement changes as long as it is being applied?Chestermiller said:Sure.
You've been holding out on me. What department? Do you follow the football team? I've become obsessed with Michigan football. Do you ever go down and hang out at the athletic campus (just to catch a glimpse of the team and Jim)?
I understand what you are asking. And, in an ideal theoretical system, what you are saying would prevail. But, the real world doesn't work like that, and the equations for an elasitic medium with mass/inertia is only an idealization. In the real world there are also dissipative forces at work, and, even though, in many situations, these force are small, they have a cumulative effect and allow the systems we are discussing to relatively rapidly reach static equilibrium.pyroknife said:Also, I was thinking, the equation of motion is a hyperbolic PDE. When you compress a rod, you sent some finite traveling in the direction of compression, which takes time to reach the other end of the rod, leading to changes/oscillations over time. So I don't understand how the response is instantaneous?
Thanks a lot for the examples. Very helpful.Chestermiller said:I understand what you are asking. And, in an ideal theoretical system, what you are saying would prevail. But, the real world doesn't work like that, and the equations for an elasitic medium with mass/inertia is only an idealization. In the real world there are also dissipative forces at work, and, even though, in many situations, these force are small, they have a cumulative effect and allow the systems we are discussing to relatively rapidly reach static equilibrium.
When I said the response of the rod is instantaneous, what I really meant was that the response is virtually instantaneous. This is because the stress disturbance travels at the speed of sound (predicted by the dynamic equations you derived, with inertia included). So for a rod only a foot or two long, the response is virtually instantaneous.
You were talking about the transient temperature effect and comparing it to the transient stress effect. First of all, the transient heat transfer equations are parabolic, rather than hyperbolic. So the disturbance front is diffuse rather than sharp. Secondly, the speed at which the diffusive front progresses is much slower than the speed of sound, because the thermal diffusivity of the material is low. So, if you put a turkey in the oven, it takes hours for the thermal wave to reach the center of the turkey.
Here are two examples to help understand what I'm saying.
Consider a mass suspended on a spring, and you pull the mass downward, and then release it. Theoretically, the mass will oscillate up and down forever at the exact same amplitude, and the stresses and strains in the spring will continue varying forever. But, in the real world, we know that the amplitude of the oscillation will fairly rapidly decrease in time, and the mass will soon come to a new static equilibrium position. This is because of air drag acting on the mass which dissipates the kinetic energy.
Next, consider the kitchen chair I am sitting on right now. This chair is wooden, and has wooden legs. The legs of the chair are in compression, but there are no dynamic oscillations occurring in the legs of the chair. The stresses and strains in the chair legs are constant. This is due to the damping provided to the chair by the kitchen floor on one end, and my fat butt on the other end. This damping caused the oscillations to die out. So the chair and its legs are in static equilibrium. How long did it take to reach static equilibrium? The response to my sitting down was virtually instantaneous.
In analyzing a complicated system such as yours, you need to know the kinds of approximations you can make to get to your answer in a finite amount of time (i.e., real life time: you don't want to be at the U of M forever, right), while at the same time getting very accurate answers. In this case, you are considering the treatment you should use for the solid deformation mechanics part of your calculation. And you already know that another group has handled the solid mechanics by assuming quasi steady state.pyroknife said:Thanks a lot for the examples. Very helpful.
The stuff I do is all computational, and so I lean more towards numerically solving the theoretical equations rather than making sense of it from a more practical aspect. This is just the base and to build on that would be to include/modify certain things in the modeling that would represent less than ideal things more representative of experiments or more realistic scenarios.There is another academic group doing similar work and in their recent paper they solved the steady state equation of motion, claiming that the structural equations reach equilibrium MUCH faster than the thermophysics, which is consistent what you said. I just didn't really understand it.
So this leads to my next question: Would it more representative of a real scenario to solve the transient equations WHILE accounting for phenomena such as damping, air resistance, etc... as opposed to solving a time-independent equation?
The code I'm currently working with doesn't solve for the external flow field, although that may be something for the future. This code only models the material itself. The full blown Navier Stokes equations are not solved for. The momentum equation is reduced to Darcy's law for the internal porous flow, which is a very simplified representation. However, we do have a validated hypersonics CFD code for solving the external flow field, which solves the solves the transient Navier Stokes equations accounting for vibrational nonequilibrium and chemical reactions. THere's one guy in my group working on adding turbulence models to the code. We're not the conventional continuum fluids group since almost everything we model is in a very high temperature regime (~10,000K), where nonequilibrium effects and chemical reactions need to be considered.The heat transfer equation (neglecting radiative heat transfer) that we solve for in the material is:Chestermiller said:In analyzing a complicated system such as yours, you need to know the kinds of approximations you can make to get to your answer in a finite amount of time (i.e., real life time: you don't want to be at the U of M forever, right), while at the same time getting very accurate answers. In this case, you are considering the treatment you should use for the solid deformation mechanics part of your calculation. And you already know that another group has handled the solid mechanics by assuming quasi steady state.
I will call your attention to the fact that you have already made numerous approximations in other parts of your development without subjecting them to the intimate scrutiny you seem to be applying to the solid mechanics. In the aerodynamics, for example, are you solving the full time dependent Navier Stokes equations, including time dependence and turbulent stresses, or are you using boundary layer approximations? In the heat transfer, are you including conduction in the thickness direction and in the flow direction, or are you neglecting heat transfer in the flow direction? Are you including transient time dependence in the heat conduction equation? In applying Darcy's law, are you including time dependence (storativity), or are you assuming quasi steady state flow?
In handling the solid mechanics, you need to consider what your main objective is in this area. Are you analyzing the solid mechanics to include its coupling with the gas dynamics, or are you analyzing it just to determine the state of stress in the solid so you can evaluate strength and failure? It might even be better to use Strength of Materials approximations to the solid mechanics rather than your current Theory of Elasticity approach. So what I'm saying is that it might be advisable to simplify the solid mechanics even further, rather than to consider getting more complicated by including viscous damping explicitly in the equations and solution. What did the other group do? Theory of Elasticity or Strength of Materials? Using a strength of material approach would make life so much easier for you.
You have a system where you have air within the pores providing small scale viscous damping, and external air providing exterior damping. In addition, the shield material itself is probably not ideally elastic, and probably includes some inherent damping of its own (i.e., viscoelasticity). So worrying about going quasi static on the shield material doesn't make sense to me. And thinking about including damping explicitly in the equations really really doesn't make sense to me.
Where is your adviser in all this? Doesn't he/she help keep you pointed in the right direction by pointing out the kinds of accurate approximations you can/should make to simplify your analysis, while at the same time getting results of the desired accuracy?
I have nothing to add, except to say simplify as much as possible. But please check that heat balance equation. It seems to be missing a derivative on the conduction term.pyroknife said:The code I'm currently working with doesn't solve for the external flow field, although that may be something for the future. This code only models the material itself. The full blown Navier Stokes equations are not solved for. The momentum equation is reduced to Darcy's law for the internal porous flow, which is a very simplified representation. However, we do have a validated hypersonics CFD code for solving the external flow field, which solves the solves the transient Navier Stokes equations accounting for vibrational nonequilibrium and chemical reactions. THere's one guy in my group working on adding turbulence models to the code. We're not the conventional continuum fluids group since almost everything we model is in a very high temperature regime (~10,000K), where nonequilibrium effects and chemical reactions need to be considered.The heat transfer equation (neglecting radiative heat transfer) that we solve for in the material is:
$$\frac{\partial (\rho e)}{\partial t} + \frac{\partial}{\partial x_i}(\phi_g\rho_g h_g u_{g_i}) - \kappa_{ij}\frac{\partial T}{\partial x_j} = 0$$
Removing that middle convection term (which accounts for energy transfer through convection of the internal gas) would yield the classical transient heat conduction equation. The materials we work with are typically anisotropic, so thermal conductivity is a 9 element tensor and accounts for the various directions.
In our modeling, Darcy's law is just the steady state form.
As far as I am aware of, there are 2 other groups trying to do thermophysics-structural modeling WITHIN the material. Both of them use the time-INDEPENDENT theory of elasticity. I think solving for the linear elastic equation is pretty simple? Although, I didn't really understand how the boundary conditions worked, but I think I have a much better idea now. Also, I'm trying to develop the software such that it can be extendable in the future for others who want to add more complex structural physics, but for me, I think my priority is just to get a stress and displacement field. The stress field is primarily for modeling the stress-dependency of the thermophysical material properties (for example conductivity in the above equation, porosity ##\phi## in the above equation) and the displacement field is primarily for moving the mesh.Also, my adviser isn't a structural person. The people providing my funding sort of steered me in this direction and these people are also not structural people. Plus, I have been working out of state this year. I am kind of on my own with the structural side.
Thank you for all the help.Chestermiller said:I have nothing to add, except to say simplify as much as possible. But please check that heat balance equation. It seems to be missing a derivative on the conduction term.