Problems writing a Heat Transfer Equation

56 views
Skip to first unread message

Maximilian Nast

unread,
Dec 4, 2023, 3:32:08 PM12/4/23
to pysph-users
Hello fellow friends of PySPH,

i've been trying to implement Monaghan Style Heat Transfer in PySPH but ran into a few Problems. I have probably misunderstood something, and would be very thankful for any hints.

I've based my Code on this Equation:
Screenshot 2023-12-04 194431.png

It ends up looking something like this:

 ```Python 
class HeatConduction(Equation):
    def __init__(self, dest, sources):
        super(HeatConduction, self).__init__(dest, sources)

    def loop(self, d_idx, s_idx, d_T, s_T, DWIJ, d_dT, s_m, d_rho, s_rho, R2IJ, XIJ):
        # i ^= d; j ^= s
        dot_product = DWIJ[0] * XIJ[0] + DWIJ[1] * XIJ[1] + DWIJ[2] * XIJ[2]

        if d_idx == s_idx:
            d_dT[d_idx] += 0
        else:
            d_dT[d_idx] += ((s_m[s_idx] / (d_rho[d_idx] * s_rho[s_idx]))
                            * ((4 * 4500 * 4500) / (4500 * 4500))
                            * (d_T[d_idx] - s_T[s_idx])
                            * (dot_product / R2IJ))

    def converged(self):
        return 1  # assume convergence for now
```

With the acccording Integrator looking something like:

 ```Python 
class HeatStep(IntegratorStep):
    def initialize(self, d_idx, d_T):
        pass

    def stage1(self, d_idx, d_T, d_dT, dt):
        pass
   
    def stage2(self, d_idx, d_T, d_dT, dt):
        d_T[d_idx] += dt * d_dT[d_idx]
```

I've included the full .py file as an attachment as well. It will run and visualize the simulation.

For the first few steps, it looks like what i would expect:

Screenshot 2023-12-04 211350.png

As you can see in the Image above, the rate of change of the Temperature "dT" after 10 timesteps is positive on the initially cold (upper) side and negative on the warm side, so the temperature should equalize (and it does for the first 8 steps or so).

However, after ~10 steps, the simulation becomes very unstable and the temperature and its rate of change start to oscillate back and forth between the two sides like this:

Screenshot 2023-12-04 211834.png

This Problem of waves of Heat travelling up and down can be seen more clearly in this 2d example:

Screenshot 2023-12-04 212625.png

I believe i must have made a mistake implementing the stepper or equation, but am neither smart nor experienced enough to figure out where.

Any hint in the right direction is highly appreciated!

Thank you in advance
Max
sim.py

A Dinesh

unread,
Dec 4, 2023, 3:50:26 PM12/4/23
to Maximilian Nast, pysph-users
Hi Max,

It is hard to debug other codes. Also, one way to verify our codes is to simulate standard benchmarks, such as Taylor-Green vortex or dam break  for fluids, cantilever beam deformation for elastic solids etc.

Similar to these, for a heat transfer (which is similar to diffusion equation) we have other standard examples. Here, I am attaching my implementation of diffusion equation applied for concentration difference (your case it is temperature). Please save the following codes and run the `mayoral_2016` file for the result. The resulting figure will be saved in the output folder. Here, I am verifying one of the SPH discretizations of the diffusion equation.

Feel free to follow up.

Cheers,

Dinesh.

--
You received this message because you are subscribed to the Google Groups "pysph-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/1c18122d-1d4f-4043-b199-c9b640137d27n%40googlegroups.com.
fluids.py
mayoral_2016_4_1_diffusion_in_aqueous_solution_1d.py

Maximilian Nast

unread,
Dec 4, 2023, 6:22:14 PM12/4/23
to pysph-users
Hi Dinesh,

Thank you so much for your input! I went through your code line by line and for the relevant parts, it's already very similar to mine. For future reference: I was just missing this single function in my Equation ...

def initialize(self, d_idx, d_aconcentration):
    d_aconcentration[d_idx] = 0.0

or in my case

def initialize(self, d_dT, d_idx):
    d_dT[d_idx] = 0

Implementing this was genuinely all it took! Now i get this result, which is what i was looking for:

ezgif.com-crop.gif

I will now, of course, check this result against an analytical solution as you suggested, but it looks good so far.

Thank you again!

Best regards
Max

A Dinesh

unread,
Dec 5, 2023, 2:00:29 AM12/5/23
to Maximilian Nast, pysph-users
Great, congrats.

Sent from Outlook for iOS

From: 'Maximilian Nast' via pysph-users <pysph...@googlegroups.com>
Sent: Monday, December 4, 2023 11:22:13 PM
To: pysph-users <pysph...@googlegroups.com>
Subject: Re: Problems writing a Heat Transfer Equation
 
Reply all
Reply to author
Forward
0 new messages