Hat Stiffened Panel Ansys

A B S T R A C T

The linearized buckling load of a stiffened panel depends on the stress stiffness its shell element assembles. We read it from exported operators against three truncations of one second variation. The classic pass of ANSYS SHELL181 carries a rotation-rotation block pairing the drilling freedom with the bending rotations and its perturbation pass does not; removing the block recovers the perturbation load factor to 0.02%. SHELL281 carries block and couplings in both passes. Abaqus S4 matches the critical mode of the complete second variation to 1.0000 on the optimized panel under cantilevered shear a 20-node continuum lies 3% to 6% above that form, S4 and SHELL281, 11% and 23% below both SHELL181 passes and 25% above Abaqus S8R, at the finest meshes. On a conventionally stiffened panel the SHELL181 passes stand 1.0% and 3.5% above the complete form, 9% and 21% at half the rib pitch; under a shear flow, on cylinders, open beams and under axial compression the three forms coincide and no pass parts by more than 0.3%.

Ntroduction

Stiffened shells carry a large part of the primary structure of aircraft, ships, pressure vessels, cranes and rail vehicles, and where the skin is thin the design is governed by stability rather than by strength. The stiffeners decide how the skin is subdivided, and therefore what buckles first and at what load, so a design method for such a structure stands or falls on the buckling analysis inside it. In nearly all of that work the analysis is a linearized eigenvalue problem: a pre-stress, a stress stiffness assembled from it, and the load factor at which the two stiffnesses cancel on some mode.

Every commercial program offers it, and most layout optimization of stiffened panels that carries a buckling constraint builds it in [1, 2, 3, 4]. The stress stiffness of a built-up section is not fixed by the continuum. A continuum element has three translations at a node and its stress stiffness follows from the second variation of the strain measure without a choice being made; a shell element assembled for a stiffened panel has six, and whether the sixth of them takes part in the shell director decides an entire block of the operator. Made one way, the drilling freedom pairs with the two bending rotations under the moments and the transverse shear forces of the pre-stress; made the other, that block does not exist. The sixth freedom cannot be dispensed with, because a rib meets a skin at a right angle and the rotation that is drilling for one of them is bending for the other. A genuine drilling rotation goes back to Allman , its variational footing to Hughes and Brezzi , and the families of element built on either are reviewed by Boutagouga ; none of the Technology Development Program Project of China Railway Eryuan Engineering Group Co., Ltd. (grant number: KSNQ253005).

Element-Dependent Buckling Of Stiffened Panels

buckling optimization studies cited above writes its stress stiffness out, and commercial documentation does not state the operator either. The usual benchmarks cannot reveal the choice: a flat plate under a uniform membrane pre-stress carries neither moments nor transverse shears, so every form passes it.

Nor can the question be settled by agreeing with a commercial program, because a commercial program is not one authority but several. Its classic eigenvalue pass and its linear perturbation pass from a nonlinear base state are two eigenproblems that need not assemble the same stress stiffness, its other shell elements may assemble yet others, and a second program brings its own. An implementation that agrees with one load factor has established which operator it shares and nothing more, and a design checked against one program has been checked against one of its operators.

The question a stiffened-panel analysis has to face is therefore not which program to trust but how far the element formulation and the eigenvalue procedure move the linearized buckling load, where in the operator the difference sits, and on what structures it appears at all.

Examining that needs three things that the literature does not supply together. The candidate stress stiffnesses of a six-freedom shell must be written out and their assembly certified independently of any program, so that a comparison is between named operators and not between numbers. The operator a commercial formulation actually assembles must be read from the formulation itself, from its exported matrices where a program exports them and from its load factors, modes and energies where it does not, rather than inferred from agreement. And the comparison must be placed against a reference whose exactness rests on none of the shell operators, and run on structures where the candidates part as well as on structures where they do not, so that the condition under which the choice matters can be stated.

This paper does the three. It derives, for a four-node shell element with a director, the three truncations of one second variation of the pre-stress work that the treatment of the sixth freedom admits, and certifies their assembly against an exact identity on that work. It reads the operators of ANSYS SHELL181 and SHELL281 from the programs’ own exports, in both passes, places Abaqus S4 and S8R by load factors, critical modes and energies on the same nodes, and adds a reference that involves no shell stress stiffness at all, a continuum model of the same structure in 20-node axial, shear and combined loads, and over cylinders, open beams, a box girder and a published strip , and asks what feature of the pre-stress and of the critical mode decides where the truncations part. Two words are used in a fixed sense from here on: a pass is one of a program’s two linearized buckling procedures, the classic eigenvalue pass or the linear perturbation pass from a geometrically nonlinear base state, and a form is one of the three truncations assembled here, membrane, block and full, defined in Section 2.2; every other term of art is defined where it first enters.

The optimized panels are produced by the stiffener layout optimizer of the companion paper on a grid of candidate lines; the optimizer is the instrument here and not the subject, its designs entering as fixed structures that both programs can analyse.

Section 2 states the element, the buckling problem and the three forms; Section 3 the structures; Section 4 the identity and the rules of the cross-code comparison; Section 5 the two operators of SHELL181, the form each commercial formulation assembles, and when the forms part; Sections 6 and 7 discuss and conclude, and Appendix C records what a repetition needs.

The Element

of freedom per node. The membrane, bending and transverse shear rigidities are 𝑡𝐃, (𝑡3∕12)𝐃and 𝑘s𝐺𝑡𝐈2 with the shear correction 𝑘s = 5∕6, where 𝑡is the plate thickness, 𝐃the plane-stress constitutive matrix of the isotropic material, 𝐺 its shear modulus and 𝐈2 the 2 × 2 identity. Membrane and bending terms use 2 × 2 Gauss quadrature. The transverse shear strain is taken from the assumed field of the four edge midpoints , integrated by the same rule, and the stress stiffness takes from that same field both the second-order shear strain its shear-weighted terms are formed on and the frozen shear force that weights them, so that on every term the two matrices are variations of one energy; that is the one choice inside the element that decides whether its stiffness and its stress stiffness are two derivatives of one energy, and it is made here so that they are.

The drilling degree of freedom is held by a spring scaled with the bending rigidity rather than the membrane rigidity, a diagonal term of 10−4 tr[(𝑡3∕12)𝐃] on each node’s drilling freedom. Scaling that spring with the membrane rigidity is admissible for a single plate, but in a built-up section the drilling axis of a stiffener wall coincides with a

Element-Dependent Buckling Of Stiffened Panels

global bending rotation of the panel, so a membrane-scaled spring absorbs energy of the global modes and inflates the load factors selectively. This is the element of every result below, and two variants of it appear beside it. The first, the control variant, takes the transverse shear of the stiffness and the frozen shear force at one central point instead of from the assumed field; the optimized designs were produced under it, and it is reported wherever it changes a reading.

The second variant is used for the tied conforming panel and for the cross-code checks. It differs from the element above in two settings: the membrane carries the two incompatible modes of Taylor et al. , condensed at element level so that the element interface is unchanged, and the drilling freedom is held by the Hughes–Brezzi penalty instead of the grounded spring. The bilinear membrane cannot represent in-plane bending without a parasitic shear strain, which the incompatible modes remove, and ANSYS SHELL181, the element whose two passes the paper sets side by side, offers the same choice through KEYOPT(3); the second variant is therefore the one nearest the program’s own membrane, and the cross-code comparison of Section 4.2 is made between matching formulations rather than across them.

The stiffener walls are meshed independently of the skin in the in-plane direction, each design segment being subdivided into 𝑛sub elements along its length and 𝑛𝑧elements through its height. The nodes at the foot of a wall that fall between two skin nodes are tied to those nodes by multipoint constraints acting on all six degrees of freedom,

(1)

with 𝒖𝑑the six freedoms of the dependent wall node, 𝒖1 and 𝒖2 those of the two skin nodes it lies between, and 𝜉its normalized position on the skin edge. The constraints are imposed by a penalty 𝛼= 𝛼𝑓𝑘ref with 𝛼𝑓= 103, the low end of a range over which the leading load factor moves by 0.07% as 𝛼𝑓is carried to 106. The wall is tied to the mid-surface of the skin, so the eccentricity of half a skin thickness between the two mid-surfaces is not carried. On the panels of this paper that offset is 5 mm against a wall height of 100 mm, and what it would add is a membrane-bending coupling of the junction; a section in which the skin thickness is a larger fraction of the wall height would need it modelled rather than tied flat.

The optimized designs come from a density-based layout optimizer, and their analysis model keeps its form: every candidate rib segment carries a densitȳ 𝜙𝑒that scales its elastic stiffness by 𝜀𝐸+(1−𝜀𝐸)̄𝜙𝑝

Bȳ 𝜙𝑝𝐺

𝑒, with 𝑝= 𝑝𝐺= 3 and 𝜀𝐸= 10−9. The designs analysed here are the optimizer’s converged fields cut at a threshold, so every element is either at full density or void, and the void elements remain in the model at 𝜀𝐸. The geometric stiffness takes the same exponent as the elastic one because taking 𝑝𝐺= 1, a common choice, makes the ratio of geometric to elastic stiffness of a void element grow as̄ 𝜙1−𝑝

𝑒

and fills the spectrum with modes localized in the void. The reference stiffness 𝑘ref used by the penalty, and by the stabilization below, is the largest full-density diagonal entry of the assembled stiffness matrix among the degrees of freedom of the skin nodes. This choice is deliberate. The skin is always solid and its elements are geometrically regular, so 𝑘ref is a property of the mesh alone and does not move with the design. The global maximum diagonal entry, which is the natural first choice, is instead set by the stiffest membrane term anywhere in the model, and a stiffener element of high aspect ratio has a short-direction membrane diagonal of order 𝐸𝑡(𝐿∕ℎ), two orders of magnitude above that of a regular element. Taken globally the scale follows the most distorted element in the model, which on one design put the spring at 15% of the skin’s bending stiffness and the load factor at nearly twice SHELL181’s on the same mesh (Appendix C).

(2)

is applied to the rotational degrees of freedom , and to those only, 𝐞𝑖being the unit vector of freedom 𝑖. A grounded spring on the translational degrees of freedom would act as an elastic foundation whose stiffness per unit area grows with mesh density, and it makes the load factors diverge under refinement. The value of 𝜀𝑟is not free either: it has to sit above the conditioning floor of the void elements’ rotations, whose stiffness is scaled by 𝜀𝐸, and below the level at which it acts as a rotational foundation on the structure itself. Appendix C gives the sweep that places it, two decades above the floor.

Element-Dependent Buckling Of Stiffened Panels

2.2. The linearized buckling problem and the block The equilibrium displacement 𝒖of the reference load 𝒇follows from 𝐊𝒖= 𝒇, and the linearized buckling problem

(3)

with 𝐊the stiffness, 𝐊g(𝒖) the stress stiffness assembled from that state, 𝜆the load factor and 𝝋the buckling mode, which we solve in the inverse form 𝐊g𝝋= −𝜇𝐊𝝋with 𝜇= 1∕𝜆, so that the critical modes are those of largest |𝜇|. The geometric stiffness of a stiffened panel needs care that a flat plate does not require. A stiffener wall stands normal to the skin, so displacements that are in-plane for the wall are out-of-plane for the panel and the reverse, and a stress stiffness assembled only from the transverse displacement gradients misses the terms that carry the local instability of the wall. We therefore retain all three displacement components. With 𝐍= [𝑁𝑥𝑥, 𝑁𝑦𝑦, 𝑁𝑥𝑦]𝖳 the membrane forces sampled at the Gauss points of the element from the equilibrium state, 𝑔𝑥

𝑎the Cartesian Shape

function gradients of node 𝑎in the local element frame, and 𝑤𝑞the quadrature weight of Gauss point 𝑞including the Jacobian, the element geometric stiffness takes the Kronecker form

(4)

the Kronecker product with the identity on the translational components expressing that the same scalar coupling 𝑀𝑒

𝑎𝑏

between nodes 𝑎and 𝑏acts on each of the three directions and nothing acts on the rotations. The form is invariant under rotation of the element frame, which is what allows one expression to serve the skin and the walls of any orientation, and it reduces to the familiar plate expression when the membrane forces of the walls vanish.

