Neqsim Fem Coupling

FEM Coupling

NeqSim solves the fluid. It reports a bulk temperature, a pressure and a set of transport properties, and a one-dimensional model turns those into a single U-value per pipe section. That works until the geometry stops being one-dimensional: a section of insulation floods, a clamp bridges the coating, a line is buried under a sloping seabed, a nozzle interrupts a shell. The heat then spreads in a direction the one-dimensional model does not have, and the U-value it was given is wrong by an amount nobody can estimate by hand.

This skill closes that gap. It takes whatever is known about the component (P&ID, STID tag register, datasheets, insulation specification, inspection report), takes the fluid side from a NeqSim flash, builds and solves the finite-element model of the solid, and converts the field back into the number the one-dimensional model needed - an effective U-value, a multiplier on it, a hot-spot factor, a no-touch time, a wall stress. When a thermal or stress finite-element study already exists, it qualifies that study instead of running a new one.

It is the solid-side companion to neqsim-cfd-coupling: same design-basis discipline, same quality gate before any number is quoted, same handoff back into a one-dimensional model. CFD resolves the fluid; FEM resolves the solid.

Every layer works on its own. NeqSim, Gmsh, scikit-fem and FEniCSx are all optional: without NeqSim the fluid state is stated directly, and without a mesher or a solver the geometry and the case are still written, with the commands needed to run them elsewhere.

When to Use

• A one-dimensional U-value is being assumed where the geometry is not one-dimensional: a local insulation defect or flooded section, a support or clamp, a buried line, a valve or flange with no insulation.
• A cooldown or no-touch time is needed and the wall's thermal inertia matters, not just the fluid's.
• A metal temperature is needed for a thermal-stress, MDMT or material-selection question, and the process-to-ambient difference is not the metal gradient.
• Species diffusion through a porous medium is being modelled and the molecular diffusivity has to come from the actual mixture.
• A thermal or thermo-mechanical finite-element report is attached to an equipment tag and its numbers are about to be used.
• A new finite-element run is being specified and needs film coefficients, a penetration depth, an element size and a time step.

When Not to Use

• As a substitute for a qualified stress or thermal analyst on a design decision. This is a screening and coupling toolkit, not a code-compliance calculation.
• For fatigue life, fracture mechanics, creep, plasticity, buckling or contact. The stress layer produces elastic thermal and pressure stresses with a category attached; it does not perform a code assessment.
• For the flow field. Velocity, wall shear and flow maldistribution belong to neqsim-cfd-coupling; this skill consumes a film coefficient, it does not resolve the boundary layer that produces it.
• To extract a number from a study that fails the quality gate.

Inputs

Outputs

Engineering Method

Design basis by precedence. Sources are ranked measurement > mechanical_datasheet > material_certificate > insulation_specification > process_datasheet > vendor > stid > inspection_report > plant_data > pid > estimate > assumption. A lower-ranked source never overwrites a higher-ranked one; if it disagrees by more than 2 % it raises a conflict, and any disagreement about a non-numeric field such as a material grade is always a conflict. The model kind determines which fields are required, so an under-specified model is reported rather than meshed.

Fluid to boundary condition. A finite-element model consumes a film coefficient, not a fluid. Properties come from the NeqSim phase after initProperties(), because thermal conductivity and viscosity are otherwise zero. The default correlation is Gnielinski with the Petukhov friction factor,

$$f = (0.79 \ln \mathrm{Re} - 1.64)^{-2}, \qquad \mathrm{Nu} = \frac{(f/8)(\mathrm{Re} - 1000),\mathrm{Pr}} {1 + 12.7\sqrt{f/8},\left(\mathrm{Pr}^{2/3} - 1\right)}$$

valid over a far wider Prandtl range than Dittus-Boelter, which matters because a dense-phase gas and a glycol are on opposite sides of that range. Below Re = 2300 the fully developed laminar value Nu = 3.66 is used and said so.

Do I even need a mesh? The Biot number over the whole conduction path, $\mathrm{Bi} = h L / k$ with $h$ the smaller of the two surface coefficients, decides it. Below 0.1 the solid is nearly isothermal and a lumped model answers the question. The skill says so rather than meshing anyway.

