Parameters in hyperelastic model curve fitting

Many tasks in finite element analysis require curve fitting. Especially in materials engineering, when only experimental test data is available, engineers must perform curve fitting to determine the parameters in a given numerical model. Common types of material curve fitting include hyperelastic, plastic, viscous, core loss materials, and more. Press enter or click to view image in full size Unlike other types of curve fitting, curve fitting for finite element material models requires both a small error and physically valid fitted parameters; furthermore, these parameters must be conducive to convergence in subsequent nonlinear finite element computation. As a result, there are high standards placed on the algorithms and details of parameter fitting. Rarely do articles discuss the parameters, initial values, and reasonable bounds for hyperelastic model parameters despite their importance when implementing hyperelastic model curve fitting and verifying the rationality of fitted parameters. From a development perspective, this article discusses how to better perform curve fitting for hyperelastic material models, with a particular focus on setting parameter bounds and initial values, as well as identifying the characteristics of reasonable parameters. Due to the large number of hyperelastic models and order-dependent variation, this article describes the computational details of each hyperelastic model in separate sections to avoid confusion. 1. neo-Hookean The neo-Hookean model is first because it is one of the most simple and classic hyperelastic models. Many other models can be simplified to the neo-Hookean model under specific parameter conditions. Its strain energy density function W is expressed as: where μ is the shear modulus. During curve fitting computation, the shear modulus must be greater than zero, i.e., μ>0. The initial value is set to μ0​=E/3, where Young’s modulus E can be obtained from the slope of the stress-strain curve near the origin (strain < 5%). Press enter or click to view image in full size To check for physical validity after parameter fitting, you can refer below for the approximate correspondence between common hyperelastic materials and μ: Material Type | Shear Modulus μ (MPa) Extra-soft Rubber | 0.3 ~ 0.5 Medium-hardness Rubber | 1.0 ~ 1.5 Hard Rubber | 3.0 ~ 5.0 Biological Soft Tissue | 0.01 ~ 0.1 A fitted parameter μ that differs significantly from the uniaxial tension data indicates that the material may exhibit a strong second strain invariant effect, in which the Mooney-Rivlin or other suitable models should be used instead. 2. Mooney-Rivlin The Mooney-Rivlin model is another one of the most commonly used hyperelastic models in engineering analysis. It is a direct extension of the neo-Hookean model and is characterized by its simple and flexible form. If the neo-Hookean is considered the entry-level model in hyperelasticity, then the Mooney-Rivlin would be the standard model for handling small-to-moderate deformations (100%-200% strain). Its simple mathematical form and linear dependence on material constants allows it to converge easily in Finite Element Analysis (FEA) calculations. Depending on the order, the Mooney-Rivlin model has several different forms with varying amounts of parameters, the most complex being the 9-parameter model. This section uses the 9-parameter Mooney-Rivlin model as an example to discuss parameter ranges and initial value settings. Its strain energy density function W includes all combinations from linear to cubic terms: Press enter or click to view image in full size Physically, the shear modulus μ>0. Therefore, μ=2(C10​+C01​)>0. The most robust constraint for practical calculations is C10​>0 and C01​≥0. In many rubber materials, C10​ is the dominant term, and C01​ is relatively small; C11​ is also usually positive. The third-order parameter C30​≥0. C03​, C21​, C12​ can capture the sharp stiffening at extremely large deformations (strain > 400%) and are generally 3–4 orders of magnitude smaller than the first-order parameters. Press enter or click to view image in full size Regarding initial conditions, C10​=0.9×(μ0​/2), C01​=0.1×(μ0​/2), with the initial shear modulus μ0​=E/3. The initial values for C11​, C20​, C02​, C30​, C21​ are often chosen to be small, typically 0.1% of C10​. The initial values for C12​ and C03​ can be set to zero. For most synthetic rubbers (e.g., NBR, EPDM, Neoprene), the parameters generally follow the order-of-magnitude rules below: Order Classification | Parameter Symbol | Typical Value Range (MPa) 1st-order (Linear) | C10​,C01 | ​0.1 ~ 2.0 2nd-order (Transition) | C20​,C11​,C02 | ​-0.2 ~ 0.2 3rd-order (Stiffening) | C30​,C21​,C12​,C03 | ​-0.05 ~ 0.1 The Mooney-Rivlin model has its limitations. When rubber is stretched to its limit (molecular chains fully extended), the stress rises sharply (the “upturn” phenomenon). The Mooney-Rivlin model cannot account for this stiffening at high strains; the Ogden or Gent model is recommended when simulating extremely large deformations (strain > 300%). If the experimental data is not comprehensive, a more robust model with fewer parameters is advised, such as the Yeoh model or the 5-parameter Mooney-Rivlin model. 3. Yeoh The Yeoh hyperelastic model is known for its simplicity and stability. It discards the second invariant I2​ (terms like C01​), which is highly susceptible to experimental error, retaining only the first invariant I1​. Therefore, it is more robust than the Mooney-Rivlin model when fitting large-deformation data. The strain energy density function W is as follows: Common Yeoh models are 1st to 3rd-order, with 1 to 3 parameters respectively. During fitting, C10​>0 is required since the material loses shear resistance and crashes the simulation immediately if C10​≤0. C20​ is usually negative; its magnitude is not too large and can be constrained to be greater than 1% of C10​. C30​ is positive and is used to capture the stiffening at large strain. Press enter or click to view image in full size The initial value settings are relatively simple: the initial value of C10​ is E/6, where Young’s modulus E is the slope of the test data at small deformations. C20​=−0.1×C10​, C30​=0.01×C10​. The parameters of the Yeoh model for common materials are as follows: Material Type/Hardness | C10​ (MPa) | C20​ (MPa) | C30​ (MPa) | Description Ultra-soft Silicone/Gel | 0.01 ~ 0.1 | –0.001 ~ -0.01 | 0.0001 ~ 0.01 | Very low initial modulus, stiffening occurs late. Natural Rubber (40A-50A) | 0.2 ~ 0.5 | –0.02 ~ -0.1 | 0.01 ~ 0.05 | Good ductility, obvious plateau in the mid-curve. Industrial-grade Rubber (60A-70A) | 0.6 ~ 1.2 | –0.1 ~ -0.5 | 0.1 ~ 0.5 | Medium stiffness, commonly used for tire sidewalls or seals. High-hardness Elastomer (80A+) | 1.5 ~ 5.0 | –0.5 ~ -2.0 | 0.5 ~ 2.0 | Stiffening occurs very early, curve is steep. 4. Polynomial The Polynomial hyperelastic model is also based on the power series expansion of the strain energy density function W with respect to the invariants I1​ and I2​. It is one of the most widely used frameworks in hyperelasticity theory, and many models (Mooney-Rivlin, Yeoh, Neo-Hookean) are essentially its special cases. The general strain energy density function W is expressed as: where N is the order of the model. Common polynomial models are 1st to 3rd-order, corresponding to 2, 5, and 9 parameters. Like the Mooney-Rivlin model, during fitting, it is necessary for C10​>0 and C10​>∣C01​∣. C30​ is also usually required to be positive. Generally, the higher the order, the smaller the magnitude of the values. The reasonable parameter ranges for hyperelastic materials are shown below. Parameter Category | Parameter Item | Typical Range (MPa) | Physical Behavior Description Base Modulus | C10 | ​0.1 ~ 0.6 | Dominates small-to-medium deformation stiffness, usually positive. Shear Correction | C01​ | 0.01 ~ 0.1 | Corrects shear behavior, usually positive. Mid-range Adjustment | C20 | ​-0.05 ~ 0.05 | Controls slope change in the mid-curve, often negative. Large-deformation Stiffening | C30 | ​0.001 ~ 0.01 | Controls the upturn stiffening at very high strains, must be positive. Coupling Terms | C11​,C02 | ​-0.01 ~ 0.01 | Cross-effects between I1​ and I2​, usually close to zero. For initial value settings, set C10​=E/6. C20​=−0.1×C10​, C30​=0.01×C10​ or smaller. Although the value of C30​ is small, its high power amplifies its impact on the end of the curve. 5. Arruda-Boyce The Arruda-Boyce model is unique in hyperelastic material simulation, and it is constructed based on statistical mechanics (molecular chain network theory). The strain energy density function W is expressed as: Press enter or click to view image in full size A key feature of this model is its low requirement for test data, and uniaxial data alone is sufficient to start. Unlike polynomial models that require multi-axial data to determine parameters, the Arruda-Boyce model’s parameters have clear physical definitions, in which relatively accurate μ and λm​ can usually be obtained from the high-quality uniaxial tension data. Become a member The value range for the initial shear modulus μ is straightforward; it must be greater than 0 and lie between 0.1~5 MPa, with common values for rubber materials ranging from 0.5~1.5 MPa. λm​ must be greater than and typically limited to [1.1, 20]. The model becomes extremely stiff when λm​ is too close to 1, leading to highly unstable numerical calculations. For most industrial rubbers, setting λm​ between 3.0 and 7.0 usually covers most working conditions. Press enter or click to view image in full size Initial value determination is also straightforward: set μ=E/3. λm​=5 is a robust starting point. Typically, the Arruda-Boyce parameters correspond to rubber materials as follows: Material | μ (MPa) | λm ​Soft Rubber (Shore 30A - 40A) | 0.3 ~ 0.6 | 5.0 ~ 8.0 Medium-hardness (Shore 50A - 60A) | 1.0 ~ 2.5 | 3.0 ~ 5.0 Hard Rubber (Shore 70A - 80A) | 3.0 ~ 8.0 | 1.5 ~ 2.5 6. Blatz-Ko The Blatz-Ko model is a simple model used to simulate foams and porous materials. The most commonly used form contains only one parameter, the shear modulus μ. Like other hyperelastic models, the shear modulus must be greater than zero. Poisson’s ratio ν is fixed as a theoretical property of this model is that its effective Poisson’s ratio is locked at approximately 0.25. This means the material is assumed to shrink significantly in volume when compressed, rather than expanding sideways like ordinary rubber. Do not use the Blatz-Ko model for solid rubber (Poisson’s ratio near 0.5, incompressible). As long as the Blatz-Ko model ensures that the shear modulus μ is positive (μ>0) and the material is indeed a compressible foam, fitting is usually straightforward. Since Poisson’s ratio ν is often fixed at 0.25, according to the linear elastic conversion relationship: μ = E/[2(1+ν)] = E/2.5 = 0.4E. The initial value can be set to 0.4E. Typically, the shear modulus corresponds to the materials as follows: Material Type | μ Value Range (MPa) | Description Ultra-light Polyurethane Foam | 0.05 ~ 0.5 | Extremely low density, highly compressible. Medium-density Foamed Rubber | 0.5 ~ 2.0 | Common for shock pads and packaging. High-hardness Closed-cell Foam | 2.0 ~ 10.0 | Structural support foam, high stiffness. Solid Compressible Rubber | 10.0+ | Special synthetic rubbers with voids. 7. Gent The Gent model is a physically based hyperelastic model, similar in behavior to the Arruda-Boyce model. Its greatest advantage is its simple form, which can effectively capture the nonlinear behavior of polymer chains reaching their limiting extension length using only two constants. The mathematical form of the Gent model is a logarithmic function, and the computational cost for curve fitting is lower than that of the Arruda-Boyce model. The Gent model has only two core parameters: μ (initial shear modulus) and Jm​ (limiting invariant). Its strain energy density function is defined as: Press enter or click to view image in full size The parameter μ must be greater than 0; Jm​ must be greater than the maximum measured value of I1​−3. If Jm​ is too small, causing (I1​−3)/Jm ​≥ 1, the logarithm becomes undefined, so the calculation will fail. The initial value of the shear modulus μ is simply determined by μ=E/3. Mathematically, Jm​ defines the upper limit for I1​−3, so its initial value must be slightly larger than the invariant value corresponding to the maximum deformation observed in the experiment. To determine this, find the maximum elongation ratio λmax​ in the experimental data and calculate the corresponding I1​. Note that the formula for I1​ varies slightly depending on the type of experiment. The initial value of Jm​ is set to 1.2×(I1_max​−3). The empirical parameters for materials of different hardness are as follows: Material Hardness/Type | Shear Modulus μ (MPa) | Limiting Constant Jm | ​Description Ultra-soft Gel/Soft Tissue | 0.01 ~ 0.15 | 0 ~ 150 | Extremely stretchable, stiffens very late. Natural Rubber (40A) | 0.3 ~ 0.6 | 30 ~ 80 | Excellent ductility, high stiffening point. Industrial Rubber (60A) | 1.0 ~ 2.0 | 10 ~ 25 | Medium stiffness, stiffens rapidly after 3–4x stretch. High-hardness Elastomer (80A) | 3.0 ~ 7.0 | 2 ~ 8 | Stiffens within a very short distance, increased brittleness. 8. Ogden The Ogden hyperelastic model is one of the most powerful and flexible hyperelastic models available today. It is constructed directly based on the principal stretches λi​. The strain energy density function of the Ogden model adopts a series superposition form expanded by order n, where n represents the model order. 1st-order Ogden: 2 parameters, suitable for small deformations and simple rubber materials, highest computational efficiency. 2nd-order Ogden: 4 parameters, balances accuracy and efficiency, mainstream for general industrial simulations. 3rd-order Ogden: 6 parameters, suitable for large-deformation, highly nonlinear complex elastomers, highest fitting accuracy. The Ogden model maintains solid predictive capability even when strains reach 700% or higher. The Ogden model has relatively more parameters, and higher-order versions require combined fitting of multiple sets of test data, such as uniaxial, pure shear, and equibiaxial, to prevent non-physical solutions and overfitting. Parameters must satisfy thermodynamic stability constraints, and reasonable bounds must be set during fitting. For example, set μ1​>∣μ2​∣>∣μ3​∣, 1<α1​<10, −10<α2​<10, −15<α3​<15. Press enter or click to view image in full size Initial values for curve fitting are often set as follows: α1​=2, α2​=2, α3​=6, μ1​=0.6μ0​, μ2​=0.1μ0​, μ3​=0.05μ0​, where the initial shear modulus μ0​=E/3. Young’s modulus E is obtained from the slope of the test data at low strains. The empirical data fitted by common Ogden hyperelastic models is as follows: Material Type | Parameter | Typical Value Range | Description Soft Rubber / Silicone | μ | 0.01 ~ 0.5 MPa | Low modulus, gentle curve. | α | -10 ~ 10 | Usually includes positive and negative exponents to fit tension/compression symmetry. Industrial Hard Rubber | μ | 1.0 ~ 10.0 MPa | High stiffness, large initial slope from combined μ. | α | 1.5 ~ 20 | Large positive exponents to simulate severe stiffening at large deformations. Biological Soft Tissue | μ | 0.001 ~ 0.1 MPa | Extremely soft, highly sensitive to small deformations. | α | 10 ~ 50 | Exponents are usually very large to simulate the "lock-up" effect of fibrous connective tissue. Conclusion Although curve fitting has relatively mature numerical methods, fitting material parameters for finite element computation requires the consideration of many physical constraints, ensuring that the parameters are conducive to convergence in subsequent finite element calculations. Oftentimes due to the lack of test data, curve fitting for hyperelastic materials becomes complex which places high demand on both developers and users. The main test data required for hyperelastic material model fitting includes uniaxial tension, equibiaxial tension, and pure shear deformation data. The author discusses this in detail in the article “Hyperelastic Material Models and Curve Fitting”. Currently, WELSIM can accurately fit hyperelastic model curves. The independent and free engineering software MatEditor and CurveFitter also feature the aforementioned curve fitting capabilities. Users can download and use this software directly.