Equation (4) is one truncation of the stress stiffness and not the only one, and the choice between them is the subject of this paper. It is the customary one for a flat shell element assembled from a membrane and a plate, whose geometric stiffness has been built from the membrane forces alone since the elements of the 1980s, as the review of Gal and Levy records . We call it the membrane form. The element carries six freedoms at a node (Fig. 1): three translations, two bending rotations that tilt the director, and a rotation about the normal, the drilling rotation, which leaves the director where it is at first order. Whether that sixth freedom takes part in the director at second order is the choice, and a built-up section does not let it be avoided: where a rib wall meets the skin at a right angle the two plates share one set of nodal rotations, and the rotation that is drilling for one is a bending rotation of the other. If the drilling rotation takes part in the shell director, a second-order term of the director carries it.

For a director 𝐝that is the reference normal 𝐞3 rotated by the exponential map of the rotation vector 𝜽= (𝜃𝑥, 𝜃𝑦, 𝜃𝑧),

(5)

and its in-plane components enter the curvature 𝜅𝛼𝛽= 1 2(𝐚𝛼⋅𝐝,𝛽+ 𝐚𝛽⋅𝐝,𝛼) and the transverse shear 𝛾𝛼= 𝐚𝛼⋅𝐝, 𝐚1 and 𝐚2 being the base vectors of the deformed mid-surface, so that curvature and transverse shear each acquire a term in which the drilling rotation multiplies a bending one. Contracting with the frozen resultants, the stress stiffness acquires

(6)

over the mid-surface 𝐴, pairing the drilling freedom with the two bending rotations and weighted by the moments 𝑀𝛼𝛽and the transverse shear forces 𝑄𝛼of the pre-stress, not by the membrane forces times 𝑡2∕12 as the continuum truncation written on displacement gradients would have it. We call Eq. (4) together with Eq. (6) the block form.

The coefficient of Eq. (6) is a property of the coordinates chosen on the rotation group at second order, and not of the director alone. Composing the same rotation as a rotation about the reference normal followed by one in the tangent plane leaves 𝐞3 untouched by 𝜃𝑧and gives the block coefficient zero; the reverse order gives it one; the exponential map, Eq. (5), gives it one half. The three agree at first order and differ by a redefinition of the bending rotations by terms 2𝜃𝑧𝜃𝛼, and a quadratic change of parameters 𝜽↦𝜽+ 𝒒(𝜽) changes the second variation at the reference configuration by the first variation contracted with 2𝒒, which is the pairing of the frozen moments and shear forces with 𝜃𝑧𝜃𝛼, Eq. (6) itself. The block and that non-invariance are one object. The bifurcation load of the geometrically exact problem is

Element-Dependent Buckling Of Stiffened Panels

unaffected, because there the second variation is taken at the buckling configuration, where the first variation of the total energy vanishes and a change of parameters adds nothing; the linearized problem takes its stress stiffness at the unstressed reference, where the first variation of the frozen work does not vanish, which is why its load factor moves.

Everything derived here is for the exponential map, on which the identity of Section 4.1 is written. That the rotational part of a geometric stiffness depends on how rotations are parameterized is established well beyond shells: for beams in space by Argyris and co-workers [13, 14] and Yang and McGuire , and in corotational formulations through the moment correction geometric stiffness [16, 17, 18, 19], whose leading term at the reference

Configuration, −1

2spin(𝒎), is the closest published counterpart of Eq. (6), although it is antisymmetric, written for nodal moments of any element and without a shear part. Simplified geometric stiffnesses that drop the rotations outright are in use for thin-walled structures . Geometrically exact theories that carry the drilling angle in the rotation tensor , or give it to a shell as a redundant micropolar freedom , carry a rotational geometric block weighted by resultants and couples , of which Eq. (6) is the structural counterpart for a facet element with frozen resultants; where the drilling freedom is supplied variationally through the skew part of the membrane strain , or appended to a five-freedom shell by a Lagrange multiplier , it enters neither the curvature nor the transverse shear and the second variation carries no term in it. Neither the term nor its absence is new; what is not on record is which of them an analysis inherits.

The same second-order director implies two couplings between the translation gradients and the rotations,

(7)

with 𝐯the translation of the mid-surface and 𝛼, 𝛽the in-plane indices, the product of the linear part of a base vector with the linear part of the director, the factor of two arising because that product carries no one half where the director’s own second-order term does. The form that carries Eq. (4), Eq. (6) and both couplings we call the full form; it is the complete second variation of the pre-stress work for the exponential-map director and the assumed shear field’s extension, and Section 4.1 certifies it on every vector. Neither coupling carries the drilling rotation, so they are terms any shell with a director may have, and they are not equally consequential: on the structures of this paper the moment-weighted one is inactive and the shear-weighted one moves the load factor by a quarter to a third, Table 5 in Section 5.2. The shear-weighted coupling is formed, like the shear part of Eq. (6), on the nonlinear extension of the assumed shear field of the four edge midpoints, which the assumed-strain method does not fix. All three forms are therefore carried to the results, and which of them a commercial formulation assembles is measured there rather than assumed.

Both Branches Of The Spectrum

Eq. (3) admits eigenvalues of both signs, 𝜆< 0 being buckling under the reversed load, and the eigensolvers used for buckling constraints are usually asked for the eigenvalues of largest magnitude of the inverse form 𝐊g𝝋= −𝜇𝐊𝝋, 𝜇= 1∕𝜆, or for those of algebraically largest 𝜇, which returns the positive branch alone. For a compressive edge load the second request is correct, the reversed load being tension. For a shear load the reverse is an equally admissible service load and the spectrum of a symmetric panel is symmetric about zero.

Every shear and combined case of this paper is therefore solved for the modes of largest |𝜇| on both branches, load factors are compared between programs on their magnitudes, and the sign pattern of the spectrum is itself one of the checks of the cross-code comparison.

Three rules follow from that and are used throughout. A load factor quoted for a design is the critical one by magnitude, whichever branch it lies on, and on shear design A that is the reversed branch on both sides. A multiplier returned by a perturbation run scales the perturbation load, so the load factor of its branch is 𝛼plus the multiplier, formed in that order and only then taken in magnitude, a reversed branch being 𝛼minus the multiplier’s magnitude.

And a mode compared with another mode is taken on one branch for both, the branch on which the model’s critical mode lies, so that a mode is never set against the mirror image of its counterpart under the reversed load.

The Design Sensitivity Of The Block

The eigenvalue sensitivity of a design-dependent structure carries three terms, the explicit derivative of the stress stiffness, the elastic term, and an adjoint term through the equilibrium state that the pre-stress depends on, the third being the one dropped when the pre-stress is treated as design-independent. What the block adds is its own adjoint load, the derivative of 𝝋𝖳𝐊g(𝒖)𝝋with respect to the state, which for the moment part is formed through the curvature operator and for the transverse shear part at the four tying points of the assumed shear field. The block’s element form

Element-Dependent Buckling Of Stiffened Panels

is written in Appendix A and both sensitivities in Appendix B, because they are what an implementation that carries the block has to get right, and the block form and the full form pass a finite-difference check over the whole chain of assembly, static solve and eigenproblem on a 12 × 4 panel, with a central step of 10−6 on six element densities; the eigenvalues checked are separated, and the derivative of a repeated eigenvalue is not defined by these expressions.

The Structures

The stiffened panels are built on the rectangular panel of Fig. 2, 𝐿𝑥×𝐿𝑦= 3×1 m, skin thickness 0.01 m, stiffener thickness 𝑡s = 0.008 m and height 𝐻𝑧= 0.1 m, with 𝐸0 = 2.1 × 1011 Pa and 𝜈= 0.3. The edge 𝑥= 0 is clamped in all six degrees of freedom in every case, and a resultant of 105 N is applied to the skin nodes of the edge 𝑥= 𝐿𝑥, in one of three directions that define the three load cases: along −𝑥, which is axial compression; along +𝑦, tangential to the edge, which is the cantilevered shear: it puts the panel under a mean shear 𝑁𝑥𝑦∼𝐹∕𝐿𝑦together with the in-plane bending that a load carried to a clamped edge implies. The two words are kept apart throughout: a cantilevered shear is this case, and a shear flow means the self-equilibrated tangential tractions on all four edges of Section 5.3, which carry no in-plane bending; or both at once, 105 N along each of the two directions, which is the combined case. The three cases are chosen for what they do to the pre-stress: compression leaves the plates in a nearly pure membrane state, shear bends the ribs sideways as flanges of a panel bent in its own plane, and the combined case carries both parts at once, which is exactly the gradation the operator question of Section 5.2 needs.

The reference load leaves every case elastic: the axial membrane stress is 10 MPa on the skin alone, 6 MPa with the stiffeners sharing it, and the in-plane bending of the cantilevered shear puts about 100 MPa at the clamped root once the stiffeners act as flanges. The axial case stays elastic up to its critical load, near 64 MPa at 𝜆≈11, well below the 235 MPa yield stress of a Q235 structural steel (S235 in EN 10025), so its buckling factor is a load the structure can elastically reach. Under shear the root of that steel would yield at 𝜆≈2.4, roughly a third of the elastic critical load, and a Q355 (S355) root at about half of it, so the shear and combined factors are read as measures of the operator gap, a ratio the load level does not move.

Four optimized designs appear in the tables (Fig. 5): one under axial compression, one under the combined load, and two under the cantilevered shear. Each is the discrete layout, cut at the threshold that conserves the material of the converged field, that the layout optimizer of the companion paper returned on a 48 × 16 grid of candidate stiffener lines at a stiffener volume of 0.4 of the fully stiffened panel. That paper reports the optimizer, the projection and the threshold, and it confines its own design study to axial compression, where the two passes agree and the choice of stress stiffness does not enter; nothing below depends on it, the designs entering here as fixed structures. The two shear designs differ only in the operator they were optimized under. Design A was produced with the rotational block carried in the stress stiffness, design B with the membrane form; their layouts are close and their load factors under a common measure differ by a fraction of a per cent, so the pair also measures how much the choice of operator moves the design itself, as distinct from the number reported for it.

The stiffener spacing of the grid is 0.0625 m, the walls are divided into 𝑛𝑧= 3 elements through the height and 𝑛sub = 2 along each design segment, and 𝑝= 𝑝𝐺= 3. The designs on this mesh, with their rib feet tied to the skin as in Section 2.1, are the tied models of the tables; they carry one element per design cell of skin and six per segment of rib, 10 368 candidate elements and 9 732 nodes in all, of which elements 3 456 are active on the conforming panel and 4 614 on design A.

Counted on that lattice, design A carries 641 rib segments, 40.06 m of rib meeting at 383 lattice nodes, of which 161 are crossings, 198 T junctions, 17 corners, 4 free ends inside the panel and 3 points where a straight rib passes through a lattice node; designs B, the combined design and the axial design carry 640 segments with 164, 164 and 109 crossings and 2, 5 and 8 free ends, the axial layout standing apart in its 90 T junctions against the shear layouts’ 190 and above.

The shear and combined cases, whose load can reverse, are solved on both branches of the spectrum, Section 2.3; the axial case, whose reverse is tension, on the positive branch alone. On the combined case the reversed branch reverses the whole load vector, compression with shear, so it is the panel pulled and sheared the other way and not the same compression with the shear reversed; the latter is a separate load case, which is not solved here, and the reversed branch is carried only so that the critical magnitude is not missed.

Further structures enter at three levels. The conforming panel is the same panel conventionally stiffened, its ribs on a 250 mm pitch in both directions under the cantilevered shear, laid out on the 48 × 16 grid with a rib on every fourth line so that its skin and ribs carry the designs’ elements and ties; in our implementation it is analysed with

Element-Dependent Buckling Of Stiffened Panels

the incompatible-modes membrane and the Hughes–Brezzi drilling penalty, the variant closest to SHELL181 with KEYOPT(3) = 2, and its load factors are marked as such; its rebuild without ties, described next, like every such rebuild and every new structure, carries the compatible membrane and the drilling spring instead. Its 448 segments, 28.00 m of rib, run edge to edge: 33 crossings, 28 T junctions at the boundary and no free end inside the panel. It is a different structure from the densely stiffened panel of Table 2, which carries a rib on every line of a 24 × 8 grid and enters only as the axial verification case.

For the second commercial program and the second element order, the conforming panel, shear design A and the axial design are rebuilt without ties, the shared-node rebuilds of the tables: skin and ribs share nodes along every rib foot and at every crossing, each design cell being one eight-node element or two by two four-node elements on the same node lattice, with the supports and the consistent edge loads of the tied models.