Discretisation from the physics. For a transient, the thermal penetration depth $\delta = \sqrt{\alpha t}$ is how far the front has moved after the time of interest; elements coarser than $\delta/4$ cannot represent it. The time step targets a mesh Fourier number $\alpha \Delta t / \Delta x^2$ of about 0.5 - an implicit scheme is stable at any step, so the step is chosen to resolve the front, not to keep the solve from diverging.

One dimension, verified. Layered walls are assembled as linear P1 elements on a one-dimensional mesh with every layer interface on a node, in cylindrical or planar coordinates, with Robin conditions on both surfaces. Temperature-dependent conductivity is handled by Picard sweeps on the layer mean temperature. The steady result is checked against the closed-form composite resistance

$$R' = \frac{1}{h_i 2\pi r_i}

• \sum_j \frac{\ln(r_{j+1}/r_j)}{2\pi k_j}
• \frac{1}{h_o 2\pi r_o}$$

and the deviation is reported. That check is the reason to start in one dimension: a two-dimensional field cannot be verified this way, so if the one-dimensional answer is not understood first, nothing downstream is.

Cooldown. Holding the bore fluid at a fixed bulk temperature models a thermal shock, not a cooldown - a shut-in line has no source keeping it warm. Supplying inner_fluid_capacity = $\rho c_p A_\text{bore}$ adds the inventory as one lumped node coupled to the wall through the inner film, and the system stays tridiagonal. NeqSim supplies the density and heat capacity of the shut-in fluid and the hydrate or wax temperature the cooldown is measured against.

Mesh. Layered geometry wants a structured tensor grid, one block per (axial segment, through-thickness layer), because an unstructured mesh dropped over a 12 mm wall next to 50 mm of insulation puts one element across the wall - and the wall gradient is what sets the thermal stress. Interfaces sit on element boundaries, physical groups carry the material assignment and the faces, and the group ids are written explicitly so a name-based reader (meshio, scikit-fem) and a tag-based reader (dolfinx) see the same partition. A local defect is represented by overriding a segment's material - water-flooded insulation conducts like seawater, which is what actually happens.

Backend screening. Python has no single dominant finite-element package, so the choice is stated rather than defaulted:

Case as data. inputs.json holds the mesh reference, the per-material properties and the boundary conditions; case.py is a fixed script that reads it. The same inputs.json drives the scikit-fem and the FEniCSx script, so a model is promoted from the light backend to the heavy one without being rebuilt.

Three dimensions, and when not to use them. revolve_deg sweeps the r-z section about the axis into a pipe or vessel wall; extrude_m sweeps a plane section or a block out of plane. Gmsh returns a sweep in a fixed order - end surface, volume, then one lateral face per curve of the loop - and the loop is written as (inner, east, outer, west), so the volumes and the six faces are tagged without any geometric search and the same physical names carry through. A partial revolve adds symmetry_start and symmetry_end, which are cut planes and take an adiabatic condition. ConductionProblem.from_mesh_spec switches the axisymmetric weighting off for a swept mesh, because a revolved model already contains the circumference and weighting it again counts it twice.

Do not revolve by default. An axisymmetric problem on a revolved mesh costs an order of magnitude more and returns the same answer as the r-z section. The reasons to go to three dimensions are a genuinely circumferential feature - a support, a nozzle, a weld that does not run all the way round, a partial-thickness defect - and presentation.

Rendering. render_field reads whatever field file the case recorded (field.vtu from scikit-fem, field.xdmf from FEniCSx) and produces off-screen PNGs. A two-dimensional mesh gets one surface view; a three-dimensional mesh gets the surface, a cut plane and a clipped view, because the outside of a solid says almost nothing about the gradient inside it. The camera defaults to the plane of the two largest extents - an isometric view of a 3.5 mm wall on a 1 m tube wastes the frame - and scale stretches the thin direction so the layers are legible. Both are presentation devices and are labelled as such in the returned message. A picture is not evidence: render after the quality gate, not instead of it.

Quality gate. Four things decide whether a finite-element number means anything.

For a transient, the mesh Fourier number and the number of time steps are checked too: an implicit scheme stays stable while smearing the front, which flatters a cooldown time and understates a thermal shock.

