Numerical methods to calculate solidification/melting fronts

[Nivee~]

So, it's summer vacations and I have an assignment. It's to write a C program to calculate the solidification/melting front of a calculation domain using the Finite Volume Method. I have no idea where to begin, to be honest.. Can someone please tell me what I should do, how I should proceed?

[kerbiloid]

Probably first try to google
"calculate the solidification/melting front"
at have a look at first 2-3 links.

[Guest]

And when that doesn't work (Google in general being completely worthless for serious research), try looking for a textbook on numerical methods. First thing to understand is always the method itself. Implementing it, even in such a primitive language as C, tends to be comparatively straightforward. I assume you already have the differential equations you need to solve, if you don't, then you'll need something on thermodynamics, as well.

Edited May 28, 2019 by Guest

[kerbiloid]

24 minutes ago, Dragon01 said:

And when that doesn't work (Google in general being completely worthless for serious research), try looking for a textbook on numerical methods.

The first mentioned link was a very large piece of text with charts, formulas, and description.
So, I guessed it's a starting point./

[Guest]

There's no guarantee it'll pop up for him, this being Google. If it decides to show you useless crap, it will only show you useless crap. Even if you disable its "search improvements", results can and will vary across countries and sometimes across computers.

Edited May 28, 2019 by Guest

[kerbiloid]

https://www.sciencedirect.com/science/article/pii/S0307904X05001460

 

[Selective Genius]

Don't use google for looking up serious research articles. Use scholar.google.com

Also, ResearchGate is another useful source of research material. It's like facebook for researchers. Do a thorough literature review of existing research, organize them and maybe, you will find a research gap! You may try to fill in the gap with your own research, should you be inclined to do so.

I would suggest you first understand the basics of using numerical methods for calculating integrals, solutions of ordinary and partial differential equations etc. Solidification/melting problems, if you also want to study the behavior of the liquid part, will incorporate a lot of partial differential equations, each having to be solved numerically.

Also, decide on what scheme you would be using. From what I have observed, Finite Volume Method is the most accurate, but difficult to implement, and Finite Difference Method is the easiest to implement.

MOST IMPORTANTLY,  learn C properly. If you know how to implement file handling in C, you can do so much more than you can imagine. However, take one step at a time.

7 hours ago, Dragon01 said:

primitive language as C

It may be decades old, but there's a reason why it's still around. I would rather call it 'vintage'.

Edited May 28, 2019 by Selective Genius
typso!

[Guest]

There's a reason, and that reason is that it's literally one step from being assembler while actually being somewhat portable. If you want performance, you go with FORTRAN (itself an ancient thing, but supercomputer people love it). If you want ease of use, you go with Wolfram (my own language of choice) or maybe C#. C, with its heisenbug-inducing pointers, linear execution and general lack of many useful abstractions, is only good for writing drivers and learning how the computers work under the hood.

Google Scholar can be handy, but I'm generally wary of using anything by Google for anything that's not exactly in line with the most popular use. Anything by Google tends to be so dumbed down that if it wasn't for sheer convenience, it wouldn't be worth using.

Edited May 28, 2019 by Guest

[IncongruousGoat]

Just now, Dragon01 said:

C, with its heisenbug-inducing pointers, linear execution and general lack of many useful abstractions, is only good for writing drivers and learning how the computers work under the hood.

Wait, what? C's execution is as linear or non-linear as you want it to be, and if pointers are causing Heisenbugs then you need to seriously re-assess your code. Also, there's plenty it's good for other than just drivers. Kernels, operating systems, servers, compilers, interpreters, VMs... Anything that needs to perform low-level interactions quickly and efficiently. It's also not half-bad for supercomputing applications.

@Nivee~ To address the original question... what part of it do you need help with? If it's the numerical analysis, I can't really help you, but if it's help with C you're after then I can answer questions all day.

[Cunjo Carl]

@Nivee~

I had to write one of these about 15 years ago, I'll do my best to remember what's up and get you going in the right direction. Come to think of it, computers are like 20 times faster now. That should be extremely handy. Unfortunately, I've long since forgotten C semantics beyond the basics!

Ok, so to get started, we need to pin down more precisely what your goal is. Then, depending on the specifics we can pick which of the 50 bajillion kinds of FEA (numerical method) we want to use for the job! I'll go over a bit of the basics just to get conversation rolling. I don't know what level you're at, but I'm assuming this is a grad course?