超弹模型曲线拟合中的参数范围

在有限元分析中，有大量需要进行曲线拟合的工作。尤其是在材料领域，当工程师们只有一些材料的离散的测试数据，在给定数值模型后，需要通过对材料实验数据的曲线拟合，才能确定模型中的参数。常见的曲线参数拟合有，超弹材料，塑性材料，粘性材料，和磁芯损耗材料等等。和其他类型的曲线拟合不同的是，这些有限元材料模型的曲线拟合不仅需要误差小，同时要求拟合出的参数物理性正确，而且在后续的非线性有限元计算中，这些参数还需要易于收敛。因此对参数拟合的算法和细节，提出了比较高的要求。 Image 很少有文章讨论超弹模型参数的取值范围，初始值设定，以及参数值的合理范围。而这些都是开发超弹模型曲线拟合，验证拟合参数合理性的关键。本文从开发的角度，讨论如何更好地实现超弹材料模型的曲线拟合。尤其是参数边界和初始值的设定，以及合理的参数应该符合怎样特征。 由于超弹模型众多，每种超弹模型还有会根据阶数的不同，有多种形式。为了避免混淆，本文按章节分别描述每种超弹模型的计算细节。 1. neo-Hookean 将neo-Hookean作为第一个讨论的模型，是因为neo-Hookean是最经典、最简单的超弹模型之一。很多其他模型，可以在某种特殊参数下， 简化为neo-Hookean模型。其应变能密度函数 W的表达式为 Image 其中μ是剪切模量。 在曲线拟合计算时，要求剪切模量必须大于零，即μ>0。初始值设置为μ0=E/3，杨氏模量E可以从应力-应变曲线在原点附近（应变 < 5%）的斜率获得。 Image 参数拟合后，可以查验是否符合真实物理，可以参考以下常见超弹材料与μ的近似对应关系。 材料类型 剪切模量 μ (MPa) 极软橡胶 0.3 ~ 0.5 中等硬度橡胶 1.0 ~ 1.5 硬橡胶 3.0 ~ 5.0 生物软组织 0.01 ~ 0.1 如果拟合出的 参数μ与单轴拉伸数据差异很大，说明材料可能表现出明显的第二应变能不变量效应，此时应改用 Mooney-Rivlin 或其他模型。 2. Mooney-Rivlin Mooney-Rivlin也是在工程分析中最常用的超弹模型之一。是neo-Hookean模型的直接扩展。特点是形式简单且灵活，如果说 neo-Hookean 是超弹界的“入门款”，那么Mooney-Rivlin就是处理中低度变形（100%-200%应变）的“标准款”。 由于其数学形式简单且线性相关于材料常数，在有限元分析（FEA）计算中容易收敛。 根据阶数的不同，Mooney-Rivlin有几种不同的形式，含有不同数量的参数，其中最高阶的9参数模型最为复杂，本节就以九参数Mooney-Rivlin模型为例，讨论参数的范围与初始值设置。其应变能密度函数W也包含了从一次项到三次项的所有组合： Image 物理上要求剪切模量μ>0。因此有μ=2(C10+C01)>0，实际计算时最稳妥的约束：C10>0 且C01>=0。在许多橡胶材料中，C10是主导项，而 C01较小。C11 通常也为正值。三阶参数C30>=0。C03, C21, C12可以捕捉极大变形（应变 > 400%）时的急剧硬化。一般比一阶参数小 3-4 个数量级。 Image 关于初始条件，C10=0.9*(μ0/2)，C01 = 0.1*(μ0/2)，初始剪切模量 μ0=E/3。C11，C20，C02，C30, C21的初始值往往选择一个很小的值，取C10的0.1%。C12和C03的初始值可以设为零。对于大多数合成橡胶（如 NBR, EPDM, 氯丁橡胶），参数大致遵循以下数量级规律： 阶数分类 参数符号 典型取值范围 一阶 (线性项) C10, C01 0.1~2.0 二阶 (过渡项) C20, C11, C02 -0.2~0.2 三阶 (硬化项) C30, C21, C12, C03 -0.05~0.1 Mooney-Rivin也有一定局限性。当橡胶拉伸到极限（分子链拉直）时，应力会急剧上升（即上翘现象）。Mooney-Rivlin 无法模拟这种高应变下的硬化。如果需要模拟极大的变形（应变 > 300%），建议使用 Ogden 或 Gent 模型。同时，如果实验数据不够全面，建议退而求其次使用参数更少、更鲁棒的模型，如Yoeh模型或5 参数 Mooney-Rivlin模型。 3. Yeoh Yeoh 超弹模型以其简洁性和稳定性著称，舍弃了受实验误差干扰较大的第二不变量 I2（即 C01 等项），仅保留第一不变量 I1，因此在拟合大变形数据时比 Mooney-Rivlin 模型更加稳健。应变能密度函数 W如下 Image 常用的Yeoh模型有1-3阶，对应1-3参数。拟合计算时，需要C10>0。如果 C10<= 0，材料将失去抗剪切能力，导致模拟计算立即崩溃。 C20 通常为负，且绝对值不会太大。拟合计算时，可以控制其绝对值大于在C10的1%。C30 为正值，用于捕捉末端硬化 。 Image 初始值设置也相对简单，C10的初始值为E/6，杨氏模量E是测试数据在小变型时的斜率。C20=-0.1*C10, C30=0.01*C10。 常见材料的Yeoh模型的参数如下： 材料类型/硬度 C10 (MPa) C20 (MPa) C30 (MPa) 特征描述 超软硅胶/凝胶 0.01~0.1 -0.001~-0.01 0.0001~0.01 初始模量极低，硬化发生很晚。 天然橡胶 (40A-50A) 0.2~0.5 -0.02~-0.1 0.01~0.05 延展性好，曲线中段有明显的平缓区。 工业级橡胶 (60A-70A) 0.6~1.2 -0.1~-0.5 0.1~0.5 刚度中等，常用于轮胎侧壁或密封圈。 高硬度弹性体 (80A+) 1.5~5.0 -0.5~-2.0 0.5~2.0 硬化发生极早，曲线陡峭。 4. 多项式 多项式超弹模型是也一种基于应变能密度函数 W 对不变量 I1 和 I2 进行幂级数展开的模型。 它是超弹性理论中应用最广泛的框架之一，许多模型（如 Mooney-Rivlin, Yeoh, Neo-Hookean）本质上都是它的特殊形式。 其通用的应变能密度函数 W表示为： Image N是模型的阶数，常用的多项式模型有1-3阶，对应2，5和9参数模型。和Mooney-Rivlin模型一样，拟合计算时，需要保证C10 > 0, 同时保证C10 > |C01|。C30通常也要求为正数。且阶数越高，数值量级通常越小。对于超弹材料，合理的参数范围如下表所示。 参数类别 参数项 典型范围 (MPa) 物理行为说明 基础模量 C10 0.1~0.6 主导中小变形刚度，通常为正值。 剪切修正 C01 0.01~0.1 修正剪切行为，通常为正值。 中段调节 C20 -0.05~0.05 控制曲线中段的斜率变化，常为负值。 大变形硬化 C30 0.001~0.01 控制极高应变下的上扬硬化，必须为正。 耦合项 C11, C02 -0.01~0.01 交叉影响 I1 和 I2通常接近于零。 初始值的设定上，令C10=E/6。C20=-0.1*C10。C30=0.01*C10或更小，C30值虽小，但在高次方放大下对曲线末端影响很大。 Image 5. Arruda-Boyce Arruda-Boyce 模型（通常被称为“八链模型”）在超弹性材料仿真中非常独特。是基于统计力学（分子链网络理论）构建的。应变能密度函数 W表示为： Image 模型最大的特点就是对测试数据要求低，单轴数据足以起步：与需要多轴数据才能锁定的多项式模型不同，Arruda-Boyce 模型的参数具有明确物理定义，通常仅靠高质量的单轴拉伸数据就能获得相对准确的μ和λm。 初始剪切模量μ的取值范围很直接，必须大于 0,值限定在0.1~5MPa，常见橡胶材料参考值是0.5~1.5MPa之间。λm必须大于1，值通常限定在[1.1, 20]。如果λm太接近 1，模型会变得极其僵硬，导致数值计算极其不稳定。对于大多数工业橡胶，λm设定在 3.0 到 7.0 之间通常能覆盖绝大多数工况。 Image 初始值的确定也很直接，令μ=E/3。λm=5是一个稳健的起点。 通常情况下，Arruda-Boyce 参数对于橡胶材料有如下对应关系： 材料 μ （MPa） λm 软橡胶 (Shore 30A – 40A) 0.3~0.6 5.0~8.0 中等硬度 (Shore 50A – 60A) 1.0~2.5 3.0~5.0 硬橡胶 (Shore 70A – 80A) 3.0~8.0 1.5~2.5 6. Blatz-Ko Blatz-Ko模型是一种形式简单的用于模拟泡沫和多孔材料的模型。最常用的形式仅包含一个参数：剪切模量μ。 Image 和其他超弹模型一样，剪切模量必须要大于零。泊松比ν是固定的： 该模型的一个理论特性是其有效泊松比被锁定在 0.25 左右。这意味着假设材料在受压时体积会显著缩小，而不是像普通橡胶那样向四周扩张。如果材料是实心橡胶（泊松比接近 0.5，不可压缩），则不要使用 Blatz-Ko 模型。 Blatz-Ko 模型只要确保剪切模量μ是正值，即μ>0且该材料确实是可压缩的泡沫类材料，拟合通常非常直接。 由于泊松比ν通常被固定为0.25。根据线弹性转换关系：μ = E/2(1+ν) = E/2.5 = 0.4E。初始值可以设定为0.4E。 材料类型 μ取值范围（MPa） 说明 超轻质聚氨酯泡沫 0.05~0.5 密度极低，极易压缩变形。 中等密度发泡橡胶 0.5~2.0 常见的减震垫、包装防护材料。 高硬度闭孔泡沫 2.0~10.0 结构支撑用泡沫，刚度较高。 实心可压缩橡胶 10.0+ 某些含气泡的特殊合成橡胶。 7. Gent Gent 模型是一种基于物理意义的超弹性模型，它与 Arruda-Boyce 模型非常相似现象。 该模型最大的优势在于形式简单，仅通过两个常数就能很好地捕捉聚合物链达到极限拉伸长度时的非线性行为。Gent 模型的数学形式是对数函数，拟合计算的开销比 Arruda-Boyce 模型要小。 Gent 模型只有两个核心参数：μ（初始剪切模量） 和Jm（极限不变量）。其应变能密度函数定义为： Image 参数μ必须严格大于 0。 Jm必须严格大于 I1-3 的最大测量值。 如果 Jm 太小，导致 (I1-3)/Jm <= 1，对数函数ln 会失去意义，计算会报错。 剪切模量μ的初始值确定比较简单，μ=E/3。Jm 在数学上定义了第一不变量I1-3的上限，所以初始值必须略大于实验中观测到的最大变形所对应的不变量值。找到实验数据中最大的伸长率λ_max，计算最大变形对应的I1，注意这里I1的计算，根据实验类型的不同，公式稍有不同。Jm初始值定为1.2 * (I1_max – 3)。 不同硬度材料的经验参数参考如下： 材料硬度/类型 剪切模量μ(MPa) 极限常数Jm 说明 超软凝胶 / 软组织 0.01~0.1 50~150 极易拉伸，极晚进入硬化阶段。 天然橡胶 (40A) 0.3~0.6 30~80 延展性极佳，硬化点较高。 工业橡胶 (60A) 1.0~2.0 10~25 刚度中等，拉伸 3-4 倍后迅速硬化。 高硬度弹性体 (80A) 3.0~7.0 2~8 极短距离内即发生硬化，脆性增加。 8. Ogden Ogden超弹模型是目前最强大、最灵活的超弹性模型之一。Ogden 模型直接基于主拉伸比（Principal Stretches）λi来构建。 Ogden 模型的应变能密度函数，采用级数叠加形式，按阶数n展开，n代表模型阶数。 Image 1 阶 Ogden：2 个参数，适合小变形、简单橡胶材料，计算效率最高。 2 阶 Ogden：4 个参数，兼顾精度与效率，是常规工业仿真的主流选择。 3 阶 Ogden：6 个参数，适合大变形、高非线性复杂弹性体，拟合精度最高。 在应变达到 700% 甚至更高的情况下，Ogden 模型依然能保持良好的预测能力。 Ogden模型参数较多，高阶版本需要单轴、纯剪切、等双轴等多组试验数据联合拟合，否则容易出现非物理解、过拟合。参数需要满足热力学稳定性约束，拟合时必须设置合理边界。如令μ1>|μ2|>|μ3|，1<α1<10，-10<α2<10, -15<α3<15。 Image 曲线拟合的初始值常做如下设置：α1=2, α2=2, α3=6, μ1=0.6*μ0, μ2=0.1*μ0, μ3=0.05*μ0, 其中初始剪切模量μ0=E/3。杨氏模量E通过低应变下的测试数据斜率得到。 常见的Odgen超弹模型拟合出的经验数据如下： 材料类型 参数符号 典型取值范围 特征说明 软橡胶 / 硅胶 μ 0.01~0.5 MPa 模量较小，曲线平缓。 α -10~10 通常包含正负指数以拟合拉伸和压缩对称性。 工业级硬橡胶 α 1.0~10.0 MPa 刚度高，μ组合后的初始斜率大。 α 1.5~20 较大的正指数用于模拟大变形下的剧烈硬化。 生物软组织 μ 0.001~0.1 MPa 极软，对小变形高度敏感。 α 10~50 指数通常很大，以模拟纤维结缔组织的“锁死”效应。 总结 虽然曲线拟合有着比较成熟的数值计算方法，但是在有限元计算的材料参数拟合上，需要考虑诸多的物理因素，同时要保证拟合出的参数在后续有限元计算中易于收敛。很多时候，由于测试数据的缺乏，使得超弹材料的曲线拟合变得复杂，对开发者和使用者都有比较高的要求。 超弹材料模型拟合所需的测试数据主要有，单轴拉伸，等双轴拉伸 ，和纯剪切变形数据。笔者在《超弹材料模型及其曲线拟合》一文已经有了详细描述。 目前，WELSIM已经具备了准确拟合超弹模型曲线的功能，同时，独立运行的免费工程软件MatEditor和CurveFitter也具有以上曲线拟合功能。用户可以直接下载并使用。