Handoff. The field is reduced to what a one-dimensional model consumes. The overall U-value follows from the boundary heat flow and the reference area, and the U-multiplier - the ratio of the finite-element heat flow to the one-dimensional heat flow over the same length - is the number a NeqSim pipeline or cooldown model carries. It is what transports a local defect, a support short-circuit or a non-radial soil path into a model that cannot represent them.

Stress category. Pressure stress is primary: it does not relax, and exceeding the allowable means the wall fails. Restrained-expansion and through-wall-gradient stresses are secondary and self-limiting, so they are assessed against a range allowable of order $3S$ and against fatigue, not against the primary membrane allowable. Thermal stresses used here are

$$\sigma_\text{axial} = E\alpha,\Delta T, \qquad \sigma_\text{biaxial} = \frac{E\alpha,\Delta T}{1-\nu}, \qquad \sigma_\text{gradient} = \frac{E\alpha,\Delta T}{2(1-\nu)}$$

with the Lame thick-wall pressure solution superposed. Comparing a thermal stress with a primary allowable is the most common way to condemn a wall that is acceptable - or to pass one that will crack in cyclic service.

Python Usage Pattern

One dimension first - and verified

from fem_coupling import (
    ConductionLayer, RadialConductionModel, build_design_basis,
    derive_thermal_conditions, film_coefficient, fluid_state_from_neqsim, material,
)

basis = build_design_basis(
    tag="20-P-001",
    model_kind="insulated_pipe",
    sources=[
        {"source": "stid", "reference": "STID line 20-P-001",
         "values": {"inside_diameter_m": 0.254, "wall_thickness_m": 0.0127,
                    "wall_material": "carbon-steel"}},
        {"source": "insulation_specification", "reference": "SPEC-INS-004 rev 1",
         "values": {"insulation_thickness_m": 0.05,
                    "insulation_material": "polyurethane-insulation"}},
        {"source": "process_datasheet", "reference": "20-DS-014 rev 2",
         "values": {"internal_temperature_c": 45.0, "external_temperature_c": 4.0,
                    "external_film_coefficient_w_per_m2k": 300.0}},
    ],
)
assert basis.ready_for_meshing, (basis.missing_fields, basis.conflicts)

from neqsim.thermo import TPflash, fluid

gas = fluid("srk")
gas.addComponent("methane", 0.85)
gas.addComponent("ethane", 0.10)
gas.addComponent("CO2", 0.05)
gas.setMixingRule("classic")
gas.setTemperature(basis.value("internal_temperature_c"), "C")
gas.setPressure(75.0, "bara")
TPflash(gas)
state = fluid_state_from_neqsim(gas, phase="gas", velocity_m_per_s=5.0)
film = film_coefficient(state, hydraulic_diameter_m=basis.value("inside_diameter_m"))

steel = material(basis.value("wall_material"))
insulation = material(basis.value("insulation_material"))
conditions = derive_thermal_conditions(
    wall_thickness_m=basis.value("wall_thickness_m")
    + basis.value("insulation_thickness_m"),
    solid_conductivity_w_per_mk=insulation.conductivity_w_per_mk,
    solid_thermal_diffusivity_m2_per_s=insulation.thermal_diffusivity_at(25.0),
    inner_film=film,
    inner_bulk_temperature_c=basis.value("internal_temperature_c"),
    outer_film_coefficient_w_per_m2k=basis.value("external_film_coefficient_w_per_m2k"),
    outer_bulk_temperature_c=basis.value("external_temperature_c"),
)

wall = RadialConductionModel(
    [ConductionLayer("steel", steel, basis.value("wall_thickness_m"), 8),
     ConductionLayer("insulation", insulation, basis.value("insulation_thickness_m"), 20)],
    inner_radius_m=basis.value("inside_diameter_m") / 2.0,
)
steady = wall.solve_steady(
    inner_film_coefficient_w_per_m2k=film.h_w_per_m2k,
    inner_bulk_temperature_c=basis.value("internal_temperature_c"),
    outer_film_coefficient_w_per_m2k=basis.value("external_film_coefficient_w_per_m2k"),
    outer_bulk_temperature_c=basis.value("external_temperature_c"),
)
assert steady.analytic_deviation_percent < 1.0   # the closed-form check

Two dimensions, when the geometry stops being radial