On the same material, three new structures of Fig. 3 are built the same way. The cylinder has radius 0.5 m, length 1.0 m and wall 5 mm, is clamped at 𝑥= 0, and carries at its free end, on the shell alone, either an axial resultant of 5 × 106 N or a torque of 106 N m; its stiffened variant adds 24 internal stringers and rings at the quarter points, blades 40 mm deep and 4 mm thick.

The I beam, 400 by 200 mm between flange mid-planes with 12 mm flanges and an 8 mm web, spans 4 m on forks that hold the web ends laterally and vertically, and carries end moments of 105 N m applied as linear edge tractions, for which the closed-form lateral-torsional factor is 4.7099.

The channel, 300 by 100 mm with 10 mm flanges and a 6 mm web, is a 2 m cantilever with 2 × 104 N along the web at its tip, away from the shear centre, so that its pre-stress carries torsion as well as bending. The cylinders carry 96 × 32 cells and the beams 25 mm cells.

The box girder, Fig. 4, is a generic girder in the proportions of a container-crane main girder, every dimension a round number: 35 m between bearing centres with 1.0 m overhangs, a 1.2 by 2.2 metre section, 14 mm flanges overhanging the webs by 0.15 m, an 8 and a 10 mm web, and thirty-four interior diaphragms of 8 mm plate at a 1.0 m pitch, each with a 0.8 by 1.8 m manhole, and end plates of the same 8 mm plate, solid, at the bearing centres. It is seated on 2.0 m bearing pads, each holding the lateral translation over its length and the vertical translation along its centre line and turning about that line, the left one also holding the axial translation at its centre, and carries two wheel loads of 200 kN, 9.0 m apart astride midspan on the top flange above the 10 mm web, together with its own 261 kN at 7850 kg/m3 and 9.81 m/s2. It is meshed at 333 × 24 × 44 divisions along its whole 37 m, across the width and through the depth, 315 of them between the bearing centres and nine in each overhang, 67 716 shell elements once the manholes are cut on grid lines, with the assumed shear field in both matrices; Section 5.3 reports what the block does to it.

Erification

Two questions are asked of the analysis and they are not the same question. Whether the operators assembled here are the ones their derivation implies is a matter of internal consistency and admits an exact answer, Section 4.1. Whether a model reproduces what a commercial program computes for the same structure is a matter of agreement and admits only a measured one, under the rules of Section 4.2. Neither answers which operator a commercial program assembles; Section 5 answers that from the programs’ own operators and modes.

4.1. The identity the stress stiffness must satisfy The element writes its internal work on the stress resultants against the strain measures of a surface carrying a director, and its stiffness is the second derivative of that work at the undeformed state. The stress stiffness must then be the second derivative of the same work, with the resultants held at the values the pre-buckling state gives them:

(8)

where 𝐯is a perturbation of the nodal freedoms, 𝑠a scalar carrying it, 𝑁𝛼𝛽, 𝑀𝛼𝛽and 𝑄𝛼the membrane forces, moments and transverse shear forces of the pre-buckling state, and 𝜀, 𝜅and 𝛾are the full nonlinear strain measures of that kinematics and not their linear parts. Whichever assembled form reproduces Eq. (8) is the consistent one, and the others are inconsistent with the very energy whose stiffness they are paired against in the eigenproblem. No commercial program enters, and no question of which form is more accurate: this is an identity or it is not.

Element-Dependent Buckling Of Stiffened Panels

Two properties make it a test rather than another approximation. 𝑊is a scalar, so nothing has to be assembled to evaluate it. And 𝑊is at most cubic in 𝐯for the kinematics the stress stiffness is written on, which is the director truncated at the second order of Eq. (5); the exponential map itself is not polynomial, and the identity is a statement about that truncation and its second variation, not about the full rotation. With the resultants frozen and each truncated strain measure at most cubic, the product of a base vector that is linear in the displacement with a director that is quadratic in the rotations, so the central second difference [ 𝑊(ℎ𝝋)−2𝑊(𝟎)+𝑊(−ℎ𝝋) ]∕ℎ2 is exact for any step, the cubic part being odd and cancelling. The usual compromise between truncation and round-off does not arise, and the step is a check rather than a parameter: evaluated at ℎ= 10−2, 10−1 and 1 the difference agrees to thirteen significant figures.

The reference is built from the kinematics alone, sharing with the assembly only the definition of the element frame, since reusing the assembled blocks would prove nothing. The transverse shear part of that kinematics has to be the element’s own. The element samples the covariant shear strains at the four edge midpoints and interpolates them , so two things in 𝑊must come from that field: the 𝛾𝛼it is written on, carried to its full nonlinear form through the same tying interpolation, and the frozen shear force 𝑄𝛼, which is the force that field produces at each Gauss point under the pre-buckling state.

Each choice is testable, because a mismatch shows as a failure of the identity on the rotations. Written with the pointwise strain measure, the reference misses the assembled block by tens of per cent; written with the assumed field but with the frozen force taken at the element centre, where a one-point rule would take it, the block form misses by 13% on the rotations alone. With both from the assumed field the identity closes to the digits of Table 1, and that is the element carried through the results.

Table 1 reports it on vectors chosen so that each isolates one part of the operator, on a small panel of the same element and the same kind of load as the designs. Its entries are relative differences from the reference column, except where the reference vanishes and the value itself is given, and the test vectors are components in the global frame, so the block form falls short by the two couplings it sets aside; with both carried, those rows read 2×10−11 and 2×10−13.

The critical mode of that panel is that panel’s, and what the block is worth on the critical modes of the designs themselves is a different number, reported with the load factors in Section 5.2. The continuum truncation of the table is Eq. (4) applied to the three rotations as well as the three translations, the membrane forces weighted by 𝑡2∕12, which is what a thickness-integrated continuum element written on displacement gradients produces.

On the translations every form is exact, the membrane block being common to all of them. On the rotations the membrane form is short by the whole of the term, having nothing there at all, and the continuum truncation carries it with the wrong sign.

The last row is the most telling. A drilling rotation alone leaves the work exactly stationary to second order, because

𝑥−𝜃2

𝑦) and every in-plane component of it carries a bending rotation as a factor; the continuum truncation nonetheless returns a drilling stress stiffness there. It is not merely an incomplete approximation, but introduces an unphysical artificial stiffness.

The two rows that mix translations with rotations are where the block form is not exact: it lacks the two couplings of Eq. (7) and falls short by them, by two tenths of a per cent on the critical mode of that small panel, while the full form is exact there to eleven digits. Two tenths of a per cent on one vector is not a measure of what a coupling does to an eigenvalue, and the shear-weighted coupling moves the load factors of the designs by a quarter to a third (Section 2.2).

The identity certifies each form as the truncation it claims to be; which truncation a structure follows is not a question it can answer.

The Commercial Programs

The first program is ANSYS Mechanical APDL with SHELL181, full integration with incompatible modes, on decks exported from the optimizer’s model element for element. Where no operator question arises the two models agree: on a densely stiffened panel, a rib on every line of a 24 × 8 grid, under axial compression, Table 2, our incompatible-modes variant reproduces the program to 0.5% on six modes and 0.01% on the compliance, and the compatible element carried through the paper sits 1.3% to 1.9% above it, inside the 4.3% by which the program’s own two integration options differ, both measured against the full-integration option as Table 2 is. On the axial design, refined until both have converged, the leading load factor agrees to 0.16% and the paired skin modes correlate to 0.9999; on the optimization mesh of the tables the offset is 2.0%, the compatible membrane’s.

Table 1

The assembled forms against the identity of Eq. (8), on a 0.25 m square panel whose stiffeners run on the interior lines of a 4×4 grid of cells, under the cantilevered shear of the designs. The column headed continuum is the truncation of Eq. (4) applied to all six freedoms, not the solid model of Section 5.2.

Table 2

The analysis model against SHELL181 on the densely stiffened panel under axial compression. Both difference columns are taken against the column headed full.

−0.0%

The second program is Abaqus/Standard, with S4, full integration with enhanced membrane strains, and S4R, reduced integration with hourglass control; ANSYS SHELL281 and Abaqus S8R are the eight-node elements. Every deck of either program is written from one mesh of ours, node for node: the tied models with their ties as constraint equations, the shared-node rebuilds and the new structures with shared nodes.

Before any operator is compared the models are checked where none can differ. Under axial compression Abaqus S4 returns 6.1494 on the conforming panel and 10.713 on the axial design against SHELL181’s 6.1453 and 10.7018, and S4R returns 10.189 against the reduced SHELL181’s 10.217, so S4 pairs with the full-integration option and S4R with the reduced one to 0.3%.

Under shear the reduced options behave as the full ones do: SHELL181 with KEYOPT(3) = 0 parts its passes by 2.39% and 14.33% on the rebuilt conforming panel and design A at one subdivision, where the full-integration option parts them by 2.39% and 13.98% on the same two meshes, and S4R returns 4.5655 and 6.8516 with its passes equal, 0.9% and 1.4% below S4 on the same meshes, so the block is not an artefact of an integration rule. Releasing the rotations from every rib-to-skin tie moves shear design A in Abaqus by 0.16%, so the ties are not where the programs could part. Both ANSYS releases used, 19.2 and 2024 R2, return every load factor of both passes to within two parts in a million (Appendix C). None of these values enters a table: together they bound what the model, the element option, the ties and the release can contribute to a cross-program difference, 1.4% at the most and a few tenths of a per cent in every other check, against the 14% to 33% the operators are about to be found to carry, so a gap of that size cannot be laid to any of them.

Three rules make these comparisons controlled. The exported element set is the one analysed, cut from the frozen design by the same threshold and the minimum of a segment’s two endpoint values; a reanalysis that updates the design before writing its deck exported a structure 144 elements away and a load factor three times off. The resultant load is identical, applied to skin nodes only, and the exporter counts the load terms it would have to skip and requires the count to be zero. And load factors are compared by magnitude on both branches, mode by mode where the modes are separated and as pairs where they are not, with the compliance taken as twice the strain energy the program reports.

Results

Three conventions are used for percentages. A value 𝑋is said to lie 𝑝% above or below a reference 𝑌with 𝑝= 100 |𝑋∕𝑌−1|, the reference named each time. The gap between two passes is (𝜆classic −𝜆pert)∕𝜆classic, and the worth of the block is its effect as a fraction of the value without it. Every value states its mesh: the tied models on the optimization mesh, or the shared-node rebuilds at a stated number of subdivisions of the design cell.

Every tabulated load factor is the critical one by magnitude on the branch stated in Section 2.3: the reversed branch on the shear and combined cases of the optimized designs and of the conventional panels, the positive branch under axial compression and on the cylinders, the beams and the shear flow, whose two branches agree to the digits printed.

Where a mode is compared with a mode, including every value of the modal assurance criterion below, both are taken on the reversed branch of the same model. A formulation is said to carry a term when the block of its exported stress stiffness that the term would occupy is not only present but acts on the structure’s critical mode, measured by the quadratic form 𝝋𝖳𝐁𝝋as a fraction of 𝝋𝖳𝐊g𝝋 and by the norm ratio ‖𝐁𝝋‖∕‖𝐊g𝝋‖; a block of any norm whose quadratic form on that mode vanishes changes no load factor through it, and both quantities are reported wherever presence and effect part company.

The Two Operators Of Shell181

The classic pass takes a linear static state 𝒖under the reference load and solves

(9)

with 𝐊0 the stiffness at the undeformed configuration, 𝐆c the stress stiffness it assembles from that state, and 𝜆and 𝝋 the load factor and the mode. The perturbation pass converges a geometrically nonlinear state 𝒖𝛼under 𝛼times that load, regenerates the element matrices there, and solves

(10)

with 𝐊T the tangent stiffness at that state, 𝐆p the stress stiffness of the perturbation load Δ𝒇, the reference load applied once, 𝛽its multiplier and 𝝍the mode, so that the load factor is 𝛼+𝛽. The solver file of each run holds its own pair, the stiffness under STIFF and the stress stiffness under MASS, and both come from one program and one model, so they share an ordering. At 𝛼= 0.01 the two states differ by a hundredth of the reference displacement; no claim is made that 𝐆p is the geometric part of 𝐊T.

Table 3 sets the two stress stiffnesses side by side, the rows labelled by the program’s own mapping file or, on the tied models, by the jump in the diagonal of the stiffness, which on both grids of the conforming panel selects exactly the rows the mapping names. The translational blocks of the two SHELL181 passes agree to one or two parts in a thousand, measured as the norm of their difference. The rotation-rotation block is present in the classic pass and absent from the perturbation pass, 10−8 to 5×10−7 against 102, on every model, designs B and the combined design repeating design A to three figures. The classic block has no diagonal and no trace, the structure Eq. (6) produces.