Generate the contact and boundary conditions for the CalculiX solver

CalculiX is a well-known open-source Finite Element Analysis software that focuses on structural mechanics (static, dynamic, nonlinear) and thermal analysis. It is suitable for scientific research, teaching, and small-to-medium engineering scenarios. The core advantages of CalculiX are as follows: Open-Source and Free. It is licensed under the GNU General Public License (GPL) which allows individuals, enterprises, and research institutions to freely use, modify, and distribute the source code. INP File Format Compatibility. It supports the INP input file format used by Abaqus. Most Abaqus models can be directly imported into CalculiX for solving, which lessens the learning curve for users. Lightweight Solver. The solver features streamlined code and requires less hardware resources compared to commercial software, making it ideal for fast computation of small-to-medium scale models. Seamless Third-Party Integration. It offers excellent compatibility with third-party preprocessors and postprocessors, facilitating the integration of pre/post-processing tools with the solver. Press enter or click to view image in full size In a previous article titled Generate CalculiX solver files using WELSIM, the author introduced methods to quickly generate CalculiX input scripts. Currently, WELSIM serves as a robust preprocessor to generate input files for CalculiX and greatly simplifies the process of performing finite element analysis with CalculiX. This article elaborates on the commands for generating contact and boundary conditions. Contact In the finite element structural analysis of multi-body systems, contact setup is an essential step. This section describes the generation of two common contact types: bonded contact and separable contact. Bonded Contact Bonded contact is a contact condition used to simulate a fully rigid connection between multiple bodies. It is applicable to scenarios where there is no relative sliding, separation, or gaps between components (e.g., welded joints, adhesive bonds, cured interference fits, and rigid bolted connections). Bonded contact enforces identical nodal displacements across the contact interface, essentially “fusing” multiple independent components into a single integrated structure. The setup for bonded contact in WELSIM is shown in the figure below. When converting to CalculiX commands, the *TIE keyword is used. The generated input commands are as follows: *Tie, Name=ID_15, Position tolerance=0.01, Adjust=No Target_Surface_15, Master_Surface_15 The surface elements defined in the mesh file are as follows: ** Box_to_Box_1 *SURFACE, TYPE=ELEMENT, NAME=Master_Surface_15 70, S2 63, S2 74, S2 61, S2 66, S1 69, S2 ** Box_to_Box_1*SURFACE, TYPE=ELEMENT, NAME=Target_Surface_15 166, S3 185, S2 160, S3 173, S2 171, S4 169, S3 Separable Contact Separable contact is a nonlinear contact condition used to simulate three interface states between multiple bodies: contact, sliding, and separation. Its core characteristic is that the contact interface only transmits normal pressure (no tensile forces) and allows free sliding or sticking in the tangential direction. It is suitable for scenarios involving dynamic contact and separation between components (e.g., gear engagement, sheet metal forming, bearing rolling, and mechanical impact). Based on friction assumptions, it can be further grouped into frictionless contact and frictional contact. The setup for seperable contact in WELSIM is shown in the figure below. When converting to CalculiX commands, the *CONTACT keyword is used. It is important to note that CalculiX requires that the master surface (also known as the independent surface) is defined using elements. The slave surface (also known as the dependent surface) can be defined using either elements or nodes. The corresponding generated CalculiX commands are as follows: ** Box_to_Box 1 *Surface interaction, Name=SurfaceInteraction_15 *Surface behavior, Pressure-overclosure=Hard *Friction 0.15 *Contact Pair, Interaction=SurfaceInteraction_15, Type=Surface to surface Target_Surface_15, Master_Surface_15 Boundary conditions in structural analysis Setting boundary conditions is a critical step in finite element analysis and a key function of preprocessing software for generating solver commands. WELSIM supports a wide range of CalculiX boundary conditions. This section demonstrates how a preprocessor generates boundary conditions for structural analysis and presents the corresponding solver commands. 1. Constraints or Displacement Displacement boundary conditions are essential (Dirichlet) boundary conditions used to constrain the nodal displacement degrees of freedom (DOFs) of a model. Their role is to restrict the rigid-body motion of the structure and simulate real-world support constraints (e.g., fixed supports, pin supports), serving as the foundational conditions for static and dynamic structural analysis. In regards to 3D models: For solid elements, translational DOFs in the X/Y/Z directions (Ux/Uy/Uz) can be constrained. For shell elements, additional rotational DOFs about the X/Y/Z axes (Rx/Ry/Rz) can be constrained. The corresponding generated CalculiX solver commands are as follows: ** BC Name: Displacement *Boundary ID_22, 1, 1, 1 *Boundary ID_22, 2, 2, 2 *Boundary ID_22, 3, 3, 3 *Boundary ID_22, 6, 6, 0.523599 2. Pressure Pressure boundary conditions are surface loads applied to structural surfaces, classified as distributed loads. Their essence is to simulate engineering scenarios involving fluid pressure, contact pressure, gas loads, etc., by specifying the normal force per unit area (e.g., internal pressure in pressure vessels, wind loads, lateral soil pressure). Pressure is one of the most common boundary conditions in structural finite element analysis. The sign convention for pressure determines the load direction. Positive Pressure: The load direction points toward the interior of the surface, causing a compressive effect on the structure (e.g., internal pressure in pressure vessels). Negative Pressure: The load direction points away from the exterior of the surface, causing a tensile effect on the structure (e.g., atmospheric pressure on the outer wall of a vacuum vessel). The pressure boundary condition setup is shown in the figure below. The corresponding generated solver commands are as follows: ** BC Name: Pressure *DLOAD ID_20, P, 123 3. Force Force boundary conditions are loads applied to geometric key points, edges, or surfaces, used to simulate localized loading scenarios in engineering (e.g., bolt preload, lifting eye tension, pin forces). Force is one of the most fundamental boundary conditions in structural mechanics. However, applying force on geometric surfaces in the finite element numerical method is non-trivial. It typically involves coupling the DOFs of multiple nodes to a single reference node, which then used to distribute a concentrated load across a group of nodes (e.g., applying the force to the reference node, which is then automatically distributed to the coupled group). As shown in the figure, in the force setup of WELSIM, the user can select surfaces, edges, or points of the geometry and specify the force magnitudes in the three coordinate directions. The generated solver commands for the force are as follows: *COUPLING, REF NODE=46, SURFACE=ID_17, CONSTRAINT NAME=c17 *KINEMATIC 1 3 ** BC Name: Force *CLOAD, OP=NEW 46, 1, 100 *CLOAD, OP=NEW 46, 2, 200 *CLOAD, OP=NEW 46, 3, 300 4. Nodal Force Nodal force boundary conditions are concentrated loads applied directly to single or multiple nodes. They represent the most fundamental form of load application after meshing. Nodal forces are more oriented toward direct post-meshing operations, making them suitable for refined load control or special scenario simulations. Become a member Nodal forces are loads applied directly to the mesh nodes. The total applied force is the superposition of forces based on the number of nodes selected. The generated solver commands are straightforward, using *CLOAD to apply loads based on the selected nodes: ** BC Name: Nodal Force *CLOAD ID_18, 1, 10 *CLOAD ID_18, 2, 20 *CLOAD ID_18, 3, 30 ** Nodal_Force *NSET, NSET=ID_18 82, 84, 86, 88, 95, 96, 103, 104, 107, 108, 111, 112, 145, 146, 147, 148, 149, 150 5. Gravity Gravity boundary conditions fall under the category of body forces, used to simulate the effect of the Earth’s gravitational field on the structure. They work by assigning an inertial force proportional to the mass to each element of the model. Gravity is applicable to almost all structural analyses involving self-weight effects (e.g., building beams and slabs, mechanical components, ground-based spacecraft conditions). Gravity is essentially a body force induced by acceleration, and its magnitude is directly proportional to the mass of the elements. The setup in WELSIM is shown in the figure below. The generated solver commands are as follows: ** DC Name: Earth Gravity *Dload Eall, Grav, 9.807, 0, 0, -1 6. Rotational Velocity Rotational velocity boundary conditions are used to simulate the scenario where a structure rotates uniformly around an axis. This condition induces centrifugal forces and Coriolis forces (which need to be considered for high-speed rotation) within the structure. Classified as inertial loads, they are widely applied in the analysis of rotating machinery (e.g., flywheels, impellers, centrifuges, turbine rotors). The input parameters include the axis of rotation, the origin of rotation, and the magnitude of the rotational velocity. The input setup is shown in the figure below. Press enter or click to view image in full size The generated solver commands are as follows: ** DC Name: Rotational Velocity *Dload ID_21, CENTRIF, 100, 5, 5, 0, 0, 0, 1 Common boundary in thermal analysis 1. Temperature In finite element thermal analysis, temperature boundary conditions are essential (Dirichlet) boundary conditions that directly define the temperature values of specific regions of the model. They are suitable for scenarios where the surface or nodal temperatures are known (e.g., the wall of a constant-temperature water tank, components in contact with a large heat source, temperature-controlled surfaces of thermal equipment). In the heat conduction governing equation, temperature boundary conditions act as enforced constraints, directly fixing the temperature DOFs of the boundary nodes. The solver prioritizes satisfying these constraints before calculating the internal temperature field distribution. The setup for WELSIM is shown in the figure below. The generated CalculiX solver commands are as follows, ** BC Name: Temperature *Boundary, op=New ID_23, 11, 11, 0 Heat Flux The heat flux boundary condition is a thermal load that directly defines the heat flow rate per unit area across a solid surface. It is suitable for scenarios where the heat input or output on a surface is known (e.g., the surface of an electric heater, solar collectors, walls subjected to high-temperature gas impingement). Heat flux is a core parameter where a positive direction indicates heat flowing into the solid, and a negative direction indicates heat flowing out of the solid. The input interface is shown below. The generated CalculiX solver commands are as follows: ** BC Name: Heat Flux *Dflux ID_24, S, 0.5 Body Heat Flux Body heat flux is a volumetric thermal load used to simulate heat generation or dissipation from internal heat sources within a solid. Unlike surface heat flux boundary conditions, body heat flux is applied directly to the volume elements of the model instead of acting on surfaces, making it suitable for temperature field analysis driven by internal heat sources. The unit of body heat flux in the SI system is W/m³, where positive values indicate heat generation, and negative values indicate heat absorption. The user input interface is shown below. The generated CalculiX solver commands are as follows: *DFLUX ID_31,BF,10. Convection Convection boundary conditions, also known as film boundary conditions, are used to simulate convective heat transfer between a solid surface and the adjacent fluid (gas or liquid). They are one of the most commonly used boundary conditions in engineering thermal analysis, with applications in electronic device cooling, pipe heat exchange, automotive engine cooling, and other scenarios. The convection heat transfer coefficient is the core parameter characterizing the intensity of convective heat transfer, while the ambient temperature parameter refers to the temperature of the mainstream fluid far from the solid surface. The input interface in WELSIM is shown below. The generated CalculiX solver commands are as follows: ** BC Name: Convection *Film, op=New ID_29_S1, F1, 80, 123 ID_29_S2, F2, 80, 123 ID_29_S3, F3, 80, 123 ID_29_S4, F4, 80, 123 Thermal radiation The radiation boundary condition is used to simulate heat transfer between a solid surface and the surrounding environment (potentially other solid surfaces) via thermal radiation. It is fundamentally governed by the Stefan–Boltzmann law and is suitable for heat exchange scenarios in media such as vacuum or gas (e.g., spacecraft thermal control, heat dissipation from high-temperature equipment). As a nonlinear boundary condition, the heat flux is proportional to the fourth power of the temperature, requiring iterative calculations for solution. Users are required to input two parameters: emissivity (a value between 0 and 1) and the ambient temperature. The input interface is shown below. The generated commands for the thermal radiation boundary condition are as follows: ** BC Name: Radiation *Radiate, op=New ID_30_S1, R1, 30, 0.9 ID_30_S2, R2, 30, 0.9 ID_30_S3, R3, 30, 0.9 ID_30_S4, R4, 30, 0.9 Conclusion In finite element analysis, there is a wide variety of boundary conditions that involve not only the parameters of the conditions themselves, but also the selected elements or nodes. For preprocessing software, generating commands for various boundary conditions is a complex undertaking. WELSIM is capable of establishing a comprehensive range of boundary conditions and can quickly generate any input files required by the CalculiX solver, which can be directly used for computation. Users can leverage WELSIM as a preprocessor for CalculiX. CalculiX also supports some additional boundary conditions such as *CONSTRAINT and *EQUATION. These were not discussed in this article as they are less frequently used in standard analyses. The input file format of CalculiX is highly similar to that of Abaqus; therefore, the content described in this article can also be applied to the Abaqus solver. Disclaimer: WelSim and the author are not affiliated with CalculiX or Abaqus, and have no direct relationship with the development teams or organizations behind CalculiX and Abaqus. The use of open-source software names and images herein is solely for reference in technical blog articles and software usage guidance.