A. Vakrushev, M. Wu, et al. Numerical Investigation of Shell Formation in Thin Slab Casting of Funnel-Type Mold

1. Pick the system we're trying to model

So if you're doing this for a class, the good news is you just need to make good looking results! It doesn't need to predict anything, you can instead massage the program to make results like how you want them to look. But we need to find out what kind of materials we're trying to model. In other words, do we want to see a very simple solidification like in plastic, a somewhat interesting solidification like water+ice, or a more complicated one like the dendritic solidification of metal castings shown above? Specifically, the things we need to nail down:

• Is our model 2D or 3D
• What size/time scale of effects do we want to see
• Do we need to model nucleation
• Is the solid crystaline or amorphous
• Do we need to model different solidification rates on different solid crystal faces
• Do we need to model different thermal conductivities in different directions of the solid
• Do we need to model dissolved materials (like solutes or dopants)
• Do we need to model bulk flow or convection in the fluid (like for injection molding or extremely high temperature differentials)

Since this is for a class if you've never heard about some of these things, then we don't need to worry about them. Hooray! Easier=better.

2. Pick mesh type

The mesh is the shape of the little cells we've divided the world into so we can compute it. We can go with 3 kinds of mesh: Fixed uniform, Fixed, or Dynamic. Then, we can go with 2 basic kinds of shapes: triangles/tetrahedrals, or squares/cubes. If at all possible, we'll go with a fixed uniform mesh of square cells, which makes life so much easier. With one of these, we divide the world into a big square grid, and all the programming and mathematics work out nice and simple. For small/simple systems these are fast to write and fast to run, but they have a really hard time with things like dendritic solidification fronts. Next in line are fixed non-uniform meshes, where you can make smaller cells (= more accurate results) in areas where you think you'll need it, like on sharp corners of a melting ice cube. Finally is the dynamic mesh where the computer automatically makes smaller cells where it detects they're needed. This is hard because it makes information tricky to find in memory- as a side note, how fast your FEA runs has much more to do with how hard it is for the computer to find things in memory than it does with how many commands the computer needs to run. If your program makes the processor 'unexpectedly' look something up from your RAM it can take 75-100 computer cycles to do it, and when it happens a lot this time delay (called cache thrashing) plagues most FEA designs including my nemesis COMSOL. Anyways, if we do simple squares in a fixed uniform grid for a mesh we get to avoid all the hassle. Hooray! Easier-better.

3. Pick governing equations

Specifically, we need an equation of state for the solidification (most simply just the enthalpy of melting), equations for heat and mass flow (most simply just thermal conduction and negligable fluid flows), and equations for materials properties like density at different temperatures (most simply just ignored and assumed constant). We'll probably pick our governing equations based very specifically on what we're trying to model and what mesh we went for. It's always tempting to get really fancy with how to do your governing equations because SCIENCE, but to be honest, it's typically fastest easiest and most effective to just precalculate whatever we need and toss it into a lookup table . Once again, Hooray! Easier=better.

It's a lot, but hopefully enough to start thinking? Always remember that bottom line. Easier=better!

Edited May 30, 2019 by Cunjo Carl

[Nivee~]

14 hours ago, Cunjo Carl said:

I'm assuming this is a grad course?

Yes, Yes I am in a grad course... and while the melting problem that I have been told to solve involves a pure metal (Gallium)(ENTHALPY-POROSITY TECHNIQUE FOR MODELING CONVECTION-DIFFUSION PHASE CHANGE: APPLICATION TO THE MELTING OF A PURE METAL), a universal method applicable for both alloys and pure metals would be great.

This is what I want my program to exhibit:

 

The calculation domain is 2D,rectangular. The material is crystalline, and is assumed to be anisotropic,so different thermophysical values in each direction. It does not have any impurities or solutes, so that's one problem out of the way. Although I would really like to go one step further and include impurities, let's go one step at a time. And yeah, convection in the calculation domain (beyond a certain value of liquid fraction) will also be required.

While the problem linked in the paper above is relatively simple, I want to go a lil further than that..

About modelling solidification process, I came across the Enthalpy Porosity Method and the Effective Heat Capacity method. Which one do you think will be more suitable??

On 5/29/2019 at 5:39 AM, IncongruousGoat said:

To address the original question... what part of it do you need help with? If it's the numerical analysis, I can't really help you, but if it's help with C you're after then I can answer questions all day.

