6. Coupled problems

6.1. One-Way Thermomechanical Coupling

Location of example: $B2VERIFICATION/coupled/oneway`

The example demonstrates how to solve a one-way coupled thermomechanical problem. One-way coupled means that solutions from a stand-alone transient heat transfer analysis are used to solve the thermally induced deformation problem.

A (non-linear) transient heat transfer analysis is obtained by a single transient heat-transfer analysis run, solving the heat transfer equation

\[\rho c_{p}(T) \frac{\partial T}{\partial t} - k_{x}(T) \frac{\partial^2 T}{\partial^2 x} - k_{y}(T) \frac{\partial^2 T}{\partial^2 y} - k_{z}(T) \frac{\partial^2 T}{\partial^2 z} = 0\]

with \(c_p\) and \(k_i\) constant. Thereafter, transient heat analysis solutions for all transient heat analysis time steps are taken one by one and a linear static temperature-induced deformation analysis for each time step is performed.

The structure is a steel plate 1 [m] by 0.1 [m] by 0.01 [m] and is modeled with HE8 solids. Material data are defined in the respective B2000++ input files heat.mdl and deform.mdl.

The plate is clamped on the face x=0, and thus produces deformations but no stresses when subject to pure thermally induced strains. Both the heat transfer and the deformation analysis meshes FE meshes are identical!

_images/examples.coupled.oneway.mesh.svg

FE Mesh (identical for heat transfer and deformation analysis).

The heat-transfer FE model boundary conditions are selected as follows:

  • A constant amount of heat per volume is applied a t=0 [sec] and kept constant during the time interval. This is the condition nbc 1.

  • The initial temperature at t=0 for all nodes is set to 0.

  • The reference temperature is 0.

  • The time interval is from 0 to 100 [sec], in constant time steps of 5 [sec].

These conditions will produce a constant temperature over the whole model for any given time step, and the temperature will rise linearly with time.

The deformation analysis FE model boundary conditions for each time step are selected as follows:

  • Constraints: \(U_x=U_y=U_z=0\) for all nodes at x=0 (clamped face). All other nodes are free to move.

  • Initial strains produced by the temperature field for each time step. The temperature field for each time step is copied from the heat analysis model database.

As mentioned before, the constant temperature for any given time step produces deformations, i.e a linear extension of the plate in x, but no stresses when subject to pure thermally induced strains.

_images/examples.coupled.oneway.deform-bc.svg

Deformation analysis: Boundary conditions.

An analytical (although one-dimensional) simple solution of the plate end deformation in the x-direction at x=1m is \(U_x = \Delta T \alpha\), see coupled.py.

The coupled analysis requires a Python script coupled.py for controlling the analysis. The script contains explanatory text:

"""One-way thermo-mechanical coupling test case execution control
   script
"""

import os
import memcom
import simples as si
import matplotlib.pyplot as plt

def copy_dof_field(db1, name, db2, cycle2=None):
    """Copy DOF field(s) name="<gname>.*.<cycle>.0.<case>" from db1 to db2
    as "<gname>.*.<cycle>.0.<case>". Copies descriptors, too! If
    cycle2=None, cycle=cycle1.

    """
    sets = db1.keys(name)
    for s in sets:
        a = s.split('.')
        if cycle2 is None:
            cy = a[2]
        else:
            cy = f'cycle'
        n1 = f'{a[0]}.{a[1]}.{a[2]}.{a[3]}.{a[4]}'
        n2 = f'{a[0]}.{a[1]}.{cy}.{a[3]}.{a[4]}'
        db2[n2] = db1[n1][:]
        db2[n2].desc = db1[n1].desc

m2mm = 1000.

timestep = 100. # from heat.mdl !

# Step1: Run transient heat transfer analysis
# -------------------------------------------

print(f'Run transient heat analysis')
os.system(f'b2000++ -define TSTEP={timestep} heat.mdl')

# Step 2: Create deformation analysis mesh database
# -------------------------------------------------

print(f'Create deformation analysis model database')
os.system('b2ip++ deform.mdl')

# Step 3: Open heat model in Simples
# ----------------------------------

heat_model = si.Model('heat.b2m')
db1 = heat_model.get_db() # MemCom database handle
heat_steps = heat_model.get_steps(case=1) # Get heat solution times

# Step 4: Solve static deformations
# ---------------------------------

for itime, time in enumerate(heat_steps):
    # Loop over hat transfer solution steps
    cycle = itime + 1
    print(f'Run deformation analysis for step {cycle}, time={time}')

    # Get temperature field for current time. Add T of point p2 to
    # graph.
    temp = si.NodeField(heat_model, 'TEMP', case=1, cycle=cycle)

    # Open the MemCom deformation analysis database in write mode - it
    # will be modified: The temperature solution DOF field will be
    # added and the case will be modified.
    db2 = memcom.db('deform.b2m/archives.mc', 'w')
    name = f'TEMP.1.{cycle}.0.1'

    # Copy temperature solution DOF field from heat to deformation
    # database.
    copy_dof_field(db1, name, db2)

    # Add equivalent to MDL CASE option "temperature_field <case>
    # <cycle>" to database. This is a hack until a bug in Simples is
    # fixed...
    copy_dof_field(db1, name, db2)
    db2['CASE.1'][0]['ID'] = name
    db2.close()

    # Make static deformation analysis for step. Note that solutions
    # are saved for static analysis (cycle=0).

    os.system('b2000++ deform.b2m')

    # Get displacement solution for step.

    # Rename linear static analysis solutions to create 'transient'
    # solutions at step <cycle> on defomation db.
    db2 = memcom.db('deform.b2m/archives.mc', 'w')
    db2[f'DISP.1.0.0.1'].name = f'DISP.1.{cycle}.0.1'
    db2[f'FORC.1.0.0.1'].name = f'FORC.1.{cycle}.0.1'
    db2[f'RCFO.1.0.0.1'].name = f'RCFO.1.{cycle}.0.1'
    db2[f'STRAIN.1.0.0.1'].name = f'STRAIN.1.{cycle}.0.1'
    db2[f'STRESS.1.0.0.1'].name = f'STRESS.1.{cycle}.0.1'

    ds = f'SOLUTION-STEP.0.{cycle}.0.1'
    db2[ds] = db1[ds][:]

    db2.close()

# Write remaining SOLUTIION sets to deformation database. Reopen.
db2 = memcom.db('deform.b2m/archives.mc', 'w')
ds = f'SOLUTION.0.0.0.1'
d = dict(db1[ds][:]) # Copy descriptor of heat analysis and replace items
d['DOF_SOL'] = 'DISP'
d['DOF_DOT_SOL'] = 'DISPD'
d['DOF_DOTDOT_SOL'] = 'DISPDD'
d['NBC_SOL'] = 'FORC'
d['RESIDUUM_SOL'] = 'RCFO'
d['SP_SOL'] = 'COOR_IP. DISP_IP,STRAIN,STRESS'
del db1[ds]
db2[ds] = d

# Rename the SP fields that are independent of cycle because of linear
# deformation analysis.
db2[f'COOR_IP.1.0.0.1'].name = f'COOR_IP.1.1.0.1'
db2[f'DISP_IP.1.0.0.1'].name = f'DISP_IP.1.1.0.1'
db2[f'ELEM_IP.1.0.0.1'].name = f'ELEM_IP.1.1.0.1'
db2.close()

# Step 5: Postprocessing
# ----------------------

t_axis = [] # Graph time
T_axis = [] # Graph T(emperature) at point p2
Ux_axis = [] # Graph Ux (displacement in global x) at point p2
Ux_analytical_axis = []

#heat = si.Model('deform.b2m') # Is used if heat gradients are accessed
deform = si.Model('deform.b2m')
steps = deform.get_steps(case=1) # Get solution times
node = deform.get_epatch_point_node(1, 2) #  Get node id of mesh patch point p2

# Extract values
for itime, time in enumerate(steps):
    cycle = itime + 1
    temp = si.NodeField(deform, 'TEMP', case=1, cycle=cycle)
    T_axis.append(temp[node][0])
    t_axis.append(time)
    Ux_axis.append(si.NodeField(deform, 'DISP', case=1, cycle=cycle)[node][0]*m2mm)
    Ux_analytical_axis.append(temp[node][0]*11.7e-6*m2mm)

# Create csv file p2.csv with values at node P2
f = open("p2.csv", "w")
f.write(f'{"time":>6s}, {"Temp":>5s}, {"Ux":>6s}, {"  Ux_analytical"}\n')
for t, T, Ux, Ux_analytical in zip(t_axis, T_axis, Ux_axis, Ux_analytical_axis):
    f.write(f'{t:>6.2f}, {T:>5.2f}, {Ux:>.5f}, {Ux_analytical:>.5f}\n')
f.close()
print('Created CSV file "p2.csv""')

# Plot values

fig, axs = plt.subplots(1, 2, figsize=(9, 3))

axs[0].plot(t_axis, T_axis)
axs[0].set_xlim(0, timestep)
axs[0].set_xlabel('Time [s]')
axs[0].set_ylabel('Temperature [degC]')
axs[0].grid(True)

axs[1].plot(t_axis, Ux_axis, label='computed')
axs[1].plot(t_axis, Ux_analytical_axis, marker = 'o', linewidth=0,
            label='analytical')
axs[1].set_xlim(0, timestep)
axs[1].set_xlabel('Time [s]')
axs[1].set_ylabel('Ux (point P2) [mm]')
axs[1].grid(True)
axs[1].legend()

fig.suptitle('One-way Coupled Thermoelastic Simulation')
fig.tight_layout()
axs[1].legend()
fig.savefig('output.svg')
print('Created B2000++ model database "heat.b2m"')
print('Created B2000++ model database "deform.b2m"')
print('Created graph "output.svg"')
python3 coupled.py
Run transient heat analysis
Create deformation analysis model database
Run deformation analysis for step 1, time=5.0
Run deformation analysis for step 2, time=10.0
Run deformation analysis for step 3, time=15.0
Run deformation analysis for step 4, time=20.0
Run deformation analysis for step 5, time=25.0
Run deformation analysis for step 6, time=30.0
Run deformation analysis for step 7, time=35.0
Run deformation analysis for step 8, time=40.0
Run deformation analysis for step 9, time=45.0
Run deformation analysis for step 10, time=50.0
Run deformation analysis for step 11, time=55.0
Run deformation analysis for step 12, time=60.0
Run deformation analysis for step 13, time=65.0
Run deformation analysis for step 14, time=70.0
Run deformation analysis for step 15, time=75.0
Run deformation analysis for step 16, time=80.0
Run deformation analysis for step 17, time=85.0
Run deformation analysis for step 18, time=90.0
Run deformation analysis for step 19, time=95.0
Run deformation analysis for step 20, time=100.0

Created CSV file "p2.csv
Created B2000++ model database "heat.b2m"
Created B2000++ model database "deform.b2m"
Created graph "output.svg"

Executing the script coupled.py creates both the heat.b2m heat transfer analysis and the deform.b2m deformation analysis B2000++ model databases. Both databases contain solutions for the synchronized time steps:

  • heat.b2m contains, for all time steps, the DOF solution sets temperatures TEMP, time derived temperatures TEMPD, heat HEAT, and residual heat RHEAT, and the sampling point (element gradient) solutions heat capacity HEAT_CAPACITY, heat flow temperatures HEAT_FLOW,

  • deform.b2m contains, for all time steps, the DOF solution sets displacements DISP, forces FORC, reaction forces RCFO, and the sampling point (element gradient) solutions stresses STRESS, and strains STRAIN. Note that all stresses are equal to 0 (see comment above).

A “Mathplotlib* graph output.svg with time-temperature and time-displacement curves for the selected FE mesh node P2 is generated by coupled.py:

_images/examples.coupled.oneway.graph.svg

Coupled analysis solution: Time-temperature and time-displacement curves for the selected FE mesh node P2.