生成CalculiX求解器的接触与边界条件

CalculiX 是一款著名的开源有限元分析（FEA）软件，主打结构力学（静力、动力、非线性）与热分析，适用于科研、教学和中小型工程场景。CalculiX的核心优势有：开源免费。个人、企业、科研机构均可自由使用、修改和分发代码，遵循 GNU GPL 许可证。兼容支持 Abaqus 的 INP 输入文件格式，多数 Abaqus 模型也可直接导入 CalculiX 求解，减少学习曲线。求解器代码精简，对硬件资源要求低于商业软件，适合中小型模型的快速计算。对第三方前后端软件支持也非常好，易于前后处理器与求解器的整合。 Image 笔者曾在《使用WELSIM生成CalculiX求解器文件》一文中，介绍了如何快速地生成CalculiX输入文件。目前，WELSIM已经可以很好的作为CalculiX的前处理器生成输入文本，极大的方便用户使用CalculiX进行有限元分析。本文以更详细的方式介绍生成接触与边界条件的命令。 接触 在多体有限元结构分析中，接触是必要的设置之一。本章节描述生成两种常见的接触类型：绑定接触和可分离接触。 绑定接触 绑定接触是一种用于模拟多体之间完全刚性连接的接触条件，适用于部件间无相对滑移、无分离、无间隙的工况（如焊接、胶粘、过盈配合后固化、螺栓紧固的刚性连接）。绑定接触会强制接触界面的节点位移完全一致，等同于将多个独立部件 “融合” 为一个整体结构。 绑定接触设置如下图所示， Image 转换成CalculiX命令时，使用*TIE 命令。生成的输入命令如下 *Tie, Name=ID_15, Position tolerance=0.01, Adjust=No Target_Surface_15, Master_Surface_15 在网格文件中，定义的表面单元如下： ** Box_to_Box_1 *SURFACE, TYPE=ELEMENT, NAME=Master_Surface_15 70, S2 63, S2 74, S2 61, S2 66, S1 69, S2 ** Box_to_Box_1 *SURFACE, TYPE=ELEMENT, NAME=Target_Surface_15 166, S3 185, S2 160, S3 173, S2 171, S4 169, S3 可分离接触 可分离接触是一种用于模拟多体接触界面允许接触、滑移、分离三种状态的非线性接触条件，核心特点是接触界面仅传递法向压力，无拉力作用，且切向可自由滑移或黏着，适用于部件间存在动态接触与分离的工况（如齿轮啮合、冲压成型、轴承滚动、机械碰撞）。根据摩擦力假设，又分为无摩擦接触和有摩擦接触。 WELSIM设置接触如下图所示， Image 转换成CalculiX命令时，使用CONTACT命令。这里值得注意的是，CalculiX规定，主接触面又称作独立接触面，必须是以单元定义的表面。被动接触面又称作是依赖接触面，可以是单元定义的表面，也可以是节点定义的表面。 生成对应的CalculiX命令为 ** Box_to_Box 1 *Surface interaction, Name=SurfaceInteraction_15 *Surface behavior, Pressure-overclosure=Hard *Friction 0.15 *Contact Pair, Interaction=SurfaceInteraction_15, Type=Surface to surface Target_Surface_15, Master_Surface_15 结构分析中的边界条件 设置边界条件是有限元分析中重要的步骤，也是前处理软件生成求解器命令的重要部分。WELSIM支持大量CalculiX的边界条件。本章节演示前处理器生成结构分析中的边界条件，并显示对应的求解器命令。 1. 固定边界或位移 位移边界条件是用于约束模型节点位移自由度的本质边界条件，其作用是限制结构的刚体运动、模拟实际支撑约束（如固定支座、铰支座），是结构静力、动力分析的基础条件。 3维模型中，对于实体单元，可以施加X/Y/Z方向的平移（Ux/Uy/Uz），对于壳单元，可以施加额外的绕 X/Y/Z 轴转动的转动（Rx/Ry/Rz）。 Image 对应生成的CalculiX求解器命令是， ** BC Name: Displacement *Boundary ID_22, 1, 1, 1 *Boundary ID_22, 2, 2, 2 *Boundary ID_22, 3, 3, 3 *Boundary ID_22, 6, 6, 0.523599 2. 压力 压力边界条件是一种作用于结构表面的面载荷，属于分布载荷，其本质是通过指定单位面积上的法向力来模拟流体压力、接触压力、气体载荷等工程场景（如压力容器内壁压力、风载荷、土壤侧压力）。 压力是结构有限元分析中常见的边界条件之一。压力的符号规则决定载荷方向：正压力：载荷方向指向表面内侧，对结构产生压缩效应（如压力容器内壁压力）；负压力：载荷方向背离表面外侧，对结构产生拉伸效应（如真空容器外壁的大气压力）。压力边界条件如图所示， Image 对应生成的求解器命令是， ** BC Name: Pressure *DLOAD ID_20, P, 123 3. 力 力边界条件是作用于几何关键点、线、或者面的载荷，用于模拟工程中局部位置的受力工况（如螺栓预紧力、吊环拉力、销钉作用力）。 力是结构力学中最常见的边界条件，然而， 对于有限元数值方法中，几何面上的力边界条件的施加并不简单，通常是将多个节点的自由度耦合到一个参考节点，常用于将集中载荷分散到节点组（如把力加在参考节点，自动分配到耦合组）。 如图所示，在前处理器的力设置中，用户可以选定几何体的面，线，或点，并给定三个方向上的力大小。 Image 得到力的求解器命令如下， *COUPLING, REF NODE=46, SURFACE=ID_17, CONSTRAINT NAME=c17 *KINEMATIC 1 3 ** BC Name: Force *CLOAD, OP=NEW 46, 1, 100 *CLOAD, OP=NEW 46, 2, 200 *CLOAD, OP=NEW 46, 3, 300 4. 节点力 节点力边界条件是直接作用于单个或多个节点的集中力载荷，是有限元离散化后最基础的载荷施加形式，与之前提到的集中力边界条件高度关联，但更偏向于网格离散后的直接操作，常用于精细化载荷控制或特殊工况模拟。 节点力是直接施加在网格节点上的力，根据节点数量大小，综合施加的力会叠加。 Image 生成的求解器命令也非常简单明了，根据所选的节点，施加CLOAD载荷。 ** BC Name: Nodal Force *CLOAD ID_18, 1, 10 *CLOAD ID_18, 2, 20 *CLOAD ID_18, 3, 30 ** Nodal_Force *NSET, NSET=ID_18 82, 84, 86, 88, 95, 96, 103, 104, 107, 108, 111, 112, 145, 146, 147, 148, 149, 150 5. 重力 重力边界条件属于体积力（Body Force） 范畴，用于模拟地球引力场对结构的作用，是通过给模型的每个单元分配与质量成正比的惯性力，适用于几乎所有包含自重影响的结构分析（如建筑梁板、机械部件、航天器地面工况）。 重力的本质是加速度场作用下的体积力，其大小与单元的质量直接相关。前处理器设置如下， Image 生成的求解器命令如下图所示。 ** DC Name: Earth Gravity *Dload Eall, Grav, 9.807, 0, 0, -1 6. 旋转角速度 旋转角速度边界条件用于模拟结构绕某一轴做匀速转动的工况，会在结构内部产生离心力和科里奥利力（高速转动时需考虑），属于惯性载荷范畴，广泛应用于旋转机械分析（如飞轮、叶轮、离心机、涡轮转子）。 输入的参数有旋转轴，旋转原点，和角速度的大小。输入界面如下图所示。 Image 生成的求解器命令如下所示。 ** DC Name: Rotational Velocity *Dload ID_21, CENTRIF, 100, 5, 5, 0, 0, 0, 1 热分析中的常见边界与初始条件 1. 温度 在有限元热分析中，温度边界条件是直接定义模型特定区域温度值的约束条件，属于本质边界条件（Dirichlet 边界条件），适用于已知表面或节点温度的场景（如恒温水箱壁面、与大热源接触的部件、温控设备的设定温度面）。 在热传导控制方程中，温度边界条件是强制约束，会直接固定边界节点的温度自由度，求解时优先满足该约束，再计算内部温度场分布。 前处理器设置如下， Image 生成的CalculiX求解器命令如下， ** BC Name: Temperature *Boundary, op=New ID_23, 11, 11, 0 2. 热流密度 热流密度（Heat Flux）边界条件是直接定义单位面积上通过固体表面的热流量的边界设置，属于热载荷的一种，适用于已知表面热流输入 / 输出的场景（如电加热器表面、太阳能集热板、高温燃气冲刷壁面）。 热流密度是核心参数，正方向为热量流入固体，负方向为热量流出。输入界面如下： Image 生成的CalculiX求解器命令如下： ** BC Name: Heat Flux *Dflux ID_24, S, 0.5 3. 体热流密度 体积热流是一种体载荷，用于模拟物体内部热源产生或消耗热量的过程，区别于作用在表面的热流密度边界条件，它直接作用于模型的体单元，适用于内热源驱动的温度场分析。 体积热流密度，SI单位为 W/m3（正号表示生成热量，负号表示消耗热量），用户输入界面如下： Image 生成的CalculiX求解器命令如下， *DFLUX ID_31,BF,10. 4. 对流 对流（Convection）边界条件又称膜（Film）边界条件。是用于模拟固体表面与相邻流体（气体或液体）之间对流换热的边界设置，是工程热分析中最常用的边界条件之一，广泛应用于电子设备散热、管道换热、汽车发动机冷却等场景。 对流换热系数是表征对流换热强度的核心参数，环境温度参数是远离固体表面的主流流体温度。前处理器中的输入界面如下， Image 生成的CalculiX求解器命令如下： ** BC Name: Convection *Film, op=New ID_29_S1, F1, 80, 123 ID_29_S2, F2, 80, 123 ID_29_S3, F3, 80, 123 ID_29_S4, F4, 80, 123 5. 热辐射边界条件 辐射边界条件是用于模拟物体表面与周围环境（或其他物体表面）之间通过热辐射方式传递热量的边界设置，其核心遵循斯蒂芬–玻尔兹曼定律，适用于真空、气体等介质的热交换场景（如航天器热控、高温设备散热）。辐射边界条件，是一种非线性边界条件，热流与温度的四次方成正比，求解时需要迭代计算。 用户需要输入的参数为发射率（Emissivity），介于0-1之间的数值，和环境温度。输入界面如下， Image 生成的热辐射边界条件命令如下， ** BC Name: Radiation *Radiate, op=New ID_30_S1, R1, 30, 0.9 ID_30_S2, R2, 30, 0.9 ID_30_S3, R3, 30, 0.9 ID_30_S4, R4, 30, 0.9 总结 有限元分析中，边界条件的类型较多，涉及了条件的参数本身，也涉及到所选择的单元，或节点。对于前处理器来说，生成各种边界条件的命令，是一项复杂的工程。 WELSIM已经可以很好的设置各种边界条件，能够快速生成CalculiX求解器所需要的输入文件，并可以直接用于计算。用户可以使用WELSIM作为CalculiX的前处理器。 CalculiX还有少量边界条件，如Constraint和Equation等, 由于常规分析中较少使用，本文未作讨论。 CalculiX与Abaqus的输入文件格式有着很大的相似性，因此，本文所描述的内容，也可以用于Abaqus求解器。 WelSim与作者不隶属于CalculiX, Abaqus。和CalculiX与Abaqus开发团队与机构没有直接关系。这里引用开源软件的名称和图片仅用作技术博客文章与软件使用的参考。