from fem_coupling import (
    BoundaryCondition, ConductionProblem, FemCase, FemCouplingModel, FemMeshSpec,
    MaterialAssignment, MeshLayer, MeshSegment, read_case_results, recommend_backend,
)

print(recommend_backend(dimension=2, physics="conduction").rationale)

spec = FemMeshSpec(
    kind="axisymmetric_section",
    inner_radius_m=0.127,
    layers=[MeshLayer("steel", "carbon-steel", 0.0127, 6),
            MeshLayer("insulation", "insulation", 0.05, 20)],
    segments=[MeshSegment("upstream", 1.5, 60),
              # A flooded section: the pores fill, so it conducts like seawater.
              MeshSegment("defect", 0.4, 32, {"insulation": "flooded-insulation"}),
              MeshSegment("downstream", 1.5, 60)],
    name="20-P-001",
)
for warning in spec.mesh_warnings(max_element_size_m=conditions.max_element_size_m):
    print(warning)

mesh = spec.generate("cases/20-P-001/mesh")     # "not_executed" without Gmsh
problem = ConductionProblem.from_mesh_spec(
    spec, name="20-P-001 defect", mesh_file=str(mesh.mesh_path),
    materials=[MaterialAssignment("carbon-steel", 45.0, 7850 * 490),
               MaterialAssignment("insulation", 0.17, 700 * 1500),
               MaterialAssignment("flooded-insulation", 0.60, 1025 * 3900)],
    boundaries=[BoundaryCondition("inner", "robin", film.h_w_per_m2k, 45.0),
                BoundaryCondition("outer", "robin", 300.0, 4.0),
                BoundaryCondition("west", "adiabatic"),
                BoundaryCondition("east", "adiabatic")],
)
case = FemCase(problem, backend="scikit-fem")
case.write("cases/20-P-001")
case.run("cases/20-P-001")                       # "not_executed" without scikit-fem
results = read_case_results("cases/20-P-001")

model = FemCouplingModel()
gate = model.assess_quality(
    element_order=1, elements_across_critical_layer=20, mesh_levels=2,
    convergence_percent=0.3, max_aspect_ratio=spec.max_aspect_ratio(),
    energy_balance_error_percent=results.energy_balance_error_percent,
    biot=conditions.biot,
)
if gate.verdict != "not_usable":
    handoff = model.evaluate_thermal_handoff(
        location="20-P-001 defect section",
        heat_flow_w=results.boundary_heat_flow_w["inner"],
        reference_area_m2=3.1416 * 0.254 * spec.total_length_m,
        inner_bulk_temperature_c=45.0, outer_bulk_temperature_c=4.0,
       
        * spec.total_length_m,
    )
    print(handoff.u_multiplier)   # what the NeqSim pipeline model should carry

Cooldown and no-touch time

from math import pi
from neqsim.thermo import hydt

hydrate_c = float(hydt(wet_gas)) - 273.15
inventory = wet_gas.getDensity("kg/m3") * wet_gas.getCp("J/kgK") * pi * 0.127**2

cooldown = wall.solve_transient(
    initial_temperature_c=steady.temperatures_c,   # the shutdown profile, not a uniform one
    duration_s=48 * 3600.0, time_step_s=120.0,
    inner_film_coefficient_w_per_m2k=50.0,         # natural convection after shutdown
    inner_bulk_temperature_c=45.0,
    outer_film_coefficient_w_per_m2k=300.0, outer_bulk_temperature_c=4.0,
    inner_fluid_capacity=inventory,                # makes it a cooldown, not a shock
)
no_touch_s = cooldown.time_to_reach(hydrate_c, location="inner_fluid")

Wall stress from the metal temperatures

from fem_coupling import evaluate_wall_stress

stress = evaluate_wall_stress(
    steel, location="20-P-001 bore",
    inner_wall_temperature_c=steady.inner_surface_temperature_c,
    outer_wall_temperature_c=dict(steady.interface_temperatures_c)["steel outer face"],
    inner_radius_m=0.127, outer_radius_m=0.1397, internal_pressure_pa=75.0e5,
)
print(stress.verdict, stress.utilisation, stress.stress_category)

Three dimensions, and rendering the result

from fem_coupling import render_field, render_mesh

