Transient conduction in a slab

The steady plane wall pinned the thermal deck’s coefficients and its sign conventions, and it could do nothing else: its exact answer is a straight line, a straight line lies inside the finite-element space, and every mesh reproduces it to round-off. It left one question open. Does refining anything in a thermal solve actually improve the answer, and which thing?

A slab that starts uniformly cold and has both faces stepped to a hot temperature answers that. Its solution is a sum of sines decaying exponentially in time. A sine is not in the space of any polynomial element and an exponential decay is not in the space of any one-step time integrator, so the difference between the run and the closed form here is discretisation error rather than bookkeeping — and there are two independent knobs to turn on it.

Separation of variables on a slab of thickness L, initially at Tᵢ throughout, with both faces held at Tₛ from t = 0. Writing the departure from the surface temperature as a Fourier sine series in x and letting each mode decay at its own rate:

    (T(x,t) - T_s) / (T_i - T_s)
        = sum over odd n of  (4/(n pi)) sin(n pi x/L) exp(-(n pi/L)^2 alpha t)

    alpha = lambda / (rho c) = 1.19427e-5 m^2/s

The even modes vanish by symmetry, which is the statement that both faces are doing the same thing. The measurement is taken at the midplane, x = L/2, where the sine is ±1 and the series is at its slowest-converging — and where, by symmetry, the temperature is the minimum anywhere in the body.

The comparison time is chosen so that the answer is a real number rather than either endpoint. At t = 104.7 s the Fourier number at the half-thickness is αt/(L/2)² = 0.5002, and the midplane has travelled 62.9 % of the way from its initial 300 K to the imposed 800 K.

The textbook first transient conduction case is a semi-infinite solid under a step surface temperature, whose profile is an error function. It cannot be read from what this product writes, and the reason is worth stating because it shaped the case.

The thermal deck’s summary is the minimum and the maximum temperature in the body. For a semi-infinite solid the maximum is the imposed surface temperature — an input — and the minimum is the far end, which is only at a measurable temperature once the thermal wave has reached it. By that point the body is no longer semi-infinite and the closed form no longer applies. The two conditions exclude each other.

Heating a slab from both faces has no such problem. The coldest point is the midplane by symmetry, so the minimum the product already reports is the quantity with a closed form behind it. Nothing on this page opens a result file the product did not offer.

Read that table down the rows rather than across, because the two rows say opposite things.

• A 113-fold refinement of the mesh changes nothing. 19 tetrahedra to 2 136, at the same time step, moves the answer from −0.910 % to −0.974 % — very slightly worse. That is the coarse mesh’s spatial error happening to offset part of the temporal error, and it is a coincidence rather than a reason to use a coarse mesh.
• A five-fold refinement of the time step cuts the error 5.06-fold. −0.974 % to −0.193 %, for a step ratio of exactly 5. That is an observed order of 1.007 — first order in time, to two figures.

The practical consequence is the opposite of what a CFD habit suggests: refining the mesh to improve a transient thermal answer does nothing here. The time step is the control that matters, and halving it halves the error rather than quartering it.

First order in time invites an obvious diagnosis and an obvious fix — that the solver is integrating with backward Euler, and that asking for Crank–Nicolson would make it second order for free. Both halves of that are wrong, and it took four more solves to find out. Same mesh, same 5.235 s step, the integration parameter θ set explicitly:

Three things follow, and none of them is the obvious diagnosis.

• The keyword unquestionably bites. θ = 0 blows the solve up, so the parameter is reaching the integrator rather than being ignored.
• The default is not backward Euler. It lands between θ = 0.5 and θ = 1.0, which is consistent with code_aster’s documented default of 0.57 — slightly over-implicit on purpose, to damp the oscillations a stepped boundary condition provokes.
• Crank–Nicolson does not recover second order here, and no keyword will. At θ = 0.5 the observed order is 1.009 — the error still falls 5.07-fold for a five-fold smaller step. The problem has a discontinuous initial condition: a body at a uniform temperature with both faces suddenly stepped has an infinitely sharp gradient at t = 0, and a one-step θ-scheme loses its formal order through non-smooth data. Damping the first few steps is the standard remedy for that and it is not a keyword.

So the first-order behaviour on this page is a property of this problem rather than a defect in the product, and the improvement on offer from switching θ is 2.7 % of the error, in exchange for the oscillations the 0.57 default exists to suppress on exactly this class of boundary condition. The product is left alone, and this page reports the order rather than treating it as a bug.

• One temporal refinement pair, not a sweep. The order comes from a single 5.235 s → 1.047 s comparison. It is corroborated by a second pair at θ = 0.5 giving 1.009, but that is still two points each.
• The 0.57 default is inferred, not read out of the solver. What is measured is that the default answer sits between the θ = 0.5 and θ = 1.0 answers. That it is exactly 0.57 is code_aster’s documented value, quoted here as the consistent explanation rather than as something this case measured.
• The attribution of the lost order to the stepped initial condition is an argument, not a measurement. It rests on the θ = 0.5 run still measuring 1.009. A run with a smooth initial condition would settle it, and has not been done.
• One material, constant properties, no source. λ, ρ and c are constant, there is no volumetric heating and there is no interface. All three are what make the series exact.
• Boundary conditions are Dirichlet only. The flux and film conditions are exercised by the steady plane wall, not here.