WELSIM 2026R1 Releases to Support Particle Generation and Enhance Electromagnetic Simulation

WELSIM, the general-purpose engineering simulation and analysis software, has released its latest version 2026R1 (Internal Version: 3.2). Compared to the previous version, 2026R1 comes with a host of new features and enhancements, enabling WELSIM to better support various types of engineering simulation CAE analyses, with notable improvements in particle and electromagnetic simulation capabilities. Support for Particle Generation The new version introduces particle generation from geometric models, along with the ability to adjust particle generation methods and particle density. This feature simplifies the setup of smoothed particle and molecular dynamics simulations in WELSIM. Press enter or click to view image in full size WELSIM is committed to supporting open data exchange, so users can also export particle data to VTK-format files for visualization or further analysis in other software tools. Enhanced Support for the OpenRadioss Solver Become a member The latest version improves compatibility with the OpenRadioss solver, allowing WELSIM to generate more types of OpenRadioss input files. The newly added input cards include: Fabric Geometric Property /PROP/TYPE16, Self-contact /INTER/TYPE24, and the SPH Cell Command SPHCEL. Meanwhile, a wide range of new material properties have been incorporated, such as LAW6 (HYD_VISC), LAW24 (CONC), LAW25 (COMPSH), LAW58 (FABR_A), LAW73, and LAW74. MatEditor, the standalone material editing software, has been updated with the same features accordingly. Press enter or click to view image in full size Enhanced Support for Palace The new version strengthens compatibility with Palace, the open-source electromagnetic solver. WELSIM now supports more Palace solver types and fully upgrades Palace and its dependent libraries to their latest versions. On the post-processing model, the Poynting Vector result type has been added. Press enter or click to view image in full size Other Enhancements and Upgrades With this release, WELSIM is more stable and user-friendly through optimizations to existing features. Disclaimer: WELSIM and its authors are not affiliated with OpenRadioss or Palace, nor do they have any connection with the respective development teams and organizations. The references to OpenRadioss and Palace herein are solely for informational purposes in technical blog posts and software usage guidance.