# The same layered section, revolved into a 90 degree wedge with two symmetry planes.
solid = FemMeshSpec(
    kind="axisymmetric_section", inner_radius_m=0.127,
    layers=[MeshLayer("steel", "carbon-steel", 0.0127, 6),
            MeshLayer("insulation", "insulation", 0.05, 20)],
    segments=[MeshSegment("upstream", 1.5, 60),
              MeshSegment("defect", 0.4, 32, {"insulation": "flooded-insulation"}),
              MeshSegment("downstream", 1.5, 60)],
    revolve_deg=90.0, circumferential_cells=12, name="20-P-001-3d",
)
mesh = solid.generate("cases/20-P-001-3d/mesh")
render_mesh(mesh.mesh_path)          # look at the geometry before solving it

# The cut planes are symmetry planes, not real boundaries.
boundaries += [BoundaryCondition("symmetry_start", "adiabatic"),
               BoundaryCondition("symmetry_end", "adiabatic")]

# ... write and run the case as above, then:
view = render_field(
    "cases/20-P-001-3d",
    scale=(1.0, 8.0, 8.0),   # a 12.7 mm wall on a 3.4 m tube needs stretching to be read
)
print(view.status, view.images)      # field_surface.png, field_slice.png, field_clip.png

Validation Checklist

• basis.ready_for_meshing is true, or every missing field and conflict is resolved and recorded.
• initProperties() ran before properties were read; density, viscosity, thermal conductivity and heat capacity are non-zero and physical.
• Material properties come from a data sheet or a certificate where one exists; a library value is marked as indicative in the assumptions register.
• The Biot number is above 0.1, or the lumped-model warning has been accepted.
• The one-dimensional steady result was computed first, and analytic_deviation_percent is below 1 %.
• Every layer carries at least three linear (or two quadratic) elements.
• mesh_warnings() is empty, or each warning has been accepted with a reason.
• The element size is no coarser than a quarter of the thermal penetration depth for the transient window of interest.
• At least two mesh levels were run, and the quantity of interest moved less than about 2 %.
• The boundary energy balance closes to better than 1 %.
• A far-field boundary sits at least three penetration depths away.
• A cooldown model carries inner_fluid_capacity; a fixed bulk temperature is a thermal-shock model, not a cooldown.
• Thermal stress was computed from the metal surface temperatures, not from the process-to-ambient difference.
• Thermal stress was assessed against a range allowable, and pressure stress separately against the primary membrane allowable.
• A three-dimensional sweep was justified by a circumferential feature or by presentation, not chosen by default over the r-z section.
• The cut planes of a partial sweep carry an adiabatic condition, and the axisymmetric weighting is off for a swept mesh.
• Any exaggerated scale or warp in a rendered figure is stated in its caption.
• The quality-gate verdict and all findings are carried into the receiving report's assumptions register.

Common Mistakes

Limitations

Linear heat conduction and species diffusion are what this skill solves, with temperature-dependent conductivity handled by Picard iteration. Radiation, convection inside the solid domain, phase change, latent heat, freezing, curing and moisture transport are outside scope. The built-in solver is one-dimensional cylindrical or planar; two-dimensional axisymmetric and plane geometry, and the three-dimensional solids swept from them, are delegated to scikit-fem or FEniCSx. The structured mesh generator produces layered sections and their sweeps only - a geometry that is not a swept layered section (a branch, a nozzle intersection, a saddle) needs an externally generated mesh, which both backends accept. A local defect is expressed as a material change over a segment, not as a change of thickness. The stress layer is linear elastic and produces membrane and gradient stresses with a category attached; it does not perform a code assessment, a fatigue evaluation or a local stress-concentration analysis, and geometry discontinuities such as nozzles, supports and welds need a local model. The thermal and mechanical problems are solved one way - a temperature field produces a stress, and the deformation does not feed back into the thermal problem. Rendering is a presentation layer with no bearing on validity: the quality gate is a screening filter, not a verification-and-validation review, a usable_with_caution verdict means any derived factor must carry an explicit uncertainty band, and a converged solve of the wrong boundary condition passes every one of its checks and still renders beautifully. Human review by a qualified thermal or stress analyst is required before a finite-element-derived number is used in a design decision. The generated scikit-fem script targets scikit-fem 8 and later, the FEniCSx script targets DOLFINx 0.8 and later, and rendering targets PyVista 0.43 and later.