Set against the block assembled here on the same rows of the conforming panel, and labelled in the global frame the mapping uses, its drilling-drilling part vanishes in both and its 𝜃𝑧-to-𝜃𝑥, 𝜃𝑦part agrees in norm to 6%, while the 𝜃𝑥, 𝜃𝑦part is 1.9 times the program’s.

The frame matters for reading those names. Equation (6) pairs each plate’s drilling rotation with its own bending rotations, and the skin lies in the 𝑥-𝑦plane, so on the skin that pairing is 𝜃𝑧against 𝜃𝑥, 𝜃𝑦; a rib wall stands normal to 𝑥or to 𝑦, so its drilling rotation is 𝜃𝑥or 𝜃𝑦and its own pairing lands inside the 𝜃𝑥, 𝜃𝑦part. That part is therefore not a bending-bending term but the ribs’ share of the same block, which is where the block does its work, three quarters of its quadratic form on the critical mode falling on wall elements away from the rib feet.

Norm and effect are not interchangeable here, and the direct test is to exchange the blocks: solving the program’s exported classic pencil with its rotation-rotation block replaced by ours on the same rows returns 9.6685 on the rebuilt design A against the program’s 9.6413 and 4.8320 on the conforming panel against 4.8342, so a block 1.9 times the program’s in norm on its rib part reproduces the program’s load factor to 0.3% and 0.05% inside the program’s own stiffness. The reverse exchange is not informative: the program’s block dropped into our pencil meets a drilling freedom held by a weak spring where the program holds it by a stiff penalty, and the indefinite block then opens a spurious mode, as it does under a normal load in Section 5.3.

The two matrices thus act on the critical mode alike, their quadratic forms on the block form’s critical mode of the conforming panel standing in the ratio 0.92, and their worths agree to half a point; The load factor is sensitive enough

Element-Dependent Buckling Of Stiffened Panels

for that agreement to mean something, scaling the program’s own part moving its worth linearly at 2.7 points per unit of scale. What the two blocks share is therefore a structure and an effect, not a matrix. Norms are not load factors, so each exported pair is solved again. Rebuilt from its matrices, the classic pass of the conforming panel returns 4.900236 against the 4.900235 printed and the perturbation pass 4.767024 against 4.767024, with residuals below 9 × 10−11, so the exports are the pairs the program solved. The classic pair is then solved with its rotation-rotation block set to zero, Table 4. The block accounts for 99.86% to 100.14% of the difference between the passes on all four tied models and leaves at most 0.020% of the classic value, and on the shared-node rebuilds of the conforming panel and design A it leaves 0.0006% and 0.020%, identically on both releases. The other two differences between the pencils are two orders smaller: giving the classic pencil the perturbation pass’s translational blocks moves its load factor by at most 0.028% and its stiffness by at most 0.24%, and substituting all three differences returns the perturbation multiplier to seven figures on every model, so nothing else separates the two operators.

Two controls tie the difference to the operators rather than to the nonlinear base state. Taken with the base load applied forward from a tenth to three times the reference, Fig. 6, the perturbation pass of design A predicts between 8.41 and 8.60, a spread of 2.2% over a thirtyfold range, and 8.75 and 9.56 at four and five times, where the in-plane bending of the cantilevered panel makes the state visibly nonlinear; the rungs above, to 8.6 times, stay far from criticality, the multiplier still 3.55 at the last, with a change of the lowest mode near seven times. The ladder runs on the forward branch, while design A’s critical value lies on the reversed one, 0.3% lower at small base load; the reversed branch is carried along its own nonlinear path in Section 5.2. Under axial compression, where the block has nothing to weight, the classic pass returns 10.7018 and the perturbation pass 10.7018 to 10.7046 over the same ladder. On shear designs A and B the gap is 14.9% and 14.0%, on the combined design 12.8% and on the conforming panel 2.7%, and it survives refinement, 14.3% on design A at twice the rib subdivision and 13.4% on its shared-node rebuild at two subdivisions.

5.2. Which stress stiffness each formulation assembles Table 5 gives the load factors of the three forms on the tied designs with the two couplings entered one at a time, and Table 6 sets the three forms beside every commercial formulation, on the tied models of the optimization mesh and on the shared-node rebuilds. The moment-weighted coupling is inactive: on the conforming panel its quadratic form on the critical mode vanishes to five decimals and it moves the load factor from 4.9081 to 4.9072. The shear-weighted coupling is not: acting on the row that already carries the moment coupling, it takes shear design A from 10.2526 to 7.1328 and the conforming panel from 4.9072 to 4.6494, a quarter to a third of the load factor on the three designs against the row it acts on. The control variant of Section 2.1, transverse shear at one central point, moves these load factors by at most 1.1% and the worth of the block from 17.9%, 16.7% and 14.8% on the two shear designs and the combined design to 17.5%, 16.3% and 14.8%.

Each SHELL181 pass sits on one truncation. The classic pass lies below the block form by 1.65%, 1.63% and 1.65% on the three designs and 0.16% on the conforming panel, and the perturbation pass below the membrane form by 1.27%, 1.30%, 1.51% and 0.20%; both are of the size and sign of the 1.97% by which our compatible membrane stands above the program on the axial design, where every form and both passes coincide, and the remaining 0.3 to 0.7 points are not accounted for.

The correspondence does not depend on that offset: under the second variant of Section 2.1, whose membrane is the program’s, design A returns 8.6002, 10.1464 and 7.0363 for the three forms, and the classic pass then lies 0.57% below the block form and the perturbation pass 0.12% below the membrane form, the block’s worth reading 17.98% against 17.90% under the element of the results and 17.44% in the program. On the shared-node rebuild of design A the classic pass lies 1.24%, 0.59% and 0.79% below the block form and the perturbation pass 0.99%, 0.66% and 1.10% below the membrane form at one, two and three subdivisions, Table 7. The block’s worth agrees with the program’s to half a point on the four tied models, 17.90%, 16.71%, 14.84% and 2.76% here against 17.44%, 16.31%, 14.68% and 2.79% there, and moves by half a point when the drilling spring is replaced by the Hughes–Brezzi penalty over four decades of its weight.

Neither SHELL181 pass carries the couplings, and the measurement takes three quantities because the norm alone would say the opposite. Its exported translation-rotation block is not empty, 1.7166 × 104 on the rebuilt design A against SHELL281’s 3.2628 × 104 on the same mesh, and in both programs that block sits entirely on entries pairing a translation with a bending rotation, none of its norm on a drilling freedom, the drilling freedom of each node being read in its own plate’s frame, 𝜃𝑧on the skin and the in-plane rotation on a rib wall.

What separates them is the action on the critical mode. On design A the quadratic form of that block is −31% of the mode’s on SHELL281 and +0.03 on SHELL181, against +35% for the two couplings assembled here on their own

Table 3

Stress stiffness of the two passes, block by block: Frobenius norms by pair of freedom types, each off-diagonal block counted once, and the load factor of each pass. A norm here says what is present, not what acts; Section 5.2 measures the action of the translation-rotation block on the critical mode.

Table 4

The classic pair of SHELL181 as exported, again with its rotation-rotation block set to zero, and the perturbation pass of the same model.

−0.0204%

share = (𝜆classic −𝜆removed)∕(𝜆classic −𝜆pert), residual = (𝜆removed −𝜆pert)∕𝜆classic, on the classic critical branch.

Table 5

The assembled forms of the stress stiffness and the critical load factor each returns on the tied designs, the two couplings entered one at a time.

+ The 𝑄term Of Eq. (7)

mode, the sign following each operator’s own branch; on the conforming panel the three read −5.0%, +0.006% and −4.7%. Zeroing the block accordingly moves the rebuilt design A by 34.5% in SHELL281 and by 0.03% in SHELL181, and the conforming panel by 5.1% and 0.01%. SHELL181’s block is therefore present and inert on these structures: its norm ratio ‖𝐁𝝋‖∕‖𝐊g𝝋‖ is 32%, so it is not small, but 𝐁𝝋stands nearly orthogonal to 𝝋and does almost no work on the mode, where SHELL281’s block and ours do. What that block contains, if not the couplings of Eq. (7), is not established here.

SHELL281 carries both the block and the couplings. Its two exported stress stiffnesses share their rotation-rotation block to 1.2 × 10−3 on design A and 8 × 10−4 on the conforming panel, so its passes cannot differ by it, and they do not; the norm of the difference of their translation-translation blocks is 9% of either without moving the load factor by a part in ten thousand. Zeroing its rotation-rotation block lowers the rebuilt design A by 14.2% of its value and the conforming panel by 2.4%, the gaps between SHELL181’s two passes on the same meshes being 14.0% and 2.4%; zeroing its translation-rotation block instead raises them by 34.5% and 5.1%, where our block form lies 38.9% and

Table 6

The three forms against the commercial formulations and the uncalibrated 20-node continuum (SOLID186 and C3D20). Of the formulations tabulated, every pass other than SHELL181’s equals its classic one to 0.06%.

–

shared-node rebuilds, one subdivision; continuum on its own meshes

Axial

† forms under the second element variant of Section 2.1, incompatible modes and Hughes–Brezzi penalty, the conforming panel being run under that variant only; ‡ coarse continuum mesh. The continuum entries of design A and the axial design are the finest of three meshes, Table 7, equal in the two programs to 0.001%; a dash marks a case not run.

4.8% above our full form on the same meshes; zeroing both leaves 7.7315 and 4.7130, against 8.2934 and 4.7185 for SHELL181’s perturbation pass and 8.3763 and 4.7840 for our membrane form. The operator of SHELL281 therefore has the structure of the full form, measured block by block. Converged at 6.5983 and 4.5896, it lies 2.3% and 1.5% below our full form, which on design A settles at 6.7531 by its third subdivision and on the conforming panel was run to its second, where SHELL281 reads the same 1.5% below it, its third subdivision moving it by 0.007%.

Table 7 carries every quantity of the two rebuilt panels by mesh level, so that each percentage of this section can be read at its own level and against a stated base. The first step of level moves the load factors themselves by up to 6.3% on design A and leaves what is compared between them: the gap between the SHELL181 passes reads 14.0% at one subdivision, 13.4% at two and 13.4% at three, the block’s worth 16.55%, 15.43% and 15.05% over the three, S4 stands 1.2%, 1.0% and 1.7% below the full form on design A over the three and 1.4% on the conforming panel at two, and SHELL281 2.3% and 1.5% below it at its finest.

The second step settles everything on design A: the full form moves from 6.7518% to 6.7531%, 0.02%, the membrane form by 0.5%, the block form by 0.8%, the two SHELL181 passes by 1.0% and 0.9%, S4 by 0.7%, and the eight-node elements by 0.04% and less for SHELL281 and 0.06% for S8R, whose sequence is not monotone; on the conforming panel the four-node values stop at two subdivisions, where they had moved by 0.2% to 0.3%.

Abaqus S4 behaves as the full form, which is inferred from its results, its stress stiffness not having been exported. The yardstick of that inference is our full form, and its stiffness can be checked on the very modes at issue against the exported SHELL181 stiffness on the same nodes: on the translations of the full form’s critical mode SHELL181’s stiffness stores 0.2% more elastic energy than ours on both the rebuilt design A and the conforming panel, and with the bending rotations included and only the rotation about each node’s own plate normal left out, 1.7% more on design A and 7% on the conforming panel, the membrane form’s mode giving the same figures; taken on all six freedoms the ratio is 14 and 21, the program’s drilling penalty acting on drilling components our spring-held element does not set the same way, the effect that also collapses the six-freedom correlation below, and at a junction node the same rotation is a bending rotation of one plate and the drilling rotation of the other, so the rotational part of the comparison cannot be freed of the penalties entirely. What can be read is that on the translations the two stiffnesses agree to 0.2% on the shear-critical modes, and on the rotations to the size of the membrane offset of Section 4.2 on design A, nowhere near the 13% to 19% of the normal-load case in Section 6.2.

Both of its passes return the same load factor, which lies below the full form by 0.28%, 1.53%, 1.52% and 1.58% on the tied conforming panel and the three designs, and on the rebuilds at two subdivisions by 1.0% and 1.4%. The values for shear design B and the combined design were run after the other two, against a band written into the batch