WELSIM发布2026R1版本，支持粒子生成，增强电磁计算

通用工程仿真分析软件WELSIM发布了最新的2026R1版本（内部版本号3.2）。相对于上一个版本，2026R1版本含有许多新的功能与增强，能够更好地支持各种类型的工程仿真CAE分析，尤其是更好的支持粒子与电磁计算。 Image 增强对粒子生成的支持 新版本增加了从几何体生成粒子的功能，同时可以调整粒子生成的方式和粒子密度。这个功能使得WELSIM在光滑粒子和分子动力学仿真上更加简单。 Image WELSIM一贯支持数据的交换与开放。用户还可以将粒子数据导出为VTK格式的文件，在其他软件上进行显示或分析。 增强对OpenRadioss求解器的支持 新版本支持增加对OpenRadioss的支持，使用WELSIM能生成更多类型的Radioss的输入文件。新增的Radioss输入卡片有：织物几何属性 /PROP/TYPE16，自接触 /INTER/TYPE24，SPH单元命令 SPHCEL。同时，新增了众多材料属性如：LAW6(HYD_VISC) ，LAW24 (CONC), LAW25 (COMPSH) ，LAW58 (FABR_A)，LAW73， LAW74等。独立材料编辑软件MatEditor，也同步增加了相同功能。 Image 增强对Palace的支持 新版本增强了对开源电磁求解器Palace的支持。支持了更多的Palace 求解类型。同时全面升级了Palace及其依赖库到最新的版本。在后处理上，增加了坡印廷矢量计算结果的显示。 Image 其他增强与升级 此外，新版本还对已有功能进行了优化与升级，使得WELSIM更稳定更加易于使用。 WelSim与作者不隶属于OpenRadioss和Palace, 和以上开发团队与机构没有关系。这里引用OpenRadioss, Palace 仅用作技术博客文章与软件使用的参考。

