Блог

Я не знаю, что я делаю не так. Я должен решить систему нелинейных уравнений, которые включают радиационную теплопередачу, но каждый раз, когда я пытаюсь ее реализовать, алгоритм не сходится, при ier = 5, что означает «Итерация не продвигается хорошо, как измеряется улучшением по сравнению с последними десятью итерациями». Я довольно новичок в Python, но я не знаю, будет ли все это иметь смысл без всего контекста. Смысл моего кода состоит в том, чтобы использовать своего рода метод конечных элементов для вычисления поля температуры в печи для остекловывания. Это 1D-модель с большим количеством гипотез, и меня не интересует решение поля скоростей, создаваемого естественной конвекцией. Печь состоит из трех элементов : перерабатываемых отходов, огнеупорного песка и металлической оболочки, поэтому я использовал классы для каждого элемента, чтобы каждый раз сохранять физические свойства каждого элемента. Каждый объект каждого слоя печи затем хранится в объекте печи. Код выглядит следующим образом :

«‘ деф Билан(self, dt):

     dim_x = self.Nbcomposants                                   
    dim_y = self.Ncouches
    SFour = self.width**2                                      
    dz = self.height/dim_y
    S_MR = self.width*dz                       
    Pc = self.puissance
    dim_T = dim_x*dim_y                                       
    A = np.zeros((dim_T, dim_T))                                #  Ax = B
    B = np.zeros(dim_T)
    T_old = np.zeros(dim_T)                                     

    for i in range(self.NombredeCouchesLiquides, dim_y):
        if self.Couches[i].composants[2].Temperature > Tfusion:
            self.Couches[i].composants[2].IsLiquid = True
            self.Couches[i].composants[2].IsSolid = False
            self.NombredeCouchesLiquides  = 1
            self.Couches[i].composants[2].proportionLiquide = 1
                        

    for x in range(0, dim_x):                             
        for y in range(0, dim_y):                               
            i = y   dim_y*x
            T_old[i] = self.Couches[y].composants[x].Temperature
            

    for x in range(0, dim_x):
        for y in range(0, dim_y):
            Layer_i = self.Couches[y]
            Main = Layer_i.composants[x]

            i = y   dim_y*x

            if x == 0:                                          # metal liner

                A[i, i] = (Main.rho*Main.Cp*Main.dl/dt   2/(Layer_i.composants[x 1].dl/Layer_i.composants[x 1].lambdath   Main.dl/Main.lambdath)   Main.hConvectionAir)*S_MR
                A[i, i dim_y] = -2*S_MR/(Layer_i.composants[x 1].dl/Layer_i.composants[x 1].lambdath   Main.dl/Main.lambdath)
                B[i] = (Main.rho*Main.Cp*Main.dl*T_old[i]/dt   Main.hConvectionAir*Tamb   sigma*(Tamb   273.15)**4)*S_MR

                if y != 0:
                    A[i, i]  = self.width*Main.dl/(dz*(1/Main.lambdath   1/self.Couches[y-1].composants[x].lambdath))
                    A[i, i-1] = -self.width*Main.dl/(dz*(1/Main.lambdath   1/self.Couches[y-1].composants[x].lambdath))

                if y != dim_y - 1:
                    A[i, i]  = self.width*Main.dl/(dz*(1/Main.lambdath   1/self.Couches[y 1].composants[x].lambdath))
                    A[i, i 1] = -self.width*Main.dl/(dz*(1/Main.lambdath   1/self.Couches[y 1].composants[x].lambdath))


            elif x == 1:                                          # refractory sand

                A[i, i] = (Main.rho*Main.Cp*Main.dl/dt   2/(Main.dl/Main.lambdath)   2/(Layer_i.composants[x-1].dl/Layer_i.composants[x-1].lambdath   Main.dl/Main.lambdath))*S_MR
                A[i, i-dim_y] = -2*S_MR/(Layer_i.composants[x-1].dl/Layer_i.composants[x-1].lambdath   Main.dl/Main.lambdath)
                A[i, i dim_y] = -2*S_MR/(Main.dl/Main.lambdath)
                B[i] = Main.rho*Main.Cp*S_MR*Main.dl*T_old[i]/dt

                if y != 0:
                    A[i, i]  = self.width*Main.dl/(dz*(1/Main.lambdath   1/self.Couches[y-1].composants[x].lambdath))
                    A[i, i-1] = -self.width*Main.dl/(dz*(1/Main.lambdath   1/self.Couches[y-1].composants[x].lambdath))

                if y != dim_y - 1:
                    A[i, i]  = self.width*Main.dl/(dz*(1/Main.lambdath   1/self.Couches[y 1].composants[x].lambdath))
                    A[i, i 1] = -self.width*Main.dl/(dz*(1/Main.lambdath   1/self.Couches[y 1].composants[x].lambdath))

            elif x == 2:                                          # Melt

                A[i, i] = Main.rho*Main.Cp*SFour*dz/dt   2*S_MR/(Layer_i.composants[x-1].dl/Layer_i.composants[x-1].lambdath)
                A[i, i-dim_y] = -2*S_MR/(Layer_i.composants[x-1].dl/Layer_i.composants[x-1].lambdath)
                B[i] = Main.rho*Main.Cp*SFour*dz*T_old[i]/dt

                if y != 0:
                    A[i, i]  = SFour/(dz*(1/Main.lambdath   1/self.Couches[y-1].composants[x].lambdath))
                    A[i, i-1] = -SFour/(dz*(1/Main.lambdath   1/self.Couches[y-1].composants[x].lambdath))
                else:
                    B[i]  = Pc

                if y != dim_y - 1:
                    A[i, i]  = SFour/(dz*(1/Main.lambdath   1/self.Couches[y 1].composants[x].lambdath))
                    A[i, i 1] = -SFour/(dz*(1/Main.lambdath   1/self.Couches[y 1].composants[x].lambdath))
                else:
                    A[i, i]  = self.Couches[dim_y][1].hconvection*SFour
                    B[i]  = SFour*(self.Couches[dim_y][1].hconvection*self.Couches[dim_y][1].Temperature) #
                    #  sigma*(FMeltToHood*(self.Couches[dim_y 1].Temperature   273.15)**4   2*FMeltToSideWall*(self.Couches[dim_y][0].Temperature   273.15)**4))

            else:
                print('Le nombre de composants dans une couche ne correspond pas')


    def CalculT(T, arg1, arg2):                   # solve Ax = B, adding radiative elements
        X = np.ravel(T)
        f = np.dot(arg1,X)
        f -= np.ravel(arg2)
        F = np.ravel(f)
        for i in range(dim_y):
            F[i]  = sigma*S_MR * (X[i]   273.15)**4           # radiation only for metal liner (0 to 99)

        F[dim_T-1]  = SFour*sigma*(FMeltToHood*(X[dim_T-1]   273.15)**4   2*FMeltToSideWall*(X[dim_T-1]   273.15)**4)
        return F


    results = fsolve(CalculT, T_old, args = (A,B), full_output = True)
    T_new = results[0]
    exitflag = results[2]

    if exitflag != 1:
        print("L'algorithme ne converge pas")
 

»’

Я действительно не знаю, понятно ли это без полного набора гипотез и уравнений, но если кто-то каким-то образом поймет, почему это не работает, я с радостью приму любой совет. Код отлично работает, если я прокомментирую детали, связанные с радиационной передачей тепла.