Table 7

The shared-node rebuilds by mesh level: one, two and three subdivisions of each design cell, the four-node and eight-node meshes sharing a node lattice at every level. The continuum runs on its own three meshes, extrapolating to 6.97 on design A at observed order 3.2 and to 10.55 on the axial design at 2.25. A dash marks a level not run.

–

file beforehand from the earlier offsets, 7.14 to 7.23 and 4.99 to 5.05; they came out marginally below it, at 7.1357 and 4.9882. Their critical modes say the same. With the modal assurance criterion taken unweighted over the translations of every node, the Abaqus mode matches the full form’s critical mode to 1.0000 on all five tied models, where the membrane form’s reaches 0.956, 0.953 and 0.997 on the three designs and the block form’s 0.574, 0.584 and 0.929; on the tied conforming panel the membrane and block forms share one mode to 0.9999, so the mode separates the forms on the designs and not there. Over all six freedoms the criterion falls to 0.16 to 0.38 on the four designs, and returns to 1.0000 once every node’s rotation about the normal of its plate is left out, a rotation each program fixes by its own penalty. Figure 7 shows what those numbers look like on design A: under the full form and under S4 the loaded end lifts at one corner, the same corner with the same skin pattern, while under the block form, whose stress stiffness the block stiffens on the wall layers with nothing to offset it, the loaded end twists, the two free corners moving opposite ways.

One reference and one check place these results, a continuum that assembles no shell stress stiffness at all and a nonlinear path that assembles each element’s own. The first is the structure itself in 20-node hexahedra: skin 10 mm thick about the shell mid-surface and ribs 8 mm thick standing on its top face up to the blade top at 0.100 m, ribs and skin sharing nodes, nothing adjusted or calibrated, run as SOLID186 in ANSYS and C3D20 in Abaqus. Under axial compression it returns 10.6348, 10.5811 and 10.5639 on its three meshes, a Richardson value of 10.55, against the shells’ 10.5137 to 10.534 at two subdivisions, Table 7, where the four formulations agree with one another to 0.2%; where no truncation matters the continuum and the shells therefore describe one structure to 0.1% to 0.5%, the span from its extrapolated value to its finest mesh against the four shells at two subdivisions, and its perturbation pass equals its classic pass to 0.02%.

The junction is the one place where the continuum is not the structure the shells idealize. A shell blade runs from the skin mid-surface to 0.100 m, so its lowest 5 mm lies inside the skin, where the continuum stands its rib on the skin’s face and has 95 mm of steel above a solid skin. Raising the rib top by half a skin thickness, so that the blade above the face is the shells’ full 100 mm, brackets that idealization from the other side and moves the continuum by 13.4% on the axial design, 12.7% on design A coarse and 12.6% on design A medium, the two programs agreeing to 0.013% on each.

Table 8 gives the values. The shift belongs to the geometry and not to the load case, and the axial design settles which of the two models is the shells’ own: as built the continuum stands 0.1% to 0.5% above SHELL181, S4, SHELL281 and S8R there, each at its finest mesh, and with the rib raised it stands 14% above them. At 2.7% per millimetre of rib height that agreement fixes the junction geometry to within a fifth of a millimetre, and the bias it can carry into the shear comparison is the same 0.1% to 0.5%, against the offsets of 3% to 6% below.

On the rebuilt design A three meshes return 7.3062, 7.0606 and 7.0058, the two programs agreeing to 0.005% on the second and 0.001% on the third, a monotone sequence whose last step is 0.8% and whose Richardson value is 6.97 at the observed order 3.2, the mesh parameter being the in-plane cell between rib faces, the through-thickness count not following it: 3% to 6% above the full form, S4 and SHELL281, 11% below SHELL181’s perturbation pass, 23%

Authors:

Peder EZ Larson 1, 2,* , Jenna ML Bernard1, James A Bankson 3, Nikolaj Bøgh 4, Robert A Bok1, Albert P. Chen 5, Charles H Cunningham 6,7, Jeremy Gordon1, Jan-Bernd Hövener 8, Christoffer Laustsen 4, Dirk Mayer 9,10, Mary A McLean11 12, Franz Schilling13, James Slater1, Jean-Luc Vanderheyden5, 14, Cornelius von Morze 15, Daniel B Vigneron1, 2, Duan Xu1, 2, and the HP 13C

94143, Usa.

Denmark. 5 GE Healthcare, Menlo Park, California, USA. 6 Physical Sciences, Sunnybrook Research Institute, Toronto, Ontario, Canada.

8 Section Biomedical Imaging, Molecular Imaging North Competence Center (MOIN CC), Medicine, Baltimore, MD, USA. Cambridge, United Kingdom.

14Jlvmi Consulting Llc, Dousman, Wi, Usa

#See Acknowledgements for a list of all HP 13C MRI Consensus Group Members This work was supported by the ISMRM Hyperpolarized Media MR Study Group, the ISMRM Hyperpolarization Methods & Equipment Study Group, and the Hyperpolarized MRI Technology Resource Center (NIH/NIBIB grant P41EB013598).

Abstract

MRI with hyperpolarized (HP) 13C agents, also known as HP 13C MRI, can measure processes such as localized metabolism that is altered in numerous cancers, liver, heart, kidney diseases, and more. It has been translated into human studies during the past 10 years, with recent rapid growth in studies largely based on increasing availability of hyperpolarized agent preparation methods suitable for use in humans. This paper aims to capture the current successful practices for HP MRI human studies with [1-13C]pyruvate - by far the most commonly used agent, which sits at a key metabolic junction in glycolysis. The paper is divided into four major topic areas: (1) HP 13C-pyruvate preparation, (2) MRI system setup and calibrations, (3) data acquisition and image reconstruction, and (4) data analysis and quantification. In each area, we identified the key components for a successful study, summarized both published studies and current practices, and discuss evidence gaps, strengths, and limitations. This paper is the output of the “HP 13C MRI Consensus Group” as well as the ISMRM Hyperpolarized Media MR and Hyperpolarized Methods & Equipment study groups. It further aims to provide a comprehensive reference for future consensus building as the field continues to advance human studies with this metabolic imaging modality.

Keywords: Hyperpolarized MRI, metabolic imaging, carbon-13, pyruvate, dissolution dynamic

Introduction

MRI with hyperpolarized 13C agents, also known as hyperpolarized (HP) 13C MRI, has shown great potential as a novel imaging modality, particularly for its ability to probe metabolic processes in real time. The first human studies with HP [1-13C]pyruvate were performed in 2011 in prostate cancer patients (1).

Since then, there have been over 60 papers published with imaging results of human subjects from 13 different sites, with applications including prostate cancer, brain tumors, breast cancer, kidney cancer, pancreatic cancer, metastatic disease, liver disease, ischemic heart disease, diabetes and cardiomyopathies. The vast majority of these studies used [1-13C]pyruvate (1–63), where [2-13C]pyruvate (64) and 13C-urea (56) have been demonstrated too.

As clinical HP 13C MRI advances, there is a growing need to build consensus for best practices, which are critical for comparing data across sites, performing multi-site trials,deploying methods to new sites, partnering with vendors, and potentially for obtaining broader regulatory approvals.

In March 2022, we initiated an effort to build consensus within the HP 13C MRI community with this opportunity in mind, and it was greeted with strong enthusiasm. The “HP 13C MRI Consensus Group”, containing over 55 members from 27 sites, identified the area of greatest need and opportunity for consensus building to be HP [1-13C]pyruvate human

●

Pyruvate is the most mature and widely used HP agent and has the most significant translational evidence emphasizing the potential clinical impact.

●

Clinical trials, particularly multi-site trials, have the strongest need for consensus methods to ensure that data can be combined across sites. This work is a Position Paper for which the goal is to describe current successful practices and study methods for HP [1-13C]pyruvate human studies along with justification to support those practices. This is divided into four major topic areas: (1) HP 13C-pyruvate preparation, (2) MRI system setup and calibrations, (3) data acquisition and image reconstruction, and (4) data analysis and quantification (Fig. 1). The current successful practices and study methods include a literature review of published peer-reviewed journal papers showing human HP [1-13C]pyruvate study data, up to September 2022 (1–63), as well as new unpublished information from surveys of HP 13C study sites. Based on this information, we also highlight the evidence gaps, strengths, and limitations of current practices which are summarized at the end of each section.

Figure 1: Illustration of the HP 13C MRI human study process, including the 4 major areas covered in this paper: Hyperpolarized 13C-pyruvate preparation, MRI system setup and calibration, Acquisition and Reconstruction, and Data Analysis and Quantification.

Figure 2: Anatomical targets of HP [1-13C]pyruvate MRI human studies published up to September 2022.

Hyperpolarized 13C-Pyruvate Preparation

This section covers the processes for creating the HP agent, 13C pyruvate, and will include many aspects and considerations that are needed to safely and effectively prepare doses for metabolic imaging studies in human subjects. These include material, personnel, equipment and facility, fluid path preparation, quality control, and release.

It is helpful to understand that the specifications of a dose of 13C pyruvate suitable for in vivo MR HP metabolic imaging were shaped in part by early preclinical studies performed by GE HealthCare summarized in Ref. (65). In short, the safety of the two novel drug components, 13C pyruvate and the electron paramagnetic agent (EPA) AH111501, were demonstrated in those studies. The more precise formulation of the dose suitable for human use was then determined from clinical studies (66) that included two Phase 1 clinical trials in young and elderly healthy volunteers without hyperpolarization of the 13C nuclei and another Phase 1/2a dose escalation and imaging feasibility study with HP 13C pyruvate in 31 prostate cancer patients at the With the exception of the first HP 13C imaging clinical trial, which utilized a prototype device in a cleanroom (1), all HP 13C studies performed in humans to date have utilized the SPINlab polarizer (manufactured by GE HealthCare). Consequently all doses of the HP 13C pyruvate delivered by SPINlab have been produced using the “SPINlab Pharmacy Kit” that serves as the container-closure system for the various drug components (13C pyruvic acid and EPA mixture, dissolution medium, and neutralization and dilution medium) during sample polarization, dissolution and quality control (QC) processes. Thus many aspects of the HP sample preparation considerations discussed below are related to the SPINlab instrument and the consumables designed to be used with it (67).

General Considerations

While more than 860 patients or healthy subjects having been injected with HP 13C pyruvate as of January 2022 without reports of any serious adverse events (68), HP 13C pyruvate injection remains an investigational MR contrast agent and can only be administered by those with Investigational New Drug (IND) exemption from the Food and Drug Administration (FDA) in the USA, a Clinical Trial Application (CTA) in Canada, approval from National Research Ethics Committee Services in the UK, or approval from the relevant local regulatory body. Thus, methods and processes involved to produce a dose should have patient safety as the first priority. Since utilizing dissolution dynamic nuclear polarization (dissolution-DNP) for human use is still a relatively new development, there are no existing published regulatory guidelines specifically for this method.

There are two major production styles that determine how various sites approach the agent preparation. In the US, the most common approach is to rely on a sterilizing filter (“Terminal Sterilization”) to ensure sterility of the final product, akin to PET tracer production, where a starting molecule with a radioisotope is processed using various other ingredients to make the final, desired and injectable contrast agent within a necessarily short amount of time (69). For these sites, sterilization of the components and accessories upstream of this filter are not required, although many of them were manufactured and tested following Good Manufacturing Practice (GMP) or Good Laboratory Practice (GLP) requirements. The filling process is usually performed under an ISO 5 laminar flow hood, but a clean room or an isolator is not required.

This approach is typically accompanied by testing the integrity of the sterilizing filter prior to release of the dose for injection. Typically, post release endotoxin and sterility tests are performed using an aliquot reserved from each released dose.

In the UK and EU, the most common approach is to more-closely follow sterile pharmaceutical compounding guidelines (70), where all components and ingredients are required to be sterile or manufactured under GMP guidelines and are assembled and filled within a clean room environment or an isolator system (“Sterile Preparation”). Typically a batch of Pharmacy Kits for HP 13C pyruvate injection are prepared together. The sterility of the final dose is also ensured by batch validation testing, in addition to the sterility of the ingredients and the sterile compounding process. The endotoxin and sterility testing are performed for the process validation but are not performed for each injected dose.

Some institutions fill and assemble the Pharmacy Kit required for a specific study on the same day or the day prior to polarization, dissolution, and patient administration, but others have also demonstrated the feasibility of preparing a batch of kits, keeping them in a -20ºC freezer and using them over a period of a few months.