Related NeqSim Functionality

• neqsim.thermodynamicoperations.ThermodynamicOperations#TPflash() and neqsim.thermo.system.SystemInterface#initProperties() - the flash and the property initialisation that supply density, viscosity, thermal conductivity and heat capacity. Reached from Python as from neqsim.thermo import TPflash, fluid, or as from neqsim import jneqsim for the full Java API.
• neqsim.thermo.phase.PhaseInterface#getThermalConductivity, #getCp, #getDensity, #getViscosity - the per-phase accessors read by fluid_state_from_neqsim.
• neqsim.physicalproperties.PhysicalPropertyHandler diffusion-coefficient models (for example Fuller-Schettler-Giddings) - the molecular diffusivities read for a species-transport model and scaled to an effective porous-medium value.
• neqsim.thermodynamicoperations.ThermodynamicOperations#hydrateFormationTemperature() (Python from neqsim.thermo import hydt) - the target a cooldown no-touch time is measured against.
• neqsim.process.equipment.pipeline.PipeBeggsAndBrills and neqsim.process.equipment.pipeline.AdiabaticTwoPhasePipe - the one-dimensional models whose U-value the finite-element multiplier corrects.
• neqsim.pvtsimulation.flowassurance.SurfCooldownAnalyzer - the flow-assurance cooldown model this skill supplies an effective U-value and a wall thermal inertia to.
• This skill does not run any NeqSim calculation itself. It reads a fluid state from a system NeqSim has already flashed, and hands U-values, multipliers and no-touch times back to NeqSim-based models.

Related Skills

• neqsim-cfd-coupling - the fluid-side companion. It resolves the boundary layer and produces the wall shear and the local velocity field; this skill consumes the film coefficient that follows and resolves the solid. The two share the same design-basis discipline and the same gate-before-quoting rule.
• neqsim-surf-cooldown-screening - consumes the no-touch time and the effective U-value from the transient model.
• neqsim-hydrate-margin-check, neqsim-wax-margin-check - supply the target temperature a cooldown is measured against.
• neqsim-pipe-wall-thickness-screening, neqsim-piping-flexibility-screening - the screening-level stress checks this skill refines with a real metal temperature gradient.
• neqsim-flow-accelerated-corrosion - consumes the metal surface temperature, which sets the film temperature that drives the corrosion rate.
• neqsim-technical-document-reading, neqsim-pid-process-operations - supply the P&ID topology, insulation specification and datasheet values that feed the design basis.

References

• Gnielinski, V. (1976). New equations for heat and mass transfer in turbulent pipe and channel flow. International Chemical Engineering, 16(2), 359-368.
• Incropera, F. P., & DeWitt, D. P. Fundamentals of Heat and Mass Transfer, Wiley
  • composite-cylinder resistance, Biot and Fourier numbers, lumped capacitance.
• Zienkiewicz, O. C., Taylor, R. L., & Zhu, J. Z. (2013). The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth-Heinemann.
• Gustafsson, T., & McBain, G. D. (2020). scikit-fem: a Python package for finite element assembly. Journal of Open Source Software, 5(52), 2369.
• Baratta, I. A., et al. (2023). DOLFINx: the next generation FEniCS problem solving environment. Zenodo. doi:10.5281/zenodo.10447666
• Geuzaine, C., & Remacle, J.-F. (2009). Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, 79(11), 1309-1331.
• Anderson, R., et al. (2021). MFEM: a modular finite element methods library. Computers & Mathematics with Applications, 81, 42-74.
• Cimrman, R., Lukes, V., & Rohan, E. (2019). Multiscale finite element calculations in Python using SfePy. Advances in Computational Mathematics, 45, 1897-1921.
• ASME Boiler and Pressure Vessel Code, Section VIII Division 2, Part 5 - classification of primary, secondary and peak stresses.
• ASME B31.3, Process Piping - displacement stress range and the allowable stress basis for restrained thermal expansion.
• Roache, P. J. (1994). Perspective: a method for uniform reporting of grid refinement studies. Journal of Fluids Engineering, 116(3), 405-413.
• ASME V&V 10-2019, Standard for Verification and Validation in Computational Solid Mechanics.