Thanks, you are too kind! But I first need to sort out the numerical foundation of this problem before I venture into C...

[Cunjo Carl]

On 5/30/2019 at 9:11 AM, Nivee~ said:

Yes, Yes I am in a grad course... and while the melting problem that I have been told to solve involves a pure metal (Gallium)(ENTHALPY-POROSITY TECHNIQUE FOR MODELING CONVECTION-DIFFUSION PHASE CHANGE: APPLICATION TO THE MELTING OF A PURE METAL), a universal method applicable for both alloys and pure metals would be great.

This is what I want my program to exhibit:

*snip*

The calculation domain is 2D,rectangular. The material is crystalline, and is assumed to be anisotropic,so different thermophysical values in each direction. It does not have any impurities or solutes, so that's one problem out of the way. Although I would really like to go one step further and include impurities, let's go one step at a time. And yeah, convection in the calculation domain (beyond a certain value of liquid fraction) will also be required.

While the problem linked in the paper above is relatively simple, I want to go a lil further than that..

About modelling solidification process, I came across the Enthalpy Porosity Method and the Effective Heat Capacity method. Which one do you think will be more suitable??

Looks good so far! Gallium is a nice choice because it has a huge energy barrier to nucleation, so we should be safe just assuming the one front. Also nice for Gallium, it doesn't change density much when it melts and doesn't form complex macroscopic shapes when solidifying, which are all big blessings. The model shown appears to be with nice low temperature differentials (10C) on large size scales (2-20cm) and long-ish times (~~1s timesteps.) Also, for larger systems, the trick in this paper became a standard- modeling all the complexities of those dendrites at the solid-liquid interface as just a thin 'mush' with scalable properties. All of this works together to make the square grid a great choice for the materials, time, temp and size scales. I like it.

I've never heard of the two techniques before. That said, I'm happy to take a look! Just looking at how they're derived, I think the enthalpy porosity model is more appropriate to us here, whereas the effective heat capacity model would be more at home with polymers or amorphous materials which have low nucleation energy barriers and wide melting temperature ranges.

I'll give your paper a quick read and see if I can learn more what it's up to.

Next is choosing what we want each loop to look like. We need to calculate a bunch of different physical processes like conduction, convection, phase change, etc and all of these can happen on different time scales. What's more, if we have sharp geometries or very nonlinear governing equations, it's possible to accidentally overshoot what should happen in a given timestep and cause crazy oscillations or divergences. To explain divergences, let's consider KSP. Here think of a rocket glitching and falling through Kerbin. When it touches the center, the model 'diverges' = miss calculates and launches our craft off to infinity at light speed despite this being a thing that wouldn't physically happen. Physically our craft should just get stuck at the center but because the timestep was too long and the physics too nonlinear, the program miscalculated what should happen. Alternatively, a bit more complicated, any time a craft wiggles itself to death while just hanging out in space, that would be an oscillating divergence. We need to deal with the same sorts of problems in our models too. The options for dealing with this are:

Single physics loop with short timesteps: We have all of our physics being run on the same timescale (in the same for loop), allowing for quick execution and clean code. However, we choose timesteps that are short enough to always avoid divergence for our model, which can be wastefully slow. It's typically not the most efficient way to go computationally, but in the cases where it's close enough it's definitely worth it for the simplicity. KSP uses this structure. Because they rely on rapid simple calculations, these setups benefit hugely from precomputed heuristics and lookup tables. My very favorite approach is to use one of these and have the setup calculate the optimal timestep on the fly.

Nested physics loops: Those parts of physics that are more prone to divergence are run with much shorter timesteps than the rest of physics. This requires nested forloops and is much slower computationally. However, if there's only one part of your physics giving divergence problems it's typically faster than slowing down everything. I tend to avoid these though, in general.