Beyond the obvious requirements that the process and the facility has to ultimately produce a dose that is safe to inject into a human, regulatory authorities will also focus on the question “Are you in control of your processes?”. To be in control of your process requires an in-depth and broad understanding of all processes involved in pre, post, and during the production process.

Personnel

It is typical and may be required to have licensed personnel involved in the production process depending on local regulations.Typically a pharmacist, radiopharmacist or other similarly qualified person (QP), in charge of the facility where the Pharmacy Kit filling and preparation is taking place, is responsible for the overall process and the release of the injectable dose.

Qualified cleanroom technicians are often involved in the Pharmacy Kit filling under the supervision of the pharmacist or QP. As is required for pharmaceutical compounding or PET tracer production, training requirements and training records for all personnel need to be maintained and available for audit by the FDA or equivalent.

Equipment And Facility

The facility and all equipment need to have standard operating procedures (SOPs) that describe how equipment is used, maintained, and calibrated to comply with relevant legislation. Currently, almost all the filling of the Pharmacy Kit takes place within a compounding laminar flow hood or isolator (typically ISO 5). At some sites, the filling is conducted within a cleanroom, while at others, it is conducted in a dedicated non-cleanroom space, reflecting differences in cleanroom approach and specifications between regulators worldwide (71). Some equipment or facilities, such as the compounding hood or cleanroom, may require external certified laboratories for testing.

Material Handling

Material handling guidelines (69,70) require SOPs detailing a system to track all of the materials involved in the HP production process for a particular patient dose, similar to current good manufacturing practice (cGMP) requirements for material handling for drug compounding. This includes acceptance standards, storage conditions, amount used in the patient dose for each ingredient and materials used in the assembly of the fluid path and Pharmacy Kit. Currently some users choose to open and inspect and sometimes modify the Pharmacy Kits upon arrival, but some users keep them in the sealed packaging until they are required for dose preparation.

Pharmacy Kit Filling And Assembling

As required by an IND or its equivalent, the preparation of the doses of HP 13C agent are detailed in the Chemistry, Manufacturing, and Control (CMC) section of an applicable regulatory submission; an example of this has been made available (72). It describes the processes of filling the Pharmacy Kit with the different components that make up the final drug product, and of assembling the final kit for either storage or immediate use in the polarizer. Special attention should be given to the laser welding process in order to satisfy installation qualification (IQ) and operational qualification (OQ). Typically, the final developed process is validated by process qualification (PQ) runs, during which 3 or more Pharmacy Kits are filled and used and the final HP 13C products are tested for endotoxin and sterility and to confirm that they meet the dose specifications for injections (usually including pyruvate concentration, residual EPA concentration, pH, liquid state polarization level and dose temperature). The data from 3 consecutive PQ runs are submitted as part of the IND submission (or its equivalent), and are often also reviewed by the Institutional Review Board (IRB) where the studies are conducted.

Quality Control And Dose Release

The quality control (QC) and dose release can be separated into two aspects: one is the QC and release of the filled Pharmacy Kit, and second is the QC and release of the HP 13C agent for injection, after polarization and dissolution. For institutions filling a batch of kits and storing them to use over a period of time, typically the batch can be released based on initial validation, environmental monitoring data from the day of kit production, and if filters are used during preparation of any of the components, filter integrity testing. But in some cases one or more kits are used for validation before the batch of kits are released for future use. For institutions that fill only the kits required for specific studies shortly before the experiment, the filled kits often do not go through separate release tests before they are used.

The quality control of the HP 13C pyruvate solution post dissolution is primarily performed to ensure that the agent meets the dose specifications (Table 1) before it is administered to the subject. These specifications target both safety (pH, residual EPA, temperature) and efficacy (pyruvate concentration, polarization, volume). Typically, the pyruvate concentration, residual EPA concentration, pH, dose temperature, dose volume, and liquid state polarization are measured by the QC accessory associated with the SPINlab polarizer. Some users perform a secondary measurement for one of the parameters, such as pH, using a different instrument or pH paper. For sites that do not go through a separate release testing process for batch filled kits, the integrity of the sterilization assurance filter, a part of the Pharmacy Kit, is typically tested as a part of the dose release. It is also common for these users to preserve an aliquot of the final HP 13C pyruvate solution for post-release endotoxin and sterility testing. This testing cannot be completed fast enough to test an individual dose prior to injection, but this is why other processes such as PQ runs and validation testing are done to minimize the chance a subject could be injected with a contaminated dose.

The Final Dose Release And Injection

should be done under the supervision of a licensed professional, based on local regulations.

Some Key Challenges

Many of the challenges associated with HP 13C pyruvate preparation can be attributed to the conditions required for the dissolution-DNP method of high magnetic field (~3-7 T) and very low temperature (~1 K) during polarization, with pressurized and superheated water necessary for the rapid dissolution event. These extreme conditions are quite challenging for the design of the container-closure and fluid path system. In particular, the cryogenic temperature in the polarizer requires special attention to any moisture or ambient (moist) air introduced into that portion of the fluid path, which can form an ice block at ~1 K. This ice can lead to flow restriction during the dissolution event and reduce the strength of the laser welded bond between the cryovial and its cap. This can ultimately produce failures in the dissolution step, including variations in final pyruvate concentration and pH that may fail to meet QC release criteria as well as fluid path ruptures that provide no available dose and result in polarizer down-time.

The polarization of the HP 13C pyruvate sample decays quickly over the span of a few minutes after dissolution, and thus the process of dissolution, QC for release, and injection should be completed as fast as possible to preserve the high polarization level achieved. Any delays in the preparation process, such as transportation time or equipment malfunction, can significantly reduce the final polarization and result in lower quality imaging data.

Current Practices

A summary of data collected from all sites performing clinical trials with HP 13C-pyruvate is shown in Fig. 3 and Table 1, including the specification of the final dose and how the quality control and release of the final dose are performed. There is a split in the Production Style, described in the General Considerations section above, with 8/13 sites using Sterile Preparation versus 5/13 using Terminal Sterilization. While many of the dose specifications show notable differences in acceptable ranges, all of these variations listed in tables have been successfully and safely been used to perform HP 13C pyruvate studies in humans. Their differences depend on the institutions’ preferences, resources and their particular regulatory situation. There is high similarity in pyruvate ranges, temperature ranges, EPA limits, and volume limits. There is modest variability in pH ranges and large variability in the endotoxin test limit. There is a 3-fold difference in acceptable polarization levels, which are measured to ensure a futile dose is not injected since the polarization is directly proportional to SNR. This reflects the decision by several sites to believe that useful data can be still be obtained with suboptimal polarizations.

Figure 3: Hyperpolarized agent preparation methods reported by sites currently performing HP

In House

Table 1: HP 13C-pyruvate preparation parameters, methods, and dose specifications used for quality control testing and release as well as validation. These were obtained from a survey of all sites performing clinical trials with HP [1-13C]pyruvate. The parameters used for product release are noted in bold text, otherwise these parameters are measured for batch validation or other QC measurements. The endotoxin and sterility testing are performed during process validation of the batch and/or post-injection, and largely depends on the agent production approach.

Summary

The overall safety record of HP 13C-pyruvate has been very strong, and the SPINlab hyperpolarizer has proven to provide high polarizations at human sized doses while meeting numerous QC and release criteria. A weakness remains the failure modes of the SPINlab Phamacy Kits (e.g. ice blocks, path ruptures), which are placed under extreme requirements particularly during dissolution. The preparation process still requires a high degree of expertise.

Therefore, there is a significant need to improve the reliability, robustness, and ease of operation for generating HP 13C-pyruvate doses for human studies. Furthermore, there is a divide between manufacturing and sterile compounding style preparation as well as other site-specific practices, resulting in variations in SOPs and justification required to relevant regulatory bodies. There have also been no comparisons between these approaches. It is also unclear what release criteria and QC parameters are truly required to ensure patient safety.

However, all of the reported methods are acceptable and approved by the appropriate regulatory authorities, and have led to the rapid expansion of successful human studies in recent years.

Mri System Setup And Calibrations

This section covers the MRI system setup, including the imaging system, RF coils, phantoms, and prescan calibration methods.

Imaging System

The main prerequisite for a given MRI scanner to be capable of supporting studies with HP 13C is its “broadband” capability to transmit and receive radiofrequency (RF) signal at the frequency of 13C, which is around 4 times lower than 1H. This does not come as a default on clinical MR devices. The transmit power of the broadband amplifier should also be sufficient to support the intended flip angle and RF pulse shape with the employed transmission RF coil(s) for 13C. Most studies to date use relatively low flip angles (< 90 degrees) for HP 13C in order to preserve polarization for time-resolved imaging. The capability to receive 13C signal on multiple channels is also desirable to increase SNR, as discussed further in the “RF coils” section.

The choice of magnetic field strength is primarily dependent on the metabolites’ frequency separation due to chemical shift dispersion and 1H imaging. High field strengths do not enhance hyperpolarized 13C signal as they do for 1H because the signal strength in a HP experiment relies on manipulating the population of quantum energy states outside of the MRI scanner.

However, the injected HP 13C-pyruvate and its metabolic products have greater frequency separation at higher fields, and it may thus be easier to separate and quantify these resonances at higher fields. This comes at the cost of a reduction in the achievable T2* and often reduced T1. As the initial polarization is independent of the imaging field strength it has been proposed that the increased T2* at 1.5T can potentially be exploited to increase SNR by adapting the acquisition bandwidth or reduce off-resonance imaging effects in cases when the decay of the transverse magnetization is dominated by T2* (73). In practice, 3T has been used in all published human 13C-pyruvate studies surveyed (Supporting Table S1), and comprises the majority of scanners currently in use for human studies (Table 3). A field strength of 3T is well-suited for 1H MRI anatomical reference and correlative imaging.

Stronger and more rapidly slewing magnetic field gradients support more rapid spatial encoding, particularly for metabolite-specific single-shot imaging using echo-planar imaging (EPI) or spiral imaging (See “Acquisition and Reconstruction”). Although the spatial resolution acquired for HP 13C imaging is typically much coarser than for 1H MRI, the factor of ~4 in gyromagnetic ratio leads to the same reduction factor in performance of the gradient system, so 13C experiments are potentially more limited by gradient hardware performance. To date, all human studies have used the commercially-available integrated gradient systems provided in clinical MRI scanners.

Optimization of scanner design has understandably focused on minimization of artifacts in 1H MRI, where devices such as room lights, the gradient amplifiers, and the motors driving the patient bed are checked to ensure that they do not produce RF interference at the 1H frequency, but artifacts may arise at other frequencies. Eddy current compensation is also not always appropriately adjusted for nuclei at other frequencies (74). In order to optimize for 13C, many sites have performed checks on phantoms for RF interference, gradient artifacts, and eddy currents (74), including the use of post-hoc gradient impulse response function characterisation and correction, and some vendors have fixed these issues as well.

Rf Coils

For HP 13C imaging studies in humans, RF coils for both 1H and 13C nuclei are needed, with 1H MRI providing an anatomical reference for registration and optional additional multiparametric MRI readouts. At the Larmor frequency of 13C nuclei, the relative contributions from coil noise compared to sample noise increase compared to 1H (73,75), although sample noise still is likely the dominant contributor for human-sized coils at 32.1MHz - the resonance frequency of 13C nuclei at 3T.

The key requirement for human 13C-pyruvate RF coils are that the coil geometry and sensitive volume must cover the volume of interest in the subject. Table 2 and Figure 4 shows coil configurations that have been used and optimized for applications in different anatomic regions.

Volume resonators are most commonly used for transmit, as they surround the subject to

Provide B1 Transmit Across The Fov (B1

+). While 1H relies on a large birdcage (“body”) coil built into the scanner, 13C transmit coils must be placed inside the bore. This takes up valuable space within the magnet, and also has led to the use of designs with relatively inhomogeneous

B1

+. Many human studies have used Helmholz pair resonators for transmit, including the “clamshell coil”, which has a notably inhomogeneous B1

+ Profile But Has Been Used Because Of

relatively easy integration into the scanner bore. B1

+ Variation Results In Variations In The Flip

angles that control the use of the hyperpolarized magnetization and creates errors in common HP metrics (9,76). The exception are head coils, where birdcage designs with highly

Homogeneous B1