Compile the open-source electromagnetic simulation solver Palace

Palace is a parallel finite element solver developed by AWS Labs for full-wave 3D electromagnetic simulation, released under the Apache 2.0 open-source license. This solver supports full-wave frequency/time-domain simulation, eigenmode analysis, and parameter extraction for electrostatics/magnetostatics. It is compatible with platforms ranging from laptops to supercomputers and supports GPU acceleration, enabling the modeling of quantum computing hardware, RF/microwave devices, antennas, and other systems. The author previously published a brief article introducing how to compile Palace on Windows, titled Compile the electromagnetic simulation solver Palace in Windows. . Expanding upon that earlier work, this article details the compilation process more in-depth, with a particular focus on the compilation of dependent libraries. The dependent libraries discussed in this article are as follows: Hypre: A large-scale linear algebra matrix computation library, using version 2.30. MUMPS: An open-source software library for solving large-scale sparse linear systems of equations. ARPACK-NG: A library that supports complex linear matrix computations and is used for eigenvalue computation. Primarily written in Fortran, it can be compiled independently without relying on PETSc. GSLib: Used for interpolation computations in high-order spectral elements. libCEED: A linear algebra computation management framework that supports parallel computing across various CPUs, GPUs, and clusters. MFEM: A flexible, efficient, and scalable finite element discretization framework used as the core solving dependency library for Palace. System Environment and Compilers Operating System: Windows 11, 64-bit Compilers: Visual Studio 2022 Community (C++17), Intel Fortran Compiler 2022, and Intel MKL 2024. Palace Version: 0.15 Compiling Hypre Hypre is an open-source high-performance parallel linear solver library developed by Lawrence Livermore National Laboratory (LLNL) designed specifically for large-scale scientific computing and engineering simulation. Its core objective is to solve large-scale sparse linear systems of equations (in the form of Ax=b) efficiently, particularly strong at handling ultra-large-scale problems in distributed memory parallel environments. Hypre natively supports MPI parallelism, enabling it to utilize the full computing power of distributed memory architectures, such as supercomputers and clusters, and easily handle linear systems with millions or even billions of degrees of freedom. This stark advantage distinguishes Hypre from smaller solvers (e.g., Eigen). Compiling Hypre on Windows is simple. You can directly use the built-in CMake configuration file to generate Visual Studio project files, which will in turn produce the static library file. When configuring CMake, set the parameters as follows: HYPRE_WITH_MPI = ON HYPRE_USING_OPENMP = ON HYPRE_WITH_CUDA = OFF HYPRE_ENABLE_SYCL = OFF HYPRE_INSTALL_PREFIX=D:/WelSimLLC/CodeDV/3rdParty/hypre/myInstall CMAKE_INSTALL_PREFIX=D:/WelSimLLC/CodeDV/3rdParty/hypre/myInstall If the compilation proceeds smoothly, then the hypre.lib static library file will be generated. Compiling MUMPS To solve Wave Port problems, the Palace solver requires one of three direct linear algebra solvers: SuperLU, STRUMPACK, or MUMPS. MUMPS is selected as the dependent library given that the author is relatively familiar with it. MUMPS is an open-source parallel sparse direct solver library co-developed by French institutions including INRIA and CNRS. Its core objective is to solve large-scale sparse linear systems of equations (Ax=b) using the direct method. Unlike Hypre, which primarily relies on iterative methods, MUMPS computes without iteration by leveraging matrix decomposition (LU/Cholesky) via direct methods. This approach delivers higher precision and more robustness, making it particularly well-suited for solving sparse matrix problems involving asymmetric, symmetric positive definite, or indefinite matrices. MUMPS depends on BLAS/LAPACK (basic linear algebra libraries) and SCALAPACK (parallel linear algebra library). In this case, the dependency is the Intel MKL Library. The original version of MUMPS does not provide a compilation method for Windows. Instead, you can opt for the scivision/mumps version available on GitHub, which offers a CMake-based compilation approach for seamless implementation. The latest version 5.8.1 was downloaded for this example. First ensure that the environment variable contains the following configuration: MKLTOOL= C:\Program Files (x86)\Intel\oneAPI\mkl\latest Under CMake configuration, set the following parameters: MUMPS_openmp = ON MUMPS_parallel = OFF Once the Visual Studio project files are generated, compile the project to produce all static libraries of MUMPS. The generated static libraries are illustrated in the figure below. Press enter or click to view image in full size Compiling ARPACK For eigenvalue-related computations (e.g., eigenmode calculations), Palace requires a complex-number eigenvalue solver such as SLEPc or ARPACK. Since SLEPc depends on PETSc, which has a rather complex compilation on Windows, ARPACK is selected as the complex-number solver for Palace. ARPACK is an open-source library for large-scale eigenvalue and eigenvector computations, co-developed by institutions including Rice University and Sandia National Laboratories. Its core objective is to efficiently compute a subset of eigenvalues and their corresponding eigenvectors for large sparse matrices. Unlike libraries such as Hypre and MUMPS, which are designed to solve linear systems of equations (Ax=b), ARPACK focuses specifically on tackling eigenvalue problems (Ax=λx or Ax=λBx). Press enter or click to view image in full size Download ARPACK-NG from GitHub and use its built-in CMake build configuration. Set the following parameters in CMake: ICE = ON Eigen = ON Build_SHARED_LIBS = OFF Meanwhile, ensure that the environment variable contains the entry MKLTOOL= C:\Program Files (x86)\Intel\oneAPI\mkl\latest; this allows CMake to automatically locate the BLAS and LAPACK libraries within MKL. After generating the Visual Studio project files, you can directly compile the project to obtain the ARPACK static library. Compiling GSLib GSLIB is an open-source parallel sparse communication library designed to efficiently support sparse data gather/scatter communication and operators in parallel numerical simulation scenarios, such as the finite element method, spectral element method, and finite difference method. It greatly simplifies the development of parallel sparse communication and lowers the barrier to distributed memory programming. Its adaptive algorithms ensure efficient communication in large-scale parallel computing, satisfying ultra-large-scale simulation requirements. Press enter or click to view image in full size Since GSLIB does not come with CMake configuration files, it is unable to automatically generate Visual Studio project files on Windows, requiring you to create the Visual Studio project files manually. Create a new static library project named gslib-palace, and add all the source files and header files to this project. In the preprocessor definitions of the Visual Studio project, add GSLIB_USE_CBLAS and GSLIB_USE_MKL. During compilation, when the following error messages appear: unary minus operator applied to unsigned type, result still unsigned You will need to modify the code, for example, changing p[map->n].el = -(uint)1; to p[map->n].el = 0-(uint)1; This allows the compilation to succeed. Become a member When the following error message appears: potentially uninitialized local pointer variable 'h' used You will need to modify the code and initialize the pointer with a NULL value. If you receive a prompt indicating that the fgs_free function has already been defined redundantly, simply comment out the fgs_setup, fgs_free, and fgs_unique functions. Compiling libCEEM libCEED is an efficient and scalable discretization library that focuses on high-performance computational discretization based on finite element and spectral element methods. It provides a lightweight algebraic interface for linear and nonlinear operators as well as preconditioners; it supports runtime selection and optimized implementations across a variety of computing devices such as CPUs and GPUs. Compiling libCEED is the most challenging part of the entire task because libCEED does not provide CMake configuration files and contains numerous options. Coupled with the fact that Visual Studio does not support C99 syntax, extensive manual modifications to the source files are required, especially modifications related to dynamic arrays. For example, change the following dynamic array: CeedInt num_points[num_elem]; to: CeedInt *num_points = (CeedInt *)malloc(num_elem * sizeof(CeedInt)); After usage, you also need to free the array from memory at the end using the free function: free(num_points);. Certain strings containing file paths also require modification. For instance, change: *(relative_file_path) = strstr(absolute_file_path, "ceed/jit-source"); to: #ifdef _WIN32 *(relative_file_path) = strstr(absolute_file_path, "ceed\\jit-source"); #else *(relative_file_path) = strstr(absolute_file_path, "ceed/jit-source"); #endif Since Visual Studio’s support for weak functions differs from that of Linux, the weak functions in the source code need to be modified as follows: Original code: #define CEED_BACKEND(name, num_prefixes, ...) \ CEED_INTERN int name(void) __attribute__((weak)); \ int name(void) { return CeedRegister_Weak(__func__, num_prefixes, __VA_ARGS__); } #include "ceed-backend-list.h" Modified code: #ifdef _WIN32 #define CEED_BACKEND(name, num_prefixes, ...) \ CEED_INTERN int name(void) __declspec(selectany); \ int name(void) { return CeedRegister_Weak(__func__, num_prefixes, __VA_ARGS__); } #else #define CEED_BACKEND(name, num_prefixes, ...) \ CEED_INTERN int name(void) __attribute__((weak)); \ int name(void) { return CeedRegister_Weak(__func__, num_prefixes, __VA_ARGS__); } #include "ceed-backend-list.h" #endif Visual Studio does not support the __restrict__ keyword and requires it to be changed to __restrict. Other minor modifications need to be made, such as replacing popen with _popen and pclose with _pclose. Additionally, there is an environment variable setting: setenv("RUST_TOOLCHAIN", "nightly", 0), which needs to be modified as follows: char env_str[1024]; snprintf(env_str, sizeof(env_str), "%s=%s", "nightly", 0); _putenv(env_str); To reduce complexity, this compilation supports only CPU parallel computing and does not require GPU parallel support. Therefore, in the ceed-backend-list.h file, only the following functionalities need to be enabled: CEED_BACKEND(CeedRegister_Opt_Blocked, 1, "/cpu/self/opt/blocked") CEED_BACKEND(CeedRegister_Opt_Serial, 1, "/cpu/self/opt/serial") CEED_BACKEND(CeedRegister_Ref, 1, "/cpu/self/ref/serial") CEED_BACKEND(CeedRegister_Ref_Blocked, 1, "/cpu/self/ref/blocked") CEED_BACKEND(CeedRegister_Xsmm_Blocked, 1, "/cpu/self/xsmm/blocked") CEED_BACKEND(CeedRegister_Xsmm_Serial, 1, "/cpu/self/xsmm/serial") All other computing functionalities can be commented out. Upon successful compilation, the static library for libCEED will be generated. Compiling MFEM MFEM is the core solver of Palace, an open-source high-order finite element library led by the Lawrence Livermore National Laboratory (LLNL). Designed specifically for large-scale scientific computing and engineering simulation, its core objective is to provide a flexible, efficient, and scalable finite element discretization framework. MFEM also serves as a mainstream open-source finite element tool for computational fluid dynamics, electromagnetic simulation, structural mechanics, nuclear physics, and other fields. Press enter or click to view image in full size Compiling MFEM is relatively straightforward. You can generate Visual Studio project files via CMake for compilation. Set the following options to ON in CMake: MFEM_USE_MPI = ON MFEM_USE_ZLIB = ON MFEM_USE_METIS = ON MFEM_USE_LAPACK = ON MFEM_USE_LEGACY_OPENMP = ON MFEM_THREAD_SAFE = ON MFEM_USE_MUMPS = ON Generate the Visual Studio project files and compile them directly to obtain the MFEM static library. Compiling Palace At this point, all complex dependent libraries have been successfully compiled. Proceed to the final step of this project: use the built-in CMake configuration file to generate Visual Studio project files and compile Palace. Note that you should select the CMakeLists.txt file located in the palace subdirectory instead of the one in the root directory, as the latter is used for generating the SuperBuild. Set the following parameters in CMake: PALACE_WITH_ARPACK = ON PALACE_WITH_OPENMP = ON PALACE_WITH_SLEPC = OFF MPI_Fortran_WORKS = C:/Program Files (x86)/Microsoft SDKs/MPI/Include/x64 Configure the following paths: nlohmann_json_DIR = D:/WelSimLLC/CodeDV/3rdParty/nlohmann_json/json-3.12.0/myInstall/share/cmake/nlohmann_json fmt_DIR = D:\WelSimLLC\CodeDV\3rdParty\fmt\fmt-11.0.2\myInstall\lib\cmake\fmt scn_DIR = D:\WelSimLLC\CodeDV\libPack\lib\scn\cmake\scn MFEM_DIR = D:/WelSimLLC/CodeDV/libPack/lib/palace/mfem-cmake arpackng_DIR = D:\WelSimLLC\CodeDV\libPack\lib\palace\arpack-cmake\arpackng Open the generated Visual Studio project, and you will see a static library project named libpalace and an executable project named palace. If the compiler reports missing header files or the linker indicates unresolved functions during compilation, you can directly add the relevant paths and static libraries to the projects. After completing the compilation, you can run a test on palace.exe. Note that you need to place all dynamic dependent libraries (such as the MKL dynamic libraries) in the same directory as palace.exe before running it. Run a simple test case. The output shown below indicates that the compilation has been successful basically. Press enter or click to view image in full size Conclusion As open-source electromagnetic field simulation solvers are presently scarce, Palace stands out as the most feature-rich one available. As a result, Palace relies on a large number of dependent libraries, and its version is still under continuous iterations. Some of these dependencies only support Linux builds, which adds to the complexity of compiling Palace on Windows. This article documents the complete compilation process of Palace and provides a resource for the open-source community on how to build Palace for Windows. There are also other dependent libraries that are relatively easy to compile, such as Eigen, fmt, Metis, and nlohmann/json. Experienced developers can often independently build these libraries on Windows, so it is not discussed in this article. The author has included the compiled palace.exe executable file in the WELSIM installation package, so you are able to directly obtain the Windows version of palace.exe from the WELSIM installer without having to compile it yourself. This article describes the configuration of Palace for CPU multi-core parallel computing. The compilation methods for GPU-accelerated versions, such as CUDA, HIP, or SYCL, will be discussed in future articles. WelSim and the author are not affiliated with Palace, Hypre, MUMPS, ARPACK, GSLIB, libCEED, or MFEM, nor do they have any direct connections with the development teams or institutions behind these projects. The references to the names and images of these open-source software in this article are solely for the purpose of technical blog documentation and software usage guidance.