Convergence loops: We full on embrace the terrible non-linear nature of our system and for each cell we have a solver (often similar to Newton's method) determine what the values of the next timestep should be. These solvers are specifically tuned (by you) approach their value slowly enough to never have each solver step diverge, but quickly enough it doesn't need to run until next week. Though they're often necessary, these are extremely cumbersome computationally.

 

Given our specific system, I can happily recommend the single physics loop for starters. And though there's a lot more to the details, if we went with the single physics loop the shape would look like:

Create two two arrays called "old" and "new" representing our system as grids of square cells, both identical and with the correct starting conditions
for (each timestep) {
  for (each cell in "old") {
    for (each neighbor) {
      (using Enthalpy/Pressure/etc. values from "old")
      calculate heat flow from conduction and add/subtract the enthalpy from the cells in "new"
      calculate heat flow from convection and add/subtract the enthalpy from the cells in "new"
      calculate heat flow from viscous friction and add/subtract the enthalpy from the cells in "new"
      etc. etc, going through all of our governing equations also including for mass flow.
      (exceptions for boundary conditions)
    }
  }
 for (each cell in "new") {
   Update physical properties (amount of phase change, the temp change from phase change, density change, etc.)
 }
 for (each cell in "new") {
   copy over cell in "old"
 }
 Display
}

Cheers!

[Nivee~]

Well, I am a lil confused.. How  do I mesh the grid? How do I let the system know that a domain has been divided into several discrete finite volumes?

[satnet]

1 hour ago, Nivee~ said:

Well, I am a lil confused.. How  do I mesh the grid? How do I let the system know that a domain has been divided into several discrete finite volumes?

Unless your instructor gave you an existing library it is up to you to create the model. I don't know quite enough thermodynamics to know all the properties you'll need to track (temperature and density are probably a given, but I doubt that is all), but in general you'll need something to store the properties of each cell, in this case probably a C struct. If you are using a simple grid a 2d array of these structs would probably make sense where each index into the array represents your spatial step unit in either the x or y dimension. This is a simple model, but it sounds like your professor gave you a problem where a simple model is sufficient. 

Then as @Cunjo Carl indicated for every temporal step you update the cells in your array representing your new model using you old model array, display it, swap it into your old model array, then start again.

[Cunjo Carl]

20 hours ago, Nivee~ said:

Well, I am a lil confused.. How  do I mesh the grid? How do I let the system know that a domain has been divided into several discrete finite volumes?

Let's start by making something simple, a program to model heat conduction in a 1mx1m uniform square of metal. It starts at some temperature (say 30C) and is coated mostly in insulator. At time=0 we suddenly apply heating and cooling a patch on top/bottom to hold the surface at 50C, and a patch on left/right to hold the surface at 20C. Let's track the temperature over time.

(A preface for others: For the sake of a good explanation, I'm going to fib a few times, and present a way that's computationally quite silly, but easy to start thinking about. Also as a note, I haven't coded in C++ for about a decade, so parts may be off... let's call it good pseudo code!)

To start, imagine our metal is made of a grid of 100x100 of little square 'cells'. Each cell is so small that only the temperature at the center of each one is important. So, at the center of each square cell we can imagine a point, and the that point will determine the properties of everything in the cell.  For this heat conduction program, the only information that needs to be stored is the Temperature at the center of each cell, and everything else (like density, heat capacity, etc.) will just be considered uniform throughout all cells in the metal block. Because our cells are so small, and our time steps are so short, we can use the simplest of heat transfer equations, the steady state heat transfer for linear temperature differences. These are shown in the spoiler and digested into a convenient form for our program.
 

Spoiler

The standard heat conduction laws applied to heating a little 'cube' volume using heat transferred across an identically sized 'cube' block of material. Notice how the area terms cancel, so the actual thickness of our volumes winds up being irrelevant. Both length terms that appear will be the same as our cell size.

Happy Linear Laws
HeatFlow (in Watts) = m*cp*deltaT/time = k*Area*(T2-T1)/(x2-x1)

Substitutions
(x2-x1) = Length
time = timeStep
m = Density*Volume = Density*Area*Length
cp = Cp/MolarMass

Plug them in
Density*Area*Length*(Cp/MolarMass)*deltaT = timeStep*k*Area*(T2-T1)/Length

Simplify
deltaT = timeStep*k*MolarMass/(Density*Cp*Length^2) * (T2-T1)
deltaT = coefficient * (T2-T1)               where coefficient   =   alpha*timeStep/Length^2     ,   and alpha =  k*MolarMass/(Density*Cp)

 

So we need our program to make a matrix of 100x100 variables that can store the temperatures at each of the 100x100 points. We'll do this using what's called a 2D array. Making one is simple in C++!

double T_old[100][100];  //Makes a 100x100 matrix of variables called an 'array'
double T_new[100][100];

/* Each value (called an element) in this array will represent the temperature at the center of one of our squares. Specifically it's stored on the computer as a double. A 'double' is a kind of computer variable that can store numbers with good accuracy, and is fast to use on newer computers. You could also use the very similar kind of variable 'float' which is half the size.  */

// The numbers stored in the arrays are random to start, so first you need to use for loops to set every value to 30 (because we start at a uniform 30C).

int x = 0;
int y = 0;
for (x=0; x<=99; x++) {
    for (y=0; y<=99; y++) {
        T_old[x][y] = 30.00;
        T_new[x][y] = 30.00;
    }
}

//Then we make the rest of our variables:

double alpha = 0.01855 0.00001855; //This is Gallium's thermal diffusivity  EDIT: Missed a few zeroes putting it into SI
double cellSize = 0.01; //The distance from the center of each cell to the one next to it is 1 cm
double timeStep = 0.0001; //Make sure timeStep is ~~ cellSize^2 or smaller to avoid divergence
double coefficient = alpha*timeStep/(cellSize*cellSize); //see spoiler

double time = 0.00;  //We make a variable to track time

for (time=0.00; time<=300.00; time+=timeStep) {  //Starting at t=0 and stepping forward .0001seconds at a time until 5 minutes
    for (x=0; x<=99; x++) {
        for (y=0; y<=99; y++) {
            Trade heat with the neighboring cells following T_new[x][y]  += coefficient * ( T_old[x neighbor][y neighbor] - T_old[x][y] )  doing this for each neighbor to add in all their contributions.
            If we're at a boundry heat/cold supply, trade heat with that boundry, again following        T_new[x][y]  += coefficient * ( BoundryTemp - T_old[x][y] ) 
        }
    }
   // Copy T_new over T_old  ,  just using another pair of for loops and  T_old[x][y]  = T_new[x][y] ;
}
//Display T_new

Edited June 4, 2019 by Cunjo Carl
** off by one coding corrections

[Flibble]

On 5/28/2019 at 10:26 PM, Dragon01 said:

There's a reason, and that reason is that it's literally one step from being assembler while actually being somewhat portable. If you want performance, you go with FORTRAN (itself an ancient thing, but supercomputer people love it). If you want ease of use, you go with Wolfram (my own language of choice) or maybe C#. C, with its heisenbug-inducing pointers, linear execution and general lack of many useful abstractions, is only good for writing drivers and learning how the computers work under the hood.

I use both at work, and would take C over C# every time. C# just feels clunky by comparison, and lets not even get started on performance.

Also heisenbugs are more common from badly constructed threading causing race conditions than pointers in my experience.

[Nivee~]

5 hours ago, Cunjo Carl said:

* incredibly helpful code snippet*

for (x=1; x<=99; x++) {
    for (y=0; y<=99; y++) {
        T_old[x][y] = 30.00;
        T_new[x][y] = 30.00;
    }
}

Hmm, but we are not providing for boundary temperatures here...

How about this?

Spoiler

// for temperatures of interior nodes

for (x=1; x<=98; x++) {
    for (y=1; y<=98; y++) {
        T_old[x][y] = 30.00;
        T_new[x][y] = 30.00;
    }
}

//For the boundaries, left wall is at 85 C, and the other walls are at 30C

for (x=0; ; ) {
    for (y=0; y<=99; y++) {
        T_old[x][y] = 85.00;
    }
}  //left wall

for (x=99; ; ) {
    for (y=0; y<=99; y++) {
        T_old[x][y] = 30.00;
    }
}  //right wall

for (y=0; ; ) {
    for (x=1; y<=98; y++) {
        T_old[x][y] = 30.00;
    }
}  //top wall

for (y=99; ; ) {
    for (x=1; x<=99; x++) {
        T_old[x][y] = 85.00;
    }
}  //bottom wall

Is this OK? Because my program is not running, strangely....

[Flibble]

3 minutes ago, Nivee~ said:

for (x=0; ; ) {

That, and similar loops aren't valid.

Try just doing something like:

x = 0;
for (y=0; y<=99; y++) {
  T_old[x][y] = 85.00;
} //left wall 

 

6 hours ago, Cunjo Carl said:

double T_old[99][99];  //Makes a 100x100 matrix of variables called an 'array'
double T_new[99][99];

Those make 99x99 arrays, you wanted:

double T_old[100][100];  //Makes a 100x100 matrix of variables called an 'array'
double T_new[100][100];

These can be indexed from 0-99 (inclusive).

[Cunjo Carl]

Thanks for the syntax correction, @Flibble, that stuff never sticks for me.

I'm also noticing my for loops on x were all starting on 1 when they should have been starting on 0. I'll just go fix that.

 

3 hours ago, Nivee~ said:

Hmm, but we are not providing for boundary temperatures here...

How about this?

*clip*

Is this OK? Because my program is not running, strangely....

I definitely left a lot to the imagination with the pseudocode! I'm not surprised it won't run by copy pasting. I can help a bit with debugging semantics problems, but it's been so long I'd need to sit down in front of C for an hour to get the syntax back in my head. Mostly, I think you should be lined up in the right direction now though as far as what you're trying to accomplish, and exactly how is now part of the fun challenge. Just as a note, by habit I'm one of those teachers with a mildly sadistic bent who'll purposefully leave hidden hurdles in the way of a problem for people to crash into. Overcoming them is where the best growth happens! When I work with people in person I tend to know how high to size them, but in this case I think I've been going too large! Sorry.

It sounds like this may be your first time coding in C++, and if so you might actually get the fastest results by doing simple silly things to start that have nothing to do with CFD. Just make it say "Hello World" 1000 times, do the Fibonacci sequence, that sort of stuff. It should help you understand the overall mechanics of it, which will be invaluable moving forward. Also, please make sure you have a debugger, at least one which can tell you why the program isn't running like "Expected ___ line 66" So you can have an idea of why it's not running. It's always something silly... In this case, it'll almost definitely be an "Error invalid index, Line XX" where the line refers to " T_new[x][y] += ...  " and when y = 99 for the first time.

Something like that should work for the boundary condition, following Flibble's correction of course. If you use this approach, you'll need to apply those temperatures to the boundary in T_new in each step of the program just before copying T_new to T_old. Otherwise, we'll accidentally overwrite our boundary conditions during the program! I was hinting that you could also just apply the effects of the boundary conditions using if statements within the main code, but this would be way slower computationally, so go with however makes sense to you! I'd recommend not applying it to the entire wall though, just a patch in the middle. Maybe 40-60?

Unfortunately, I've really overdone typing in the last few days, so I'll be needing to mostly tune out for a week. Best of luck moving forward!

 

Edited June 4, 2019 by Cunjo Carl

[Flibble]

7 minutes ago, Cunjo Carl said:

It's always something silly... In this case, it'll almost definitely be an "Error invalid index, Line XX" where the line refers to " T_new[x][y] += ...  " and when y = 99 for the first time.

Using built-in arrays like that you get no bounds checking. So the program will just run, but will corrupt adjacent data in the array (also it might crash due to accessing off the end of the array, but YMMV there).

If you use something with safeties, like std::array then it will throw an exception when you try to access off the end. This is probably a good thing while learning!

I code C++ for a living, so if you have questions post them up and I'll try to answer.

[FleshJeb]

For impurities, would you just have a grid of alphas that you reference by index? (Impurities have a different alpha). You can scatter those randomly (or in whatever proportion) after setting the alpha array uniformly.

You can do the same trick for any other relevant physical property.

How in the hell do you deal with the liquid state, where the impurities would shift cells? (After googling a bit: https://www.cc.gatech.edu/~turk/my_papers/melt.pdf)

[Cunjo Carl]

1 hour ago, Flibble said:

So the program will just run, but will corrupt adjacent data in the array ...

  That's amazing. God I miss C, why did I ever leave? (Joking but also actually serious)

 

[Nivee~]

Assuming the following to be the calculation domain, with blue dots being the boundary nodes, and the yellow dots being the interior nodes,

Well, here's the final code:

Spoiler

#include<stdio.h>

int main()
{
    double T_old[5][5];
    double T_new[5][5];
    double Tin;
    int x = 0;
    int y = 0;
    printf("\n Enter the Initial Temperature");
    scanf("%d",Tin);
    for (x = 1; x < 4; x++){
        for (y = 1; y < 4; y++){
            T_old[x][y] = Tin;
            T_new[x][y] = Tin;
        }
    }

    //Boundary Conditions
    x = 0;
    for (y = 0; y <= 4; y++){
        T_old[x][y] = 85.00;
    }//Left Wall

    x = 4;
    for (y = 0; y <= 4; y++){
        T_old[x][y] = 30.00;
    }//Right Wall

    y = 0;
    for (x = 1; x <= 4; x++){
        T_old[x][y] = 30.00;
    }//Top Wall

    y = 4;
    for (x = 1; x <= 4; x++){
        T_old[x][y] = 30.00;
    }//Bottom Wall

    double alpha = 0.000137;
    double cellSize = 0.01;
    double timeStep = 0.0001;
    double coefficient = alpha*timeStep/(cellSize*cellSize);

    double time = 0.00;

    for (time = 0.00; time <= 300.00; time += timeStep){
        for (x = 0; x <= 4; x++ ){
            for (y = 0; y <= 4; y++){
                T_new[x][y] += coefficient*(T_old[x][y-1]-T_old[x][y]);
            }
        }
        for (x = 0; x <= 4; x++){
            for (y = 0; y <= 4; y++){
                T_new[x][y] = T_old[x][y];
            }
        }
    }
    for (x = 0; x <= 4; x++){
        for (y = 0; y <= 4; y++){
            printf("%d\n",T_new[x][y]);
        }
    }

    return 0;

}

Now, what's wrong with this? Can someone please tell me?

[FleshJeb]

Spoiler

#include<stdio.h>

int main()
{
    double T_old[5][5];
    double T_new[5][5];
    double Tin;
    int x = 0;
    int y = 0;
    printf("\n Enter the Initial Temperature");
    scanf("%d",Tin);
    for (x = 1; x < 4; x++){
        for (y = 1; y < 4; y++){
            T_old[x][y] = Tin;
            T_new[x][y] = Tin;
        }
    }

    //Boundary Conditions
    x = 0;
    for (y = 0; y <= 4; y++){
        T_old[x][y] = 85.00;
    }//Left Wall

    x = 4;
    for (y = 0; y <= 4; y++){
        T_old[x][y] = 30.00;
    }//Right Wall

    y = 0;
    for (x = 1; x <= 4; x++){
        T_old[x][y] = 30.00;
    }//Top Wall


    y = 4;
    for (x = 1; x <= 4; x++){
        T_old[x][y] = 30.00;
    }//Bottom Wall

    double alpha = 0.000137;
    double cellSize = 0.01;
    double timeStep = 0.0001;
    double coefficient = alpha*timeStep/(cellSize*cellSize);

    double time = 0.00;

    for (time = 0.00; time <= 300.00; time += timeStep){
        for (x = 0; x <= 4; x++ ){
            for (y = 0; y <= 4; y++){
                T_new[x][y] += coefficient*(T_old[x][y-1]-T_old[x][y]);
            }
        }
        for (x = 0; x <= 4; x++){
            for (y = 0; y <= 4; y++){
                T_new[x][y] = T_old[x][y];
            }
        }
    }
    for (x = 0; x <= 4; x++){
        for (y = 0; y <= 4; y++){
            printf("%d\n",T_new[x][y]);
        }
    }

    return 0;

}

 

Code view helps ^^^^

This is just from a brief glance, I could be misunderstanding:

• The hot and cold wall resets aren't inside the loop? Shouldn't they be inputting continuously? I don't think the top and bottom walls have to be set more than once, since they're insulators. Just consider 1,0 to 3,0 and 1,4 to 3,4 to be interior nodes.
• It should be T_old = T_new, yours is backwards.
• The [y-1] index goes out of bounds when y=0.

[Cunjo Carl]

4 hours ago, Nivee~ said:

Now, what's wrong with this? Can someone please tell me?

Value of alpha is actually .00001855 m^2/s . Screwed up my conversion- looks like it got you too.   values: 40.6 * 69.7 / (5.91*(10^2)^3 * 25.86)   units W/mK * g/mol  / (g/cm^3 * (cm/m)^3 * J/molK)  -> m^2/s

This line:                 T_new[x][y] += coefficient*(T_old[x][y-1]-T_old[x][y]);   Will try to access an out of bounds variable at the index T_old[0][-1] when y=0 . For starters, use an if statement to prevent.

We also need to do each of its neighbors, not just the one above it, so a line for x+1, x-1, y+1, y-1

Presently walls not maintained at boundary conditions, just initially set to them. Is that intended?

This line:                  T_new[x][y] = T_old[x][y];   Is backwards.

Good luck

Edited June 4, 2019 by Cunjo Carl
Sniped, just gonna leave it.