+ can be placed around the head while easily fitting inside the bore. As with 1H MRI, higher SNR can typically be achieved by smaller receive coil elements, such as surface coils or phased arrays, and the majority of 13C receive coils used have layouts similar to 1H phased arrays.

RF coil quality control is important to ensure proper functioning of the coils to provide consistent imaging quality, especially with limited natural abundance 13C signal in vivo. It typically involves 1) a physical integrity check of the coil cables and connectors and 2) phantom SNR tests to check the coil’s performance and to monitor it over time (see Phantoms below). An useful reference for RF coil quality control is outlined in the MRI accreditation program of the American College of Radiology (77) and can be adapted for 13C coils.

Notably, configurations for brain and prostate studies used dual-tuned 1H/13C coil designs, which greatly simplify workflow and registration of 1H and 13C images, as no switching of coils is needed.

(1)

Table 2: RF coil configurations reported for human HP [1-13C]pyruvate studies.

Tx = Transmit

coil, RX = receive coil. The commonly used “clamshell” TX coil is a Helmholz pair design. For 1H RF configurations, all used the Body coil for TX unless otherwise noted, and “repositioned” indicates the 13C coil was removed for 1H imaging. One representative reference is listed for each configuration. The RF coil configurations reported in the reviewed papers are shown in Supporting Table S1.

Figure 4: Examples of RF coil configurations used for human HP [1-13C]pyruvate brain studies. (A,B) 13C Clamshell TX (Helmholz pair) and 2× 4-channel paddle RX arrays. (C) 13C Birdcage volume TX and 32-channel RX array (RX array slides into TX coil). (D) 13C Birdcage volume TX and 24-channel RX array, combined with a 1H 8-channel RX array. Image reproduced with permission from Ref (16).

Phantoms

Since hyperpolarized magnetization is non-renewable, phantoms containing 13C nuclei are important to: 1) test the multi-nuclear capabilities of the imaging system, including all parts of the signal excitation and receive chain; 2) perform calibration measurements before a scan with hyperpolarized nuclei; and 3) perform necessary pre-scan adjustments (see “Prescan Calibration” section). The phantoms currently in use are listed in Table 3. Their composition must provide sufficient 13C signal, with additional considerations of conductivity, stability, chemical shift(s) present, potential for dynamic imaging, and cost. The phantom geometries are typically either compact, in order to be used alongside the subject during a HP scan, or large enough to mimic the inner volume of a RF coil for system testing.

One popular compact design contains enriched 13C-urea at high concentration, typically 8 M, which provides a single resonance, placed inside a small container ~1 mL. The most common recipe mixes 13C-urea in a 90% water/10% glycerol solution, with glycerol used to increase the urea solubility and doping with a Gd-based contrast agent to shorten T1 which increases the potential SNR per unit time. For example, when Dotarem is added at a 3:1000 volume ratio the 13C-urea T1 is around 500 ms and T2 is around 100 ms. However, when testing pulse sequences influenced by T1 and T2, doping should be used carefully. This phantom is suitable for frequency calibration, transmit gain calibration, sequence testing, and as a fiducial marker when placed next to a patient. However, enriched 13C-urea has a relatively high cost compared to natural abundance compounds.

For larger volumes (>100 ml), the phantoms most often used contain undiluted ethylene glycol, glycerol, or dimethyl silicone. These compounds have sufficiently high carbon concentrations to provide sufficient 13C signal even with the 1.1% natural abundance of 13C. These larger phantoms matching the inner volume of an RF coil are useful for coil testing, including transmit

+) And Receive (B1

-) coil profile mapping, as well as to mimic acquisitions using in vivo FOV requirements. In this case, size and conductivity should match the expected subject size in order to mimic coil loading and get a realistic estimation of B1+. Large-volume natural abundance urea phantoms have also been used by some sites, but suffer from higher conductivity compared to biological tissues. Typically, it is easier to increase the conductivity and hence coil loading of the non-conductive phantom by adding NaCl to match physiological loading (16,78).

Dynamic phantoms that aim to mimic metabolite kinetics have also been developed (79–81), and have the potential to more closely mimic the HP experiment, but so far these are not widely used.

Prescan Calibration

Prior to performing an MRI acquisition, the so-called prescan procedure is used to set the shim parameters to maximize B0 homogeneity over the field of view (FOV) or a specific region of interest (ROI), the scanner center frequency (CF), the RF transmit gain, and the receiver gain.

While this calibration procedure is usually automated for 1H, the lack of sufficient natural abundance 13C signal prevents use of automated methods. (Although natural abundance 13C lipid signal has been detected, there are so far no reports on using this signal for prescan.) Table 3 shows current practices across sites.

Maximizing B0 homogeneity is independent of the nucleus and is therefore performed prior to 13C imaging using the 1H water signal and existing shimming tools, such as by a standard automated process (“Auto Shimming”) or using high order shimming routines. Similarly, the 13C CF can be calculated from the 1H CF using a predetermined scaling factor that depends on the target chemical shift (82). Another common approach used is to have a small, high-concentration 13C phantom, e.g. 8M 13C-urea, integrated in the RF coil or placed next to the scan subject (1). The reference frequency can also be based on real-time measurements after the HP injection but prior to imaging (83). Both the CF and B0 shimming are critical when using spectrally-selective RF pulses, as inmetabolite-specific imaging methods, where the desired excitation bandwidths are typically very narrow and frequency offsets can lead to a failure mode that is only apparent after injection.

The calibration of the RF transmit power is typically performed on a small, high-concentration 13C phantom placed near the region of interest during the scan or on a large 13C phantom of similar size and coil loading as the subject, prior to the subject scan. Reference power is often done by sweeping the power in a pulse-acquire sequence (53,62), or the Bloch-Siegert method (52,84). When using a small phantom, the location of the phantom, B1

+ Inhomogeneity As Well

as any shielding effects, e.g., when the phantom is integrated into a coil (1), may degrade the accuracy. Other methods include real-time Bloch-Siegert method measurements after the HP injection (83), and using the stronger natural abundance 23Na signal that is close enough to the 13C resonance frequency to be detected by 13C coils (82).

The receiver gain is predetermined, either systematically based on independent phantom measurements and assuming the dose and polarization of the HP compound is known prior to injection, or based on past HP imaging studies.

Power [Kw]

Phantom(s) - during study Phantom(s) - before study 13C Frequency

8

13C-bicarbonate doped with dimethyl silicone, various

Power [Kw]

Phantom(s) - during study Phantom(s) - before study 13C Frequency

Maximum Values

Table 3: Summary of the imaging systems, phantoms, and prescan procedures used at sites currently performing HP 13C-pyruvate human studies. These were obtained from a survey of all sites performing clinical trials with HP [1-13C]pyruvate. *Previously performed studies with a Siemens 3T Tim Trio. The imaging systems, phantoms, and prescan procedures reported in the reviewed papers are shown in Supporting Table S1.

Summary

Commercially available 3T MRI systems are by far the most commonly used for human HP 13C-pyruvate studies, although a systematic investigation of the impact of B0 has only recently been investigated (73). The multi-nuclear RF transmit and receive chain has proven sufficient for current acquisition strategies, although many sites have observed artifacts due to RF interference, gradient interference, and residual eddy currents when operating at the 13C frequency. A variety of 13C RF coils, tailored for numerous anatomical targets, have been successfully demonstrated, with the main limitation that most transmit coils take up a lot of additional space inside the bore and provide relatively inhomogeneous B1

+ Profiles. The

phantoms used have converged into generally 2 categories - small phantoms containing 13C-enriched compounds that can be used during the study and human-sized phantoms containing compounds with high carbon concentrations but without 13C enrichment that are used to test and calibrate the coils. There are no standardized compositions or geometry, and dynamic phantoms that recapitulate in vivo kinetics would be desirable but are still an emerging area. Prescan calibration procedures were not well defined in most publications, so we surveyed individual sites to determine current practices. Calibration procedures for the B0 field (13C CF and shimming) for most sites take advantage of 1H signal and methods, while methods

For Calibration Of B1

+ is more variable across sites, likely a reflection of remaining challenges in how to perform this calibration. Standardization of both phantoms and calibration procedures would synergistically improve the robustness and reproducibility of HP 13C studies.

Acquisition And Reconstruction

Data acquisition strategies in human HP [1-13C]pyruvate MRI studies must account for multiple chemical shifts, efficiently utilize the non-renewable HP magnetization, and acquire data quickly relative to metabolism and relaxation decay processes. These studies require spectral encoding to separate metabolites, necessitating pulse sequences that efficiently encode up to 5D data (3 spatial + 1 spectral + 1 temporal dimension). RF pulses must efficiently sample without immediately saturating the non-renewable HP magnetization, and sequences must acquire data quickly and be robust to both experimental and physiologic variation (e.g. B1

+ Inhomogeneity,

variation in perfusion) to ensure reproducibility and minimize scan-to-scan variability. This section covers current successful practices for data acquisition in human [1-13C]pyruvate studies, and accompanying 1H imaging, from different anatomic regions, including scan parameters and image reconstruction.

Acquisition And Reconstruction Methods

The acquisition methods used in human [1-13C]pyruvate studies can be classified into 3 categories: 1) MR spectroscopy or MR spectroscopic imaging (“MRS/I”), 2) chemical shift encoding methods, and 3) metabolite-specific imaging (Fig. 5).

Mrs/I Methods Specifically

resolve a spectrum that can be analyzed to extract expected as well as unexpected resonances, making this approach very robust. It was used in many initial studies (1).

Chemical Shift

encoding methods, most commonly the Iterative Decomposition of water and fat with Echo Asymmetry and Least-squares estimation (IDEAL) method, use imaging sequences acquired with multiple TEs and rely on a model-based separation of expected chemical shifts (85).

Metabolite-specific imaging methods use specialized RF pulses that are spatially and spectrally selective to excite individual metabolites which are then typically imaged with fast k-space trajectories such as echo planar imaging (EPI) or spirals (86).

Their Application To Different

organ systems is described below. The image reconstruction methods used in human [1-13C]pyruvate studies have typically been conventional methods (e.g. FFT, non-uniform FFT, or equivalent). The incorporation of accelerated imaging and advanced reconstruction methods including parallel imaging (4,57,87) and compressed sensing (7) has also been applied in human studies for improved spatial resolution, temporal resolution and coverage, but have the potential for additional artifacts as well as SNR losses due to ill-conditioning of the reconstruction (e.g. g-factor).

The Majority Of

published studies do not use accelerated imaging indicating the resolution and coverage achievable without acceleration is currently adequate for successful data collection. Performing coil combination, even with fully sampled data has also been shown to have specific challenges for HP human images: using naive sum-of-squares methods suffer from high noise amplification in the relatively low SNR regime of HP [1-13C]pyruvate (compared to 1H), motivating several HP 13C-specific methods that include data-driven coil sensitivity estimation which have shown obvious improvements over sum-of-squares (11).

More recently denoising techniques have been applied as post-processing of human HP data(41,42,44). The techniques applied are based on spatial-temporal singular value decomposition for unsupervised estimation of signal and noise components. They have shown improvements in apparent SNR in the brain and liver, while care must be taken to choose parameters such as the rank threshold to avoid oversmoothing and overfitting to the estimated signal components.

Prostate Studies

Prostate cancer was the first human application of HP [1-13C]pyruvate (1), and data was acquired with MRS/I methods: 1D dynamic MRS, single-slice 2D dynamic echo-planar spectroscopic imaging (EPSI), and single time point 3D EPSI. Advances in imaging strategies led to the development and application of new acquisition schemes, including undersampled 3D EPSI with compressed-sensing (7), model-based chemical shift encoding methods that use a priori information (47,59), and metabolite-specific EPI (10), all of which can provide volumetric whole-organ coverage and dynamic acquisitions.

The pyruvate bolus arrival in the prostate can vary by ± 10 s between patients, necessitating dynamic imaging to reliably and consistently capture the pyruvate bolus (18). For this reason, all currently ongoing studies acquire dynamic data. While MRS/I, chemical shift encoding, and metabolite-specific imaging can all achieve dynamic imaging, chemical shift encoding and metabolite-specific imaging provide greater dynamic and volumetric coverage (85). For scan prescriptions, the FOV is designed to provide full prostate coverage and typically to match the orientation of the anatomic imaging used for registration. Flip angles used in current studies are constant through time, as quantification with a variable-through-time flip scheme is highly sensitive to bolus timing (8) and errors in the RF transmit (B1 +) field (76).

Heart Studies

