material change after Ekill-trigger during solution

[julia00.huber]

Hello, 

I'm doing a thermal analysis of timber with gypsum boards during fire. I want to simulate the fall-off of the gypsum board after T_timber has reached a certain temperature. My code is running and calculating correctly, but now I want to change the material property data of timber after the fall-off. I implemented the switch, but now I get a jump in my temperature-time-curve in all depths of the timber after the change even if the temperature at that point is below 100°C (the change in the data only affects the data above 250°C). I have the same temperature points and the material property data is implemented correctly. Can someone help me, why this jump occurs and how I can fix it? It doesn't make sense that the temperature at a depth of 36 mm rises all of the sudden, only because the surface is now exposed. For the surface it makes sense of course.

 This is my solution code: I read the material properties before as material 3 and switch them after the trigger ( the fall off of the first layer can be ignored, currently I only use one layer and kill the second one emediately). I also get the same problem if I read the material data directly during the solution process. 

def SolutionWithGypsumFalloff_AnsysComputation(mapdl, total_time,time_step,t_initial,T_ambient,h_fire,Gypsum_Node,Timber_Node,Tcrit_1st, Tcrit_timber, PD_gypsum_1, PD_gypsum_2,l_timber,save=True):

    #Entering solution process

    mapdl.slashsolu()

    #Setting Initial temperature

    mapdl.allsel()

    mapdl.ic("ALL","TEMP",t_initial)

    #Time stepping definition

    #----------------------------------------------------------------------------------------------

    mapdl.kbc(0) #ramped loads

    mapdl.autots("OFF") #no automatic time steps

    nsteps =int(total_time / time_step)

    mapdl.nsubst(1,1,1) #no further substeps

    #----------------------------------------------------------------------------------------------

    depth_nodes = node_ids_at_depths(mapdl, PD_gypsum_1, PD_gypsum_2)

    gypsum_1_killed = False

    event_substeps = 50         

    event_steps_after = 24      

    event_steps_left = 0  

    gypsum_2_killed = False

   

    times, Gypsum2_surface, Timber_surface,xnode_temp= [], [], [], [] #Collecting data

    temps_by_depth = {y: [] for y in depth_nodes.keys()}

   

    for istep in range(1, nsteps+1):

        current_time = istep * time_step

        mapdl.time(current_time)

        #Setting Loads on Extra-Nodes

        mapdl.run("D,100000,TEMP,%T_bulk%")

        mapdl.d(100001, "TEMP", T_ambient)

        # Set substeps: fine integration around the switch

        if event_steps_left > 0:

            mapdl.nsubst(event_substeps, event_substeps, event_substeps)

            event_steps_left -= 1

        else:

            mapdl.nsubst(1, 1, 1)

        #Solution for time-step

        mapdl.solve()

        #Collect temperature for checking

        t_fire=mapdl.get_value("NODE",100000,"TEMP") #Temperature on Extra_node fire exposed side

        t_cold=mapdl.get_value("NODE",100001,"TEMP") #Temperature on unexposed side correct?

        T_gypsum2=mapdl.get_value("Node", Gypsum_Node, "TEMP") #Temperature gypsum surface

        T_timber= mapdl.get_value("Node",Timber_Node,"TEMP") #Temperature timber surface

        #Temperatures at depths

        for y, nid in depth_nodes.items():

            temps_by_depth[y].append(mapdl.get_value("NODE",nid,"TEMP"))

        #Collecting data

        times.append(current_time)

        Gypsum2_surface.append(T_gypsum2)

        Timber_surface.append(T_timber)

        xnode_temp.append(t_fire)

        #Check: Calculation during Loop correct?

        #-------------------------------------------------------------------------------

        print(f"Time: {current_time}s | Tfire(node100000): {t_fire:.2f}°C |T(100001)={t_cold:.2f}°C | Gypsum2 surface: {T_gypsum2:.2f} °C |Timber surface: {T_timber:.2f}°C")

   

       

        #Check: Fall-off criterion met?

        #-------------------------------------------------------------------------------

        if (not gypsum_1_killed) and (not gypsum_2_killed)and (T_gypsum2 >= Tcrit_1st): #(T_timber>=Tcrit_timber): #

            print("1st gypsum layer failure triggered. Switching exposure to timber.")

            #1) Kill gypsum board

            mapdl.cmsel("S", "gypsum1_board_elem")

            mapdl.ekill("ALL")

            #2) Kill surf element on gypsum surface

            mapdl.cmsel("S", "Fire_surf_gypsum1")

            mapdl.ekill("ALL")

            #3) Activate surface element on timber (now exposed)

            mapdl.cmsel("S", "Fire_surf_gypsum2")

            mapdl.ealive("ALL")

            #COnvection and radiation on timber element

            mapdl.cmsel("S","Fire_surf_gypsum2")

            mapdl.sfe("ALL",1,"CONV",0,h_fire) #Define convection on surface element using film coefficient

            mapdl.allsel()

   

            gypsum_1_killed = True

       

        if (gypsum_1_killed) and (not gypsum_2_killed) and (T_timber >= Tcrit_timber):

            print("2nd gypsum layer failure triggered. Switching exposure to timber.")

            #1) Kill gypsum board

            mapdl.cmsel("S", "gypsum2_board_elem")

            mapdl.ekill("ALL")

            #2) Kill surf element on gypsum surface

            mapdl.cmsel("S", "Fire_surf_gypsum2")

            mapdl.ekill("ALL")

            #Change wood material properties to Kleinhenz values

            #mapdl.prep7()

            mapdl.finish()

           

            # 2) change material number for currently alive timber solid elements

           

            mapdl.prep7()

            mapdl.allsel()

            mapdl.esel("S","MAT","",2)     # timber (old material)

            mapdl.esel("R","TYPE","",2)    # only PLANE292 solids (avoid SURF151)

            mapdl.emodif("ALL","MAT",3)    # switch to Kleinhenz material

            mapdl.allsel()

           

            mapdl.esel("S","TYPE","",2)

            mapdl.esel("R","MAT","",3)

            print("Alive elems with MAT=3:", mapdl.mesh.enum.size)

            mapdl.allsel()

           

            # 3) back to SOLU

            mapdl.finish()

            mapdl.slashsolu()

     

            mapdl.antype("TRANS")

            # 4) Transient-Optionen/Controls zur Sicherheit nochmal explizit setzen

            mapdl.trnopt("FULL")

            mapdl.nropt("FULL")

            mapdl.kbc(0)

            mapdl.autots("OFF")

            #mapdl.nsubst(1, 1, 1)

            #3) Activate surface element on timber (now exposed)

            mapdl.cmsel("S", "Fire_surf_timber")

            mapdl.ealive("ALL")

            #COnvection and radiation on timber element

            mapdl.cmsel("S","Fire_surf_timber")

            mapdl.sfe("ALL",1,"CONV",0,h_fire) #Define convection on surface element using film coefficient

            mapdl.allsel()

   

            gypsum_2_killed = True

            event_steps_left = event_steps_after

    print("Simulation done")        

    if save:

        mapdl.save()

    mapdl.finish()

    return mapdl, times, Timber_surface, xnode_temp, temps_by_depth, depth_nodes

I would be very grateful, if someone could help me!

 

[julia00.huber]

This is the thermal property data before and after switch

 

[ErKo]

 

Hello

 

PyANSYS/PyMAPDL are open source -> supported via the github channels – please post your question there.

General question:

https://github.com/ansys/pymapdl/discussions

 

Of if an issue/bug:

https://github.com/ansys/pymapdl/issues

 

All the best

Erko