Data acquisition methods for 13C imaging in the heart must be designed to meet the demands of significant cardiac motion and blood flow. To cope with the periodic cardiac motion, most human heart studies to date used gating to the diastolic window, the longest cardiac cycle interval, which has reduced motion (2,22,28,30,35,36,38,45,52). The duration of the diastolic window limits the available data sampling time, making cardiac acquisitions the most time-constrained of the HP 13C MRI applications. The most common acquisition approach is metabolite-specific imaging with spiral k-space trajectories (2). Their single-shot imaging capability makes these methods particularly robust to motion effects. Furthermore, spiral k-space trajectories provide rapid k-space coverage and relatively benign flow and motion artifacts. The majority of studies have used 2D multi-slice acquisitions, but 3D encoding has also been used successfully (35).

Brain Studies

For HP 13C MRI of the human brain, the majority of studies have also used 2D (slice selective) acquisitions (10–12,14,16,28,33,40,41,44,51,53,60), with a trend toward volumetric coverage using 2D multi-slice metabolite-specific imaging. 3D metabolite-specific imaging of the whole brain, with phase encoding of the slice direction (34,57), has been shown to provide similar SNR efficiency (88) compared with multislice imaging. A number of studies have employed MRS/I (5,6,29,31–33,50,55) resulting in a spectrum from each voxel, which has the advantage of not requiring a priori information about which peaks to encode. This was important in early brain studies when it was not known which peaks would be detectable. Chemical shift encoding, using a set of images with different echo times and an iterative reconstruction of the individual resonances (i.e. the IDEAL approach (85)), has also been used (12,49,54), with the drawback that coverage in the slice direction was limited due to the time required to acquire multiple echo time images.

Abdomen And Breast Studies

The fundamental approaches to data acquisition and reconstruction in the abdomen and breast are largely similar to the aforementioned applications, but demand attention to particular challenges associated with these anatomic regions, especially relating to respiratory motion.

Although it has been shown that a basic 2D MRSI approach based on phase encoding and FID readout can be successfully applied for HP 13C imaging in breast (15) and kidney (13), major advantages in terms of spatiotemporal resolution and coverage have been realized using tailored approaches based on metabolite-specific imaging (43,62) and chemical shift encoding (43), which have facilitated multi-slice or 3D dynamic acquisitions over large FOVs in the abdomen (4,37,46).

The significant respiratory motion encountered in these regions can directly blur 13C images, and has further favored these rapid acquisition strategies. Motion also degrades B0 homogeneity, which can shift frequency-selective excitation profiles and introduce artifacts into rapid imaging readouts. This makes accurate determination of the acquisition center frequency and shimming essential in these regions which often cover large FOVs. (See “Prescan Calibration” section for more information). In some studies, breath-holding was used to minimize motion effects and enforce frame-to-frame data consistency (42). A pragmatic and reasonably effective approach for dealing with respiratory motion during 13C data acquisition is an initial breath-hold (as long as can be tolerated), followed by free-breathing (46,62).

1H Imaging

Collection of 1H imaging data is essential both for prescribing the 13C acquisition and for interpretation of the resulting 13C data. Multi-planar 1H scouts are acquired prior to 13C acquisition to enable graphical prescription of the 13C imaging region. All human HP 13C-pyruvate imaging studies acquire conventional MRI scans (e.g. T1- and T2-weighted volumes) for anatomic reference, aiming to cover at least the full 13C FOV. Acquiring these anatomic scans as close as possible to the time of 13C imaging (immediately before or after) minimizes potential misregistration between the data sets. Depending on the application, other advanced 1H sequences are also acquired (e.g. diffusion-weighted imaging for cancer imaging).

When contrast-enhanced data is acquired, it is done after 13C imaging, as paramagnetic contrast agents will accelerate 13C relaxation.

Reported Study Parameters

Figures 5 and 6, and Supporting Table S2 shows the reported acquisition study parameters for human HP [1-13C]pyruvate studies published as of September 2022. Figure 5 shows a mixture of MRS/I, metabolite-specific imaging, and chemical shift encoding methods have been successfully used, where spectroscopy-based methods have become less prevalent in recent studies. Figure 6 shows the acquisition timing, including the important start time and interval/temporal resolution, is quite variable across studies.

Figure 5: Acquisition methods used in published HP [1-13C]pyruvate human studies published up to September 2022, classified into: MR spectroscopy and spectroscopy imaging (MRS/I); chemical shift encoding methods, such as IDEAL, that use multiple TEs and model-based reconstructions; and metabolite-specific imaging methods that use spectrally-selective excitation to image a single resonance at a time.

Figure 6: Temporal acquisition characteristics reported in HP [1-13C]pyruvate human studies published up to September 2022. (a) Reported referencing of acquisition start times.

(B)

Acquisition start times reported when using dynamic imaging and when timing was reported relative to the end of the injection. (c) Temporal resolutions. “Not Applicable” indicates dynamic imaging was not used.

Summary

Three general categories of acquisition strategies have been used successfully for human HP 13C-pyruvate studies: MRS/I, model-based chemical shift encoding (e.g. IDEAL) methods, and metabolite-specific imaging methods. These have enabled successful studies in the prostate, heart, brain, abdomen, and breast. Recent studies increasingly have used the imaging-based strategies of metabolite-specific imaging and chemical shift encoding which are the fastest methods, although a heads-to–head comparison between techniques has not been performed.

Metabolite-specific imaging is quite popular because of its speed and compatibility with single-shot imaging, but is sensitive to B0 field variations and thus requires careful calibrations. Nearly all studies surveyed acquired data dynamically, allowing measurement of the bolus and metabolite kinetics. The exact timings and associated flip angles vary quite widely across reported studies, with no consensus yet as to how to choose these parameters. Image reconstruction is typically done directly using Fourier Transform methods, and accelerated imaging strategies are uncommon.

Data Analysis And Quantification

This section covers the analysis of data from human HP [1-13C]pyruvate studies, including modeling and metrics, visualization, as well as considerations for how to store data and metadata. Depending on study design, the analysis may need to give quantitative or semi-quantitative output reflecting a biological process or may just reflect a contrast between different regions of interest for quantitative evaluation.

Metrics

Figure 7: HP [1-13C]pyruvate raw data (A) have typically been quantified using four categories of metrics depending on the acquisition. Data acquired as a single time point are often quantified using normalized metabolite images or metabolite ratios (B). Dynamic data can be quantified using normalized metabolite images or metabolite ratios (B), or with metabolite timings such as time-to-peak (TTP) or pharmacokinetic (PK) models (C). The latter two require the data to be time-resolved. [1-13C]alanine and 13C-bicarbonate are analyzed similarly to [1-13C]lactate but omitted here for display.

Metabolite images are commonly used as summary metrics for HP MRI data, often including some form of normalization as well as summed over time as an area under the time curve (AUC) (17). These are analogous to the visual evaluation that is most used for routine clinical work (89,90). In these metabolite images, we expect that the [1-13C]pyruvate AUC signal is predominantly weighted towards perfusion and uptake, while [1-13C]lactate, [1-13C]alanine and 13C-bicarbonate AUCs represent metabolic conversion. The strength of this approach lies in its simplicity and relatively few underlying assumptions. Limitations to the use of single-metabolite images or AUCs include sensitivity to inhomogeneous coil profiles (57,87,91), the acquisition strategy and acquisition parameters, pyruvate polarization and concentration level, and signal relaxation rates (92). Further, the reader must be careful to interpret all the images in conjunction to better understand the underlying biology; for example, increased [1-13C]lactate in the presence of decreased [1-13C]pyruvate delivery can have a very different meaning compared to increased [1-13C]lactate with increased [1-13C]pyruvate delivery.

In an attempt to address variations in coil sensitivity, polarization level, and pyruvate delivery, AUC images are often computed by normalizing to a specified parameter, such as the maximum pyruvate or average lactate signals, or presented as a ratio such as lactate/pyruvate or divided by “total Carbon” - the sum total of HP 13C signal observed across all metabolites. The AUC ratios between metabolites and pyruvate are proportional to the corresponding forward kinetic rates (81,93), but are not directly comparable to rate constants when magnetization loss rates (e.g. relaxation and losses due to signal excitation) differ between studies. Similarly, the ratios between the produced metabolites (e.g. bicarbonate/lactate) can reflect the balance between downstream metabolic pathways (12,55). Care must be taken to consider how AUC images are calculated and normalized before comparing values between studies.

To further quantify the interpretation, pharmacokinetic (PK) modeling approaches were developed to compute the apparent kinetics of pyruvate-to-metabolite exchange (92,94–99). These yield semi-quantitative to quantitative apparent rate constants, given in s-1. Some models require a vascular input function, while others avoid this requirement (95). PK models can explicitly account for acquisition-specific details such as excitation angle and repetition time, and thus may reduce the effects of these details on quantification. An input-less model, provided in the Hyperpolarized-MRI-Toolbox (https://github.com/LarsonLab/hyperpolarized-mri-toolbox) (100) and thus frequently employed for human data, has been shown to fit well and robustly to prostate and brain data (8,20). PK models are quantitative in nature, arguably provide more relevant biological information (8,20), and appear to be reproducible across sites (51). However, rate constants derived from PK models are still apparent rates, and likely do not reflect a single biological characteristic.

Some additional considerations include whether complex or magnitude data is used, as the noise behaviors will impact the analysis differently. Additionally, cut-off thresholds or other criteria may be used to identify and avoid voxels with insufficient SNR before analysis to improve robustness (20,41).

Regardless of the analysis approach, the underlying biology is not always clearly represented by the data; instead, the metrics may be influenced by perfusion, barrier permeability, intercellular shuttles, enzyme activities, co-substrate concentrations, or combinations thereof, depending on the organ and disease of interest (19,43,94,101–103). This may be addressed by incorporating complementary information. As an example, HP 13C pyruvate data is influenced by perfusion, and thus addition of perfusion MRI could be important for interpretation (98,104,105).

All the methods outlined above have been explored in clinical studies, described in Supporting Table 3 and summarized in Figure 8. As of September 2022, approximately 52% of studies involving human subjects report rate constants derived from a PK model with a few different models reported. A nearly equal fraction (51%) of the studies report AUC ratio values.

Approximately 66% of these studies report metabolite-specific images or AUC values. About 40% report SNR values; this metric is particularly frequent in manuscripts that describe technical developments for clinical HP MRI. Approximately 16% of these studies summarize model-free metrics, and 10% report measurements from a single timepoint. Most studies report a combination of quantities.

Figure 8: Reported metrics used for analysis in HP [1-13C]pyruvate human studies published up to September 2022.

Visualization

A wide variety of approaches have been used for visualizing data from human HP 13C-MRI studies. The challenges and practical considerations are: 1) choosing the appropriate metrics to display, 2) how to encode the parameters (e.g. the colormap), and 3) choosing how to provide anatomical context and other multi-parametric data. The choice of visualization also depends on the goal which could be for diagnostic interpretation, but also quality control, reproducibility among readers and publication.

Metrics

The choice of HP 13C metrics is described in detail above. At this stage in HP 13C development where there is no standardized metric, often a combination of metabolite images and ratios or PK model parameters are shown.

Parameter Encoding

The mapping function chosen should provide an adequate, often quantitative, impression of the parameter mapped. There is a consensus in the visualization field that perceptually uniform maps are best suited to visualize continuous parameters, like the greyscale typically used by radiologists as well as other monochrome (black to blue) and color ranges (fire-type, rainbow-type) (106,107). Multi-color heatmaps have been the most frequently employed method for HP 13C data, while greyscale has infrequently been used but it ensures there is no coloring-based bias as well as facilitating later reuse (Fig. 9a). Among the color schemes employed in the clinical HP 13C literature, fire-type scheme seems to be the most common [similar to “Plasma” or “Inferno” in matplotlib.org]. Next most commonly employed is the rainbow-type scheme [similar to “Rainbow” in matplotlib.org].

Anatomical Context

HP MRI faces the challenge that it does not necessarily depict the anatomical features, similar to PET, and thus requires an anatomical reference. Most often, a grayscale anatomical image is overlaid with a HP colormap (Fig. 9c,d). This approach is very intuitive, but can skew perception as the grey-scale anatomical reference may affect the brightness of the HP data (e.g. signal in the skull). This bias does not occur when showing adjacent maps (Fig. 9a, b). Here, anatomical outlines may help to provide reference (Fig. 9b).

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.