Yukun’s WARP3D(Version 17.8.7) Reference Note¶
Version 17.9.2: didriv.f: Domain Integral DRIVer¶
Description¶
- domain driver:
drive the computation of j-integral values using the domain integral technique.
or:
drive the computation of mixed-mode stress intensity factors and t-stress using the interaction integral.
didriv.f Declaration Tree¶
didriv.f Calls Map¶
Procedure¶
1. check validity of domain before starting computations.¶
在开始计算前检查积分域
2. perform exhaustive check on consistency of 1) number of front nodes, 2) front interpolation order, and 3) domain type.¶
对一致性进行彻底检查：1）裂尖节点数，2）裂尖插值阶数，3）积分域类型。
call dickdm:
c ***************************************************************
c * *
c * domain_check - exhaustive testing of domain defintiion for *
c * consistency *
c * *
c ***************************************************************
3. allocate space for a q-value at each structure node. allocate space for expanded lists of coincident front nodes at each front location.¶
为q_values（q值函数）和expanded_front_nodes（裂尖重合节点）分配内存
4. output header information for domain¶
输出积分域的头信息
call diheadr:
c ***************************************************************
c * *
c * domain_header - output info at start of domain computations *
c * *
c ***************************************************************
5. when q-values are given by the user, expand the compressed list of q-values into a list of length number of structure nodes.¶
当用户给出q值时，将压缩的q值列表展开为长度为结构节点数的列表。
6. set up the domain. for automatic domains, set q-values at front nodes. then find coincident front nodes and set their q-values.¶
设置q值函数。对于自动q值积分域，在裂尖节点设置q值。然后找到重合的裂尖节点并设置它们的q值。
call distup.f:
c ***************************************************************
c * *
c * set_up_domain -- set up the domain for processing *
c * find coincident front nodes and set their *
c * q-values *
c * *
c ***************************************************************
Output:
j_data.q_valuesREAL (:) ALLOCATABLE SAVE
7a. create vector with list of elements on the crack front for this domain¶
创建积分域中位于裂尖的单元列表的数组
call di_cf_elem:
c **********************************************************************
c * *
c * di_cf_elem - create a logical vector whose entries are .true. for *
c * elements incident on the crack tip, and .false. for *
c * those that are not. dicmj will use this info to set a *
c * flag for each element that is analyzed by dielem. if *
c * a user includes the domain integral command *
c * 'omit crack front elements for fgms yes', the flag *
c * will cause terms7 and 8 to be set to zero. *
c * *
c **********************************************************************
Output:
j_data.crack_front_elemLOGICAL (:) ALLOCATABLE SAVE
7b. use crack front element to set if J computations will use small or large strain formulation. all elements in domain must match be the same formulation.
依据裂尖单元设置计算J积分使用小应变或大应变形式。积分域中的所有单元必须使用相同的变形形式。
8. at point on front where integral is being computed, build the global->crack rotation matrix. gather coordinates and displacements of crack-front nodes, and rotate them to local crack-front system.¶
在将要计算积分的裂尖点处，建立全局至裂尖旋转矩阵。收集裂尖节点的坐标和位移，并将它们旋转到裂尖局部坐标系。
call dimrot.f:
c **********************************************************************
c * *
c * dimrot - compute the 3x3 global -> crack front local rotation *
c * *
c **********************************************************************
Output:
j_data.domain_originINTEGERj_data.domain_rot(3,3)DOUBLE PRECISION (3,3)
8c. calculate strain e33 at node at domain origin. this is for T-stress calculations using the interaction integral
计算域起源节点处的应变e33。这是用于使用相互作用积分的T-应力计算。 仅线性弹性分析
call di_calc_e33:
c *******************************************************************
c * *
c * calculate strain e33 at domain origin for T-stress calcs. *
c * calculate strain e33 as the difference between the *
c * deformed and undeformed crack-front lengths delta_L / L *
c * *
c *******************************************************************
8c. calculate properties of a curve passing through the front nodes. these will be used to compute distance ‘r’ from integration points to a curved crack front.
计算通过前节点的曲线的属性。这将用于计算从积分点到弯曲裂纹前沿的距离‘r’。交互积分的线弹性分析
call di_calc_curvature from
c *******************************************************************
c * *
c * calculate coefficients of curve described by crack front *
c * nodes. *
c * *
c *******************************************************************
9. compute area under the q-function over that part of crack front for this domain. the area must be >0 else fatal error in domain (user forgot to set q-values on front nodes)¶
计算在该积分域中裂纹前沿的q函数下的面积。面积必须大于0
call difrar.f:
c **********************************************************************
c * *
c * di_front_q_area - compute area under q-function along front for *
c * this domain *
c * *
c **********************************************************************
10. formerly set logical flags to indicate if the nodal velocities and accelerations are all zero for this load step.¶
no longer used
以前设置逻辑标志以指示此加载步骤的节点速度和加速度是否都为零。
11.¶
Need to do this when temperatures have been specified in model. Possibly leading to computation of J_6
Not needed if the user has input initial stresses or set an initial state for J-processing to accommodate prior thermo-mechanical processing. In such cases those effects, FGMs, thermal, … will be included thru J_7, J_8
Build the node average value of thermal expansion coefficient. for temperature-dependent material properties, also build the node average value of young’s modulus and poisson’s ratio. For temperature-independent material properties, values of e and nu are obtained within dicmj.f
nodal properties are needed for domain integral computations to compute spatial derivatives within the domain.
在模型中指定温度时需要执行此操作。可能导致J_6的计算
如果用户输入初始应力或设置J处理的初始状态以适应先前的热机械处理，则不执行此操作。在这种情况下，这些效果（FGM，热，等）将包括在J_7，J_8中
建立热膨胀系数的节点平均值。对于温度依赖的材料特性，还建立杨氏模量和泊松比的节点平均值。对于与温度无关的材料属性，e和nu的值在dicmj.f中获得。
区域积分的计算需要节点属性来计算域内的空间导数。
call di_node_props_setup:
c **********************************************************************
c * *
c * di_node_props_setup - obtain alpha values at nodes. for *
c * temperature-dependent properties, also *
c * compute e and nu values at nodes. this *
c * routine replaces di_expan_coeff_setup, *
c * and the routine it calls, di_node_props, *
c * replaces di_node_expan_coeff. *
c **********************************************************************
12.¶
Needed when FGM nodal properties are actually used for elements in model, when the user has defined an initial state, or temperature dependent stress-strain curves are used (effectively makes the material an FGM), or input initial stresses.
Needed for I integrals with temperatures in model.
Even if not needed, the code builds small, dummy arrays to satisfy checks.
Build the nodal averages of strain energy density (W) and nodal values of displacement gradient (3x3) at t_n relative to coordinates at t=0. These terms are used to calculate the derivative of W wrt crack local X and derivatives of the gradients wrt crack local X – both at integration points for J(7) and J(8).
当1）FGM节点属性实际用于模型中的单元，2）用户定义了初始状态，3）使用温度相关的应力-应变曲线（有效地使材料成为FGM），4）输入初始应力时，需要需要执行此操作。
模型中有温度时需要执行此操作。
即使不需要执行，代码也会构建小的虚拟数组以满足检查要求。
在t_n处相对于坐标t=0处建立应变能密度（W）的节点平均值和位移梯度（3×3）的节点值。这些项用于计算应变能密度相对于裂尖局部坐标X的导数和位移梯度相对于裂尖局部坐标X的导数。这两个量在积分点处的值都将用在J（7）和J（8）中。
call di_setup_J7_J8:
c **********************************************************************
c * *
c * allocate data structures and drive computations *
c * to enable subsequent computations of J(7) and J(8) terms. *
c * *
c * if computations not actually required, create small, dummy arrays *
c * to simplify checking *
c * *
c * get a stress work value (W) at all model nodes by extrapoaltion *
c * from element integration point values, then node averaging *
c * *
c * also get the displacement gradiaents partial(u_i/x_j) (3x3) *
c * at model nodes by extrapolation from element integration points *
c * *
c **********************************************************************
12b. at point on front where integral is being computed, collect young’s modulus and poisson’s ratio. this assumes that all elements connected to this crack-front node have identical, homogeneous material properties, or that fgm material properties have been assigned to the model. for homogeneous material, “props” contains material data. for fgms, read data from “fgm_node_values.” for temperature-dependent properties, segmental data arrays contain the properties.
提取在计算积分的裂尖点处的杨氏模量和泊松比。这里假设连接到此裂尖节点的所有单元具有相同的均匀材料属性，或者已将梯度材料属性分配给模型。对于均匀材料，变量“props”储存了材料数据。对于梯度材料，从变量“fgm_node_values”读取数据。对于与温度相关的材料属性，分段数据数组包含属性。
12c. at point on front where integral is being computed, props(7,elemno) and props(8,elemno) may be zero if the element has been killed. in this case, use the first nonzero values of e and nu.
在计算积分的裂尖点上，如果单元被杀死，则props（7，elemno）和props（8，elemno）可能为零。 在这种情况下，使用第一个e和nu的非零值。
12d. output info for domain about temperature at crack front and alpha values. Also E, nu at front
输出积分域关于裂尖的温度，热膨胀系数，杨氏模量和泊松比的信息。
13a. init various tracking values for for J and I¶
初始化用于J积分和I积分的数组
13b. separate drivers for comput values on 1 domain for user-defined q-values, and for defining-computing over automatic domains
分离在一个积分域上使用用户定义的q值进行J、I积分，或在一个积分域上使用自动q值进行J、I积分
14. user wants automatic construction of domains.¶
14a. get last ring at which output will be printed. domains are always generated starting at ring 1 but j-values may not be computed for every domain.
14b. allocate arrays needed to support construction/definition of the domains. for type 4, we need only a nodal bit map. for types 1-3, we need two sets of nodal bit maps. each set has 3 maps of length to record all structure nodes. element list stored in the common vector.
14c. set up to accumulate statistics for computed domain values.
14d. loop over all domains. construct definition of the domain (q-values, element list). call driver to actually calculate value for domain.
call diexp4.f:
c ***************************************************************
c * *
c * domain expand 4 - expand type 4 automatic domain *
c * domain expand 13 - expand type 1-3 automatic domain *
c * *
c ***************************************************************
call dicmj.f:
c ***************************************************************
c * *
c * domain_compute - drive execution of element routine to *
c * compute j and i-integrals for a single *
c * domain *
c * *
c ***************************************************************
15a. write j-integral and i-integral data to standard output¶
将j积分和i积分数据写入标准输出
15b. write j-integral and i-integral data to packets¶
将j积分和i积分数据写入数据包
16. release arrays used for both user defined and automatic domains. call routines to get releases on other than rank 0 for MPI¶
释放用于用户定义积分域和自动域的数组
Version 17.9.2: dicmj.f & dielem.f Declaration Tree and Calls Map¶
dicmj.f Declaration Tree¶
dicmj.f Calls Map¶
dielem.f Calls Map¶
didriv.f: Domain Integral DRIVer¶
Description¶
- domain driver:
drive the computation of j-integral values using the domain integral technique.
or:
drive the computation of mixed-mode stress intensity factors and t-stress using the interaction integral.
Calling Tree¶
c ***************************************************************
c * *
c * didriv.f *
c * -dickdm *
c * -diheadr *
c * -distup.f *
c * -di_cf_elem *
c * -dimrot.f *
c * -dimrot_coord_displ *
c * -di_trans_nodvals *
c * -dimrott *
c * -di1dsf.f *
c * -difrts *
c * -di_calc_e33 *
c * -di_calc_curvature *
c * -di_calc_coefficients *
c * -difrar.f (di_front_q_area) *
c * -di1dsf.f *
c * -di_node_props_setup *
c * -di_node_props *
c * -di_fgm_alphas *
c * -di_constant_alphas *
c * -di_seg_alpha_e_nu *
c * -di_fgm_setup *
c * -di_nod_vals *
c * -di_extrap_to_nodes *
c * -ndpts1.f *
c * -oulg1.f (oulgf) *
c * -diexp4.f *
c * -diexp4 *
c * -diadit *
c * -dibmck *
c * -diexp13 *
c * -diadit *
c * -dibmck *
c * -digete *
c * -dicmj.f *
c * -di_trans_nodvals *
c * -vec_ops (zero_vector.f) *
c * -dielem *
c * -dielem_a.f *
c * -dielem_b.f *
c * -dielem_c.f *
c * -di_write_std_out *
c * -di_write_packets *
c * *
c ***************************************************************
Procedure¶
- check validity of domain before starting computations.
- perform exhaustive check on consistency of 1) number of front nodes, 2) front interpolation order, and 3) domain type.
call dickdm:
c ***************************************************************
c * *
c * domain_check - exhaustive testing of domain defintiion for *
c * consistency *
c * *
c ***************************************************************
- allocate space for a q-value at each structure node. allocate space for expanded lists of coincident front nodes at each front location.
- output header information for domain
call diheadr:
c ***************************************************************
c * *
c * domain_header - output info at start of domain computations *
c * *
c ***************************************************************
- when q-values are given by the user, expand the compressed list of q-values into a list of length number of structure nodes.
- set up the domain. for automatic domains, set q-values at front nodes. then find coincident front nodes and set their q-values.
call distup.f:
c ***************************************************************
c * *
c * set_up_domain -- set up the domain for processing *
c * find coincident front nodes and set their *
c * q-values *
c * *
c ***************************************************************
Output:
j_data.q_valuesREAL (:) ALLOCATABLE SAVE
- allocate a vector of logicals and assign .true. for each element connected to a crack front node.
call di_cf_elem:
c **********************************************************************
c * *
c * di_cf_elem - create a logical vector whose entries are .true. for *
c * elements incident on the crack tip, and .false. for *
c * those that are not. dicmj will use this info to set a *
c * flag for each element that is analyzed by dielem. if *
c * a user includes the domain integral command *
c * 'omit crack front elements for fgms yes', the flag *
c * will cause terms7 and 8 to be set to zero. *
c * *
c **********************************************************************
Output:
j_data.crack_front_elemLOGICAL (:) ALLOCATABLE SAVE
- at point on front where integral is being computed, build the global->crack rotation matrix. gather coordinates and displacements of crack-front nodes, and rotate them to local crack-front system.
call dimrot.f:
c **********************************************************************
c * *
c * dimrot - compute the 3x3 global -> crack front local rotation *
c * *
c **********************************************************************
Output:
j_data.domain_originINTEGERj_data.domain_rot(3,3)DOUBLE PRECISION (3,3)
8c. calculate strain e33 at node at domain origin. this is for T-stress calculations using the interaction integral
call di_calc_e33:
c *******************************************************************
c * *
c * calculate strain e33 at domain origin for T-stress calcs. *
c * calculate strain e33 as the difference between the *
c * deformed and undeformed crack-front lengths delta_L / L *
c * *
c *******************************************************************
8c. calculate properties of a curve passing through the front nodes. these will be used to compute distance ‘r’ from integration points to a curved crack front.
call di_calc_curvature from
c *******************************************************************
c * *
c * calculate coefficients of curve described by crack front *
c * nodes. *
c * *
c *******************************************************************
- compute area under the q-function over that part of crack front for this domain. the area must be >0 else fatal error in domain (user forgot to set q-values on front nodes)
call difrar.f:
c **********************************************************************
c * *
c * di_front_q_area - compute area under q-function along front for *
c * this domain *
c * *
c **********************************************************************
- set logical flags to indicate if the nodal velocities and accelerations are all zero for this load step. if so, some later computations can be skipped.
- Build the node average value of thermal expansion coefficient. for temperature-dependent material properties, also build the node average value of young’s modulus and poisson’s ratio. for temperature-independent material properties, values of e and nu are obtained within dicmj.f. nodal properties are needed for domain integral computations to compute spatial derivatives within the domain.
call di_node_props_setup:
c **********************************************************************
c * *
c * di_node_props_setup - obtain alpha values at nodes. for *
c * temperature-dependent properties, also *
c * compute e and nu values at nodes. this *
c * routine replaces di_expan_coeff_setup, *
c * and the routine it calls, di_node_props, *
c * replaces di_node_expan_coeff. *
c **********************************************************************
- Build the nodal averages of strain energy density (stress work density) and strains. These terms are used to calculate the derivative of the strain energy density, which appears in the domain integral when material properties vary spatially (e.g. fgms). The nodal values calculated in di_fgm_setup will be used to compute their spatial derivatives at integration points.
call di_fgm_setup:
c **********************************************************************
c * *
c * di_fgm_setup - allocate data structures for two terms used *
c * in the calculation of the derivative of the stress work density. *
c * these are: nodal values of stress work density and strain. *
c * *
c **********************************************************************
12b. at point on front where integral is being computed, collect young’s modulus and poisson’s ratio. this assumes that all elements connected to this crack-front node have identical, homogeneous material properties, or that fgm material properties have been assigned to the model. for homogeneous material, “props” contains material data. for fgms, read data from “fgm_node_values.” for temperature-dependent properties, segmental data arrays contain the properties.
- if q-values given by user, we compute domain integral right now. otherwise, we set up a loop to generate q-values for automatic domains and their computation.
- user wants automatic construction of domains.
14a. get last ring at which output will be printed. domains are always generated starting at ring 1 but j-values may not be computed for every domain.
14b. allocate arrays needed to support construction/definition of the domains. for type 4, we need only a nodal bit map. for types 1-3, we need two sets of nodal bit maps. each set has 3 maps of length to record all structure nodes. element list stored in the common vector.
14c. set up to accumulate statistics for computed domain values.
14d. loop over all domains. construct definition of the domain (q-values, element list). call driver to actually calculate value for domain.
call diexp4.f:
c ***************************************************************
c * *
c * domain expand 4 - expand type 4 automatic domain *
c * domain expand 13 - expand type 1-3 automatic domain *
c * *
c ***************************************************************
call dicmj.f:
c ***************************************************************
c * *
c * domain_compute - drive execution of element routine to *
c * compute j and i-integrals for a single *
c * domain *
c * *
c ***************************************************************
14e. release allocatable arrays for automatic domains
15a. write j-integral and i-integral data to standard output
15b. write j-integral and i-integral data to packets
- release arrays used for both user defined and automatic domains
distup.f: Domain Integral SeT UP¶
Description¶
call distup( c, crdmap, nonode, debug_driver, out )
set_up_domain:
set up the domain for processing find coincident front nodes and set their q-values
Procedure¶
1 determine if the user has specified the crack front nodes using node sets. this is done when the model is not collapsed at the ct. For a user defined crack tip, load expanded_front_nodes with node_set specified by the user. check to be sure that each front position has the same number of crack tip nodes.
确定用户是否使用节点集指定了裂尖节点。当模型未在裂尖处塌缩时完成此操作。对于用户定义的裂尖，使用用户指定的节点集加载expanded_front_nodes数组。检查以确保每个裂尖位置具有相同数量的裂尖节点
- for automatic domain definitions, set the q-values at user specified front nodes based on the domain type and the user specified front interpolation order.
确定用户是否使用节点集指定了裂尖节点。当模型未在裂尖处塌缩时完成此操作。对于用户定义的裂尖，使用用户指定的节点集加载expanded_front_nodes数组。检查以确保每个裂尖位置具有相同数量的裂尖节点
- if linear interpolation of q-values over the domain is on, we also need to interpolate q-values along front for quadratic elements.
对于自动积分域定义，根据区域类型和用户指定的插值阶数在用户指定的裂尖点处设置q值。
- (skip this section when the user has defined the ct.) look at first 1 or 2 nodes specified to be on the crack front. compute a distance to be used for locating all other nodes coincident with each specified crack front node. set the “box” size for proximity checking as a fraction of the distance measure.
- loop over each specified front node. locate all other nodes in model that are coincident with the front node. (use the box test). set the q-value of each coincident node to be same as user specified value for the specified front node. dump verification list if user requested it.
dimrot.f: Domain Integral ROTation¶
Description¶
call dimrot( c, crdmap, u, dstmap, debug_driver, out )
compute the 3x3 global -> crack front local rotation
Output:
j_data.domain_originINTEGER
j_data.domain_rot(3,3)DOUBLE PRECISION (3,3)
Calling Tree¶
c ***************************************************************
c * *
c * dimrot.f *
c * -dimrot_coord_displ *
c * -di_trans_nodvals *
c * -dimrott *
c * -di1dsf.f *
c * -difrts *
c * *
c ***************************************************************
Procedure¶
- gather coordinates and displacements of crack-front nodes
2.1 if only one node on the crack front. we have the simple 2-D case. rotation matrix determined only by direction cosines. global-z and crack-z are assumed to coincide.
call dimrot_coord_displ
2.2 if more than one node on front specified. we have a 3-D situation. get the position of the front-node to be used as the domain origin. for domain types ‘a’, ‘b’ and ‘c’, (types 1, 2 and 3, respectively), take the domain origin as the location where the crack-front advance is maximum (where the q-value equals 1.0). for domain type ‘d’ (uniform crack advance) take the domain origin as the first front node.
- based on the number of front nodes, interpolation order, and q-function type construct a unit tangent vector to the front that lies in the crack plane (crack-z direction). crack-y is given by direction cosines of normal to crack plane. construct crack-x and final 3x3 structure-to-crack rotation by cross-product.
3a. for q-function type ‘d’ (=4), we just use the first two nodes on the front to find the tangent vector. this case is for a uniform crack advance along the full front.
alternatively, the user may have defined components of the crack front tangent vector rather than using the computed value. this feature helps, for example, when the mesh has the crack front intersecting a symmetry plane at other than 90-degrees. A correct user-defined tangent vector restores domain independent J-values.
call dimrott
- get crack front (unit) tangent vector (crack-z). compute crack-x unit vector from z and y. fill in global-to-crack rotation matrix.
call dimrot_coord_displ
Call dimrot_coord_displ¶
rotate coordinates of crack-front nodes to local system.
if necessary, transform nodal displacement values from constraint-compatible coordinates to global coordinates. then rotate displacements from global to local coordinates
call di_trans_nodvals
subroutine to transform a 3x1 vector of constraint-compatible-coordinate-system values to global-coordinate-system values. the routine premultiplies the incoming constraint-compatible values by the transpose of the stored global-to-constraint coordinate system rotation matrix.
Call dimrott¶
unit tangent vector to front directed along crack local-Z direction
Output:
tvec(3)
c compute or load unit tangent vector to the crack
c front according to the situation indicated by the
c case number or the user-defined flag.
c for cases 2,5,6 we average two tangents
c computed for adjacent line segments. for cases 1,3,4 we
c compute a single tangent at start or end of line segment.
c fnodes is list of structure nodes defining the segment(s).
c
c tangent_vector_defined = .T. --> the user specified the
c tangent vector to use.
c components are in
c crack_front_tangent
c skip computations.
c
c
c case no. situation
c -------- ---------
c 1 linear; 2-nodes on front
c 2 linear; 3-nodes on front
c 3 quadratic; 3-nodes on front
c 4 cubic; 4-nodes on front
c 5 quadratic; 5-nodes on front
c 6 cubic; 7-nodes on front
c
c fnctyp: 1-4, used for cases 3,4 to determine if
c tangent is needed for first or last node
c on front. = 1 (first node), =3 (last node).
c
c status: returned value = 0 (ok) = 1 (bad data)
c
c case 1,3,4: 2, 3 or four nodes on front but only a single
c element (linear, quadratic or cubic). compute
c tangent at first or last node (2,3 or 4). use
c utility routine to get derivatives of the 1-D
c isoparametric shape functions.
c
c case 2,5,6: two linear, quadratic or cubic segments
c connecting 3, 5 or 6 front nodes.
c compute average tangent at center (common)
c node. compute two unit tangents, average
c components, then restore unit length.
Call di1dsf.f¶
calculates shape function values and derivatives for 1-dimension
Output:
sf(3)- shape function values
dsf(3)- derivatives of shape function
di_calc_curvature: Domain Integral CALculate CURVATURE¶
Description¶
- Domain Integral CALculate CURVATURE:
- calculate coefficients of curve described by crack front nodes.
c determine if crack front is curved.
c find coefficients of a quadratic expression defined by
c three crack front nodes: x = az^2 + bz + c
c if 'a' equals zero +/- tolerance, consider the front
c nodes to be colinear, and exit.
c
c for a circular crack, find center and radius of a
c circle defined by the first three crack-front nodes.
Output:
j_data.crack_curvatureDOUBLE PRECISION (7)
c contents of array 'crack_curvature':
c
c crack_curvature(1) = 0 for straight front
c = 1 for curved front
c crack_curvature(2) = 0 for unknown curvature
c = 1 for for a circular curve
c crack_curvature(3) = local x coordinate for circle
c = coefficient 'a' for polynomial
c crack_curvature(4) = local z coordinate for circle
c = coefficient 'b' for polynomial
c crack_curvature(5) = circle radius for circle
c = coefficient 'c' for polynomial
c crack_curvature(6) = coefficient 'd' for polynomial
c crack_curvature(7) = coefficient 'e' for polynomial
Calling Tree¶
c ***************************************************************
c * *
c * di_calc_curvature *
c * -di_calc_coefficients *
c * *
c ***************************************************************
curvature:
- q = ( x3 - x2 ) / ( ( z3 - z2 ) * ( z3 - z1 ) )
- & - ( x1 - x2 ) / ( ( z1 - z2 ) * ( z3 - z1 ) )
circular curve:
- x_center = ( ( x1*x1 + z1*z1 ) * ( z3 - z2 )
- & + ( x2*x2 + z2*z2 ) * ( z1 - z3 ) & + ( x3*x3 + z3*z3 ) * ( z2 - z1 ) ) & / ( two * ( x1 * ( z3 - z2 ) & + x2 * ( z1 - z3 ) & + x3 * ( z2 - z1 ) ) )
- z_center = ( ( x1*x1 + z1*z1 ) * ( x3 - x2 )
- & + ( x2*x2 + z2*z2 ) * ( x1 - x3 ) & + ( x3*x3 + z3*z3 ) * ( x2 - x1 ) ) & / ( two * ( z1 * ( x3 - x2 ) & + z2 * ( x1 - x3 ) & + z3 * ( x2 - x1 ) ) )
- circle_radius = ( (x1 - x_center)**two
- & + (z1 - z_center)**two )**half
3-node curve:
- a = ( x3 - x2 ) / ( ( z3 - z2 ) * ( z3 - z1 ) )
- & - ( x1 - x2 ) / ( ( z1 - z2 ) * ( z3 - z1 ) )
b = ( x1 - x2 + a * ( z2*z2 - z1*z1 ) ) / ( z1 - z2 )
c = x1 - a * z1*z1 - b * z1
5-node curve:
call least-squares regression subroutinedi_calc_coefficientsto compute coefficients of polynomial.
Call di_calc_coefficients¶
subroutine to compute coefficients of a fourth-order polynomial by a linear least-squares regression of the five data points of the five crack-front nodes.
or:
solve a linear set of simultaneous equations using gauss elimination with partial pivoting. one rhs is permitted. the coefficient matrix may be non-symmetric. the coefficient matrix is replaced by it’s triangulated form.
solves: [a] {x} = {b}:
c * a -- coefficient matrix of *
c * x -- vector to receive computed solution *
c * b -- vector containing the right hand side *
x = [e d c b a]
difrar.f: Domain Integral FRont ARea¶
Description¶
di_front_q_area:
compute area under q-function along front for this domain
Output:
j_data.front_q_areaDOUBLE PRECISION
j_data.front_lengthDOUBLE PRECISION
di_node_props_setup: Domain Integral NODE PROPertieS SETUP¶
Description¶
call di_node_props_setup (1, nonode, nelblk )
di_node_props_setup:
obtain alpha values at nodes. for temperature-dependent properties, also compute e and nu values at nodes. this routine replaces di_expan_coeff_setup, and the routine it calls, di_node_props, replaces di_node_expan_coeff.
Output:
j_data.count_alphaINTEGER (nonode) ALLOCATABLE SAVE
j_data.snode_alpha_ijREAL (nonode,6) ALLOCATABLE SAVE
j_data.seg_snode_eREAL (nonode) ALLOCATABLE SAVE
j_data.seg_snode_nuREAL (nonode) ALLOCATABLE SAVE
j_data.block_seg_curvesLOGICAL (nellk) ALLOCATABLE SAVE
j_data.process_temperaturesLOGICAL
Calling Tree¶
c ***************************************************************
c * *
c * -di_node_props_setup *
c * -di_node_props *
c * -di_fgm_alphas *
c * -di_constant_alphas *
c * -di_seg_alpha_e_nu *
c * *
c ***************************************************************
Call di_node_props¶
di_node_props:
build average nodal values of material properties at each node in the model. all values are single-precision reals because double-precision accuracy is unnecessary.
loop through element blocks. check type of property assignment in each block.
property assignment: element nodes (fgm) = 1
constant in element = 2
temperature dependent = 3
- check for fgm alphas. the identifier ‘fgm_mark’ for fgms assigned in inmat.f is -99.0. elprp.f then assigns this value to props(9,*), props(13,*) and props(34,*)
call di_fgm_alphas( span, first_elem, count_alpha, snode_alpha_ij )
c **********************************************************************
c * *
c * di_fgm_alphas - retrieve nodal values of alpha assigned through *
c * 'fgm' input. *
c * *
c **********************************************************************
- check for alphas that are constant throughout element. (this is check #2 because fgm assignment of alpha places ‘-99’ in the props array.)
call di_constant_alphas( span, first_elem, count_alpha, snode_alpha_ij )
c **********************************************************************
c * *
c * di_constant_alphas - retrieve element values of alpha assigned *
c * through element definitions, and add values *
c * to nodal sums. the average nodal value is *
c * then computed in the calling routine. *
c * *
c **********************************************************************
- check for temperature-dependent material properties assigned using segmental curves. curve_set_type
curve_set_type =0, temperature-independent properties
=1, temperature-dependent properties
=2, strain-rate dependent properties
call di_seg_alpha_e_nu( curve_set, num_curves_in_set,
& first_curve, curve_set_type, span, first_elem,
& count_alpha, snode_alpha_ij,
& seg_snode_e, seg_snode_nu )
c **********************************************************************
c * *
c * di_seg_alpha_e_nu - retrieve nodal values of alpha, e and nu *
c * assigned through temperature-dependent *
c * segmental curves. to simplify logic, nodal *
c * alpha values are added to sum and averaged *
c * in calling routine. e and nu values at nodes *
c * are not averaged. *
c * *
c **********************************************************************
di_fgm_setup: Domain Integral FGM SETUP¶
Description¶
call di_fgm_setup ( 1, nonode, out, temperatures )
di_fgm_setup:
allocate data structures for two terms used in the calculation of the derivative of the stress work density. these are: nodal values of stress work density and strain.
Output:
j_data.extrap_countsINTEGER (nonode) ALLOCATABLE SAVE
j_data.swd_at_nodes(nonode) ALLOCATABLE SAVE
j_data.strain_at_nodesDOUBLE PRECISION (6,nonode) ALLOCATABLE SAVE
j_data.fgm_eLOGICAL
j_data.fgm_nuLOGICAL
Calling Tree¶
c ***************************************************************
c * *
c * -di_fgm_setup *
c * -di_nod_vals *
c * -di_extrap_to_nodes *
c * -ndpts1.f *
c * -oulg1.f (oulgf) *
c * *
c ***************************************************************
Call di_nod_vals¶
call di_nod_vals( extrap_counts, swd_at_nodes, strain_at_nodes )
build average nodal values of the strain and stress work density for all nodes in the model.
c build average nodal values of the stress work density
c (swd) and strains in the model. for each element,
c shape functions will be used to extrapolate integration-
c point values to each of the element's nodes where they
c are then averaged.
Call di_extrap_to_nodes¶
subroutine to extrapolate strain and stress work density values from integration points to nodes for one element.
c loop over element nodes. for 8 and 20-noded
c hex elements using 2x2x2 integration, we
c extrapolate to the nodes using the lagrangian
c polynomials. otherwise we just average the
c integration-point values and use that value
c at every node. warp3d uses the following
c arguments to describe elements and integration
c order:
c
c etype = 1: 20-noded brick
c etype = 2: 8-noded brick
c etype = 3: 12-noded brick
c etype = 4: 15-noded brick
c etype = 5: 9-noded brick
c
c etype = 1, int_order = 1: 27 point rule (not used)
c etype = 1, int_order = 8: 8 point rule
c etype = 1, int_order = 9: 14 point rule
c etype = 2, int_order = 1: 8 point rule
c etype = 2, int_order = 2: 6 point rule (not used)
diexp4.f: Domain Integral EXPand¶
Calling Tree¶
c ***************************************************************
c * *
c * diexp4.f *
c * -diexp4 *
c * -diadit *
c * -dibmck *
c * -diexp13 *
c * -diadit *
c * -dibmck *
c * -digete *
c * *
c ***************************************************************
c *******************************************************************
c * *
c * diadit -- add to bit map from specified item *
c * *
c * dibmck -- test for entry presence in a bit map *
c * *
c * digete -- get element edge data *
c * *
c *******************************************************************
diexp4: Domain Integral EXPand domain_type 4¶
c ***************************************************************
c * *
c * domain expand 4 - expand type 4 automatic domain *
c * *
c ***************************************************************
diexp13: Domain Integral EXPand domain_type 1-3¶
call diexp13( ring, nonode, noelem, incmap, iprops, incid,
& out, bits, q_new_map, q_old_map, q_map_len )
c ***************************************************************
c * *
c * domain expand 13 - expand type 1-3 automatic domain *
c * *
c ***************************************************************
convert list of coincident nodes at each corner node on the front to a bit map. there are 2 or three corner nodes. for 2 corner nodes we just leave one of the maps all zero which makes code logic simpler in later steps. set which positions in the list of front nodes are the corner nodes.
loop over all edges of this element. get the two structure nodes corresponding to the two corner nodes for the edge. classify the edge based on it’s connectivity to nodes listed in the three node maps. build 3 logical flags for each corner node, marking it’s appearance in the 3 node maps. based on connectivity of edge nodes on the 3 nodes lists, add a corner node to the lists.
logic tree:
0 = both corner nodes appear in the 3 node lists
0 = neither corner node appears in the 3 node lists
if corner node 1 only appears in the 3 lists,
add corner node 2 to the node list in which corner
node 1 appears.
if corner node 2 only appears in the 3 lists,
add corner node 1 to the node list in which corner
node 2 appears.
dicmj.f: Domain Integral CoMpute J¶
Description¶
domain_compute:
drive execution of element routine to compute j and i-integrals for a single domain
Calling Tree¶
c ***************************************************************
c * *
c * -dicmj.f *
c * -di_trans_nodvals *
c * -vec_ops (zero_vector.f) *
c * -dielem *
c * -dielem_a.f *
c * -dielem_b.f *
c * -dielem_c.f *
c * *
c ***************************************************************
Procedure¶
- loop over all structure elements. if element is involved in integral computations, call element dependent routine.
build copy of element coordinates, displacements, velocities, accelerations, q-values from global data. zero load terms that correspond to constrained dof. if we have temperature loadings, get temps for element nodes (node values + element uniform value). rotate displacements, velocities and accelerations from constraint-compatible to global coordinates. obtain element nodal values of stress work density and strain from structure node-value arrays.
e_coord(3,num_enodes): element nodal coordinatessnodes(num_enodes): element nodal numbersq_element(num_enodes): element nodal q-valuese_displ(3,num_enodes): element nodal displacementse_vel(3,num_enodes): element nodal velocitiese_accel(3,num_enodes): element nodal accelerationse(num_enodes): element nodal Young’s modulusnu(num_enodes): element nodal Poisson’s ratioe_force(num_enodes): element nodal forcee_node_temps(num_enodes): element nodal temperaturese_alpha_ij(6,num_enodes): element nodal thermal expansion coefficientselem_nod_swd(num_enodes): element nodal values of stress work densityelem_nod_strains(6,num_enodes): element nodal strains
if necessary, transform nodal values from constraint-compatible coordinates to global coordinates.
call di_trans_nodvals( e_displ(1,enode), trnmat(snode)%mat(1,1))
call di_trans_nodvals( e_vel(1,enode), trnmat(snode)%mat(1,1))
call di_trans_nodvals( e_accel(1,enode), trnmat(snode)%mat(1,1))
c***************************************************************
c *
c subroutine to transform a 3x1 vector of constraint- *
c compatible-coordinate-system values to global-coordinate- *
c system values. the routine premultiplies the incoming *
c constraint-compatible values by the transpose of the stored *
c global-to-constraint coordinate system rotation matrix. *
c *
c***************************************************************
- for problems with temperature loads, pull out temperatures of element nodes and set the thermal expansion coefficients at element nodes. remove reference temperatures of model nodes.
- crack-face tractions. the contribution to J arising from crack-face tractions is calculated element-by-element using the equivalent nodal loads on each element face computed during setup of load step for applied pressures, surface tractions and body forces. here we just get the element equiv. nodal forces. they are passed to element J routine which sorts out the loaded face. a similar procedure is used to compute the contribution to I.
call vec_ops( e_force,
& eq_node_forces(eq_node_force_indexes(elemno)),
& dummy, 3*num_enodes, 5 )
c ****************************************************************
c * *
c * perform simple vector-vector operations on single or *
c * double precision vectors (based on the port) *
c * *
c ****************************************************************
subroutine vec_ops( veca, vecb, vecc, n, opcode )
c opcode v: a = b
determine if we really need to process element based on q-values, nodal velocities, accelerations and specified temperature loadings.
gather element stresses, histories, properties and 3x3 rotation matrices from polar decompositions at each gauss point of element for geonl
e_stress(nstrs * num_gpts): element stresses at gauss pointselem_hist(hist_size * num_gpts): element histories at gauss pointse_props(mxelpr): element propertiese_rots(9*num_gpts): rotation matrices
gather element strains at integration points. and store in tensor form. convert engineering (total) strains of off-diagonal terms to tensor strains.
e_strain(9,num_gpts): element strains at gauss points in tensor form
for crack-face loading, make a copy of the tractions input by the user in the domain definition. if none was input, dielem_c.f automatically uses the equivalent nodal loads used to solve the boundary value problem.
cf_load(1) = cf_tractions(1)
cf_load(2) = cf_tractions(2)
cf_load(3) = cf_tractions(3)
- call the element routine to compute contribution to the j-integral and or i-integral for domain
call dielem( e_coord, q_element, e_displ, e_vel, e_accel,
& e_force, e_stress, elem_hist, e_props, e_props,
& e_props, e_rots, hist_size, domain_rot, iout,
& e_jresults, e_node_temps, elem_temps, e_alpha_ij,
& ierr, element, debug_elements, one_point_rule,
& geonl, numrow_sig, snodes, elem_nod_swd,
& elem_nod_strains, omit_ct_elem, fgm_e, fgm_nu,
& e_strain, e, e_front, nu, nu_front, front_nodes,
& num_front_nodes, front_coords, domain_origin,
& domain_type, front_order, e_iresults, comput_i,
& comput_j, cf_traction_flags, cf_load,
& front_elem_flag, expanded_front_nodes,
& myid, numprocs, crack_curvature, face_loading,
& seg_curves_flag, process_temperatures,
& max_exp_front )
- print values as required based on user specified flags.
2a. J-integral results: calculate stress intensity factor K, from J
2b. interaction-integral results: calculate K from relationship with interaction integral
c the 8 terms for I-integral results are stored
c in array e_iresults(8,7) as follows:
c
c value auxiliary field
c
c iiterms(i,1): KI plane stress
c iiterms(i,2): KI plane stress
c iiterms(i,3): KII plane stress
c iiterms(i,4): KII plane stress
c iiterms(i,5): KIII anti-plane shear
c iiterms(i,6): T11,T33 plane stress
c iiterms(i,7): T11,T33 plane strain
c iiterms(i,8): T13 anti-plane strain
2c. compute J-values from K-values
- print sum of values from all elements in domain, and J, I for domain
dielem: Domain Integral ELEMent¶
c defined in dielem_a.f
subroutine dielem ( coord, qvals, edispl, evel, eaccel,
& feload, estress, history,
& props, iprops, lprops, erots, nrow_hist,
& rotate, iout, e_jresults, e_node_temps,
& temperatures, enode_alpha_ij, ierr, elemno,
& gdebug, one_point, geonl, numrows_stress,
& snodes, elem_nod_swd, elem_nod_strains,
& omit_ct_elem, fgm_e, fgm_nu, e_strain,
& ym_nodes, ym_front_node, nu_nodes,
& nu_front_node, front_nodes, num_front_nodes,
& front_coords, domain_origin, domain_type,
& front_order, e_iresults, comput_i, comput_j,
& cf_traction_flags, cf_tractions,
& front_elem_flag, expanded_front_nodes,
& myid, numprocs, crack_curvature, face_loading,
& seg_curves_flag, process_temperatures,
& max_exp_front )
c called by dicmj
call dielem( e_coord, q_element, e_displ, e_vel, e_accel,
& e_force, e_stress, elem_hist, e_props, e_props,
& e_props, e_rots, hist_size, domain_rot, iout,
& e_jresults, e_node_temps, elem_temps, e_alpha_ij,
& ierr, element, debug_elements, one_point_rule,
& geonl, numrow_sig, snodes, elem_nod_swd,
& elem_nod_strains, omit_ct_elem, fgm_e, fgm_nu,
& e_strain, e, e_front, nu, nu_front, front_nodes,
& num_front_nodes, front_coords, domain_origin,
& domain_type, front_order, e_iresults, comput_i,
& comput_j, cf_traction_flags, cf_load,
& front_elem_flag, expanded_front_nodes,
& myid, numprocs, crack_curvature, face_loading,
& seg_curves_flag, process_temperatures,
& max_exp_front )
Description¶
- domain integral for 3-d isoparametric elements. supports:
- finite strains-rotations
- body forces (including inertia)
- crack-face tractions
- temperature loads
- kinetic energy terms
- anisotropic thermal expansion coefficients
- nonhomogeneous material properties
- interaction integral for 3-d isoparametric elements. supports:
- linear-elastic material behavior
- crack-face tractions
- nonhomogeneous material properties
c ****************************************************************
c * *
c * element name type no. description *
c * ------------ -------- ----------- *
c * *
c * q3disop 1 20 node brick(*) *
c * l3dsiop 2 8 node brick *
c * ts12isiop 3 12 node brick *
c * ts15isiop 4 15 node brick *
c * ts9isiop 5 9 node brick *
c * *
c ****************************************************************
c ****************************************************************
c * *
c * strain-stress ordering in warp3d vectors: *
c * ---------------------------------------- *
c * *
c * eps-x, eps-y, eps-z, gam-xy, gam-yz, gam-xz *
c * sig-x, sig-y, sig-z, tau-xy, tau-yz, tau-xz *
c * *
c ****************************************************************
Inputs¶
coord:e_coord(3,num_enodes): element nodal coordinatesqvals:q_element(num_enodes): element nodal q-valuesedispl:e_displ(3,num_enodes): element nodal displacementsevel:e_vel(3,num_enodes): element nodal velocitieseaccel:e_accel(3,num_enodes): element nodal accelerationsfeload:e_force(num_enodes): element nodal forceestress:e_stress(nstrs * num_gpts): element stresses at gauss pointshistory:elem_hist(hist_size * num_gpts): element histories at gauss pointsprops:e_props(mxelpr)(=props(mxelpr,elemno)): element propertiesiprops:e_props(mxelpr)(=props(mxelpr,elemno)): element propertieslprops:e_props(mxelpr)(=props(mxelpr,elemno)): element propertieserots:e_rots(9*num_gpts): 3x3 rotation matrices from polar decompositions at each gauss point of element for geonlnrow_hist:hist_size: history size at each gauss pointrotate:domain_rot(3,3): global -> crack front local rotation matrixiout:iout: output file handlee_node_temps:e_node_tempselement nodal temperaturestemperatures:elem_tempselement temperaturesenode_alpha_ij:e_alpha_ij(6,num_enodes): element nodal thermal expansion coefficientsierr:ierrerror codeelemno:elementcurrent element numbergdebug:debug_elements: debug flags for domain integral computationsone_point:one_point_rule: use one point rule flag.true.for use 1 point integration for 8-node elementsgeonl:geonl(=lprops(18,elemno)) geometric nonlinearity flagnumrows_stress:numrow_sig(=nstrs=9) number of stress values per strain point in arrays passed by warpsnodes:snodes(num_enodes): element nodal numberselem_nod_swd:elem_nod_swd(num_enodes): element nodal values of stress work densityelem_nod_strains:elem_nod_strains(6,num_enodes): element nodal strainsomit_ct_elem:omit_ct_elemomit current element flagfgm_e:fgm_e‘fgm’ values of e flag.true.for ‘fgm’ values of e have been assigned at the nodesfgm_nu:fgm_nu‘fgm’ values of nu flag.true.for ‘fgm’ values of nu have been assigned at the nodese_strain:e_strain(9,num_gpts): element strains at gauss points in tensor formym_nodes:e(num_enodes): element nodal Young’s modulusym_front_node:e_frontyoung’s modulus at point on front for homogeneous materialnu_nodes:nu(num_enodes): element nodal Poisson’s rationu_front_node:nu_frontpoisson’s ratio at point on front for homogeneous materialfront_nodes:front_nodes(num_front_nodes)numbers of front nodesnum_front_nodes:num_front_nodesnumbers of front nodesfront_coords:front_coords(3,num_front_nodes)coordinates of crack-front nodesdomain_origin:domain_originposition of domain origin in the listfront_nodesdomain_type:domain_type: domain type of q functions a-dfront_order:front_order: order of crack front interpolationcomput_i:comput_i: compute interaction inetgral flag.true.for compute interaction inetgralcomput_j:comput_j: compute j inetgral flag.true.for compute j inetgralcf_traction_flags:cf_traction_flags: crack face traction flag.true.for compute crack face load termscf_tractions:cf_load(3): crack face loadfront_elem_flag:front_elem_flag: front element flag.true.for current element is front elementexpanded_front_nodes:expanded_front_nodes:myid:myid:numprocs:numprocs:crack_curvature:crack_curvature(7): coefficients of curve described by crack front nodes. (see subroutinedi_calc_curvature)face_loading:face_loading: flagseg_curves_flag:seg_curves_flag: temperature loading flag.true.for have temperature loadingprocess_temperatures:process_temperatures: process temperatures flag.true.for temperature depended processmax_exp_front:max_exp_front(=200): maximum number of front nodes list
Outputs¶
e_jresults:e_jresults(8)e_iresults:e_iresults(8,8)
Calling Tree¶
c ****************************************************************
c * *
c * (dielem_a.f) *
c * -dielem *
c * -digetr *
c * -diqmp1 *
c * -getgpts (getgpts.f) *
c * -derivs (derivs.f) *
c * -dielcj *
c * -dieler *
c * -digrad *
c * -dielrv *
c * -dieliq *
c * -dielrt *
c * -dippie *
c * -shapef (shapef.f) *
c * -dielav *
c * -di_calc_j_terms *
c * *
c * (dielem_b.f) *
c * -di_calc_r_theta *
c * -di_calc_distance *
c * -di_calc_constitutive *
c * -di_calc_aux_fields_k *
c * -di_calc_aux_fields_t *
c * -di_calc_i_terms *
c * *
c * (dielem_c.f) *
c * -di_calc_surface_integrals *
c * -dielwf *
c * -dielrl *
c * -di_calc_surface_j *
c * -di_calc_surface_i *
c * -di_reorder_nodes *
c * *
c ****************************************************************
Procedure¶
- set up basic element properties. load all six thermal expansion coefficients from material associated with the element. if they are all zero, we skip processing of element for temperature loading. all temperature terms in domain involve derivatives of temperature at gauss points and expansion coefficients at gauss points. if specified coefficients for the element are zero, we skip temperature processing of elements.
- for ‘fgm’ alpha values, elem_alpha(1)-elem_alpha(3) will receive a value of -99.0. in the loop over integration points, subroutine dielav changes elem_alpha to the value of alpha at the integration point by interpolating nodal alpha values. to be consistent with the solution of the boundary-value problem, for linear-displacement elements, dielav will assign to elem_alpha the average of nodal alpha values.
- set integration order info for volume integrals.
- set number of stress values per strain point in arrays passed by warp
- set up data needed at all integration points and nodes needed for the volume integrals and crack face traction integrals:
- transform unrotated Cauchy stresses to Cauchy stresses for geonl formulation. estress on input contains unrotated stresses. first convert to Cauchy stresses then to 1 st PK stresses. these are in global coordinates. note: 1 PK stresses are non-symmetric so we keep the 3x3 tensor stored as a vector. position 10 is for the work density. the work density is scaled by det F to refer to t=0 configuration.
- for small-displacement formulation, copy the input (symmetric) stresses (6x1) into 3x3 tensor form.
do ptno = 1, ngpts
loc = ( ptno-1 ) * numrow
stresses(1,ptno) = estress(loc+1)
stresses(2,ptno) = estress(loc+4)
stresses(3,ptno) = estress(loc+6)
stresses(4,ptno) = estress(loc+4)
stresses(5,ptno) = estress(loc+2)
stresses(6,ptno) = estress(loc+5)
stresses(7,ptno) = estress(loc+6)
stresses(8,ptno) = estress(loc+5)
stresses(9,ptno) = estress(loc+3)
stresses(10,ptno) = estress(loc+7)
end do
- copy engineering strains into 3x3 tensor form. change gamma strains to tensor strains for off-diagonal terms. this formulation is not correct for large strains.
do enode = 1, nnode
strains(1,enode) = elem_nod_strains(1,enode)
strains(2,enode) = elem_nod_strains(4,enode)/2.0
strains(3,enode) = elem_nod_strains(6,enode)/2.0
strains(4,enode) = elem_nod_strains(4,enode)/2.0
strains(5,enode) = elem_nod_strains(2,enode)
strains(6,enode) = elem_nod_strains(5,enode)/2.0
strains(7,enode) = elem_nod_strains(6,enode)/2.0
strains(8,enode) = elem_nod_strains(5,enode)/2.0
strains(9,enode) = elem_nod_strains(3,enode)
end do
- rotate nodal coordinates, displacements, velocities and accelerations to crack front normal coordinate system.
call dielrv: Domain Integral ELement Rotate Vector
call dielrv( coord, cdispl, cvel, caccel, edispl, evel, eaccel,
& rotate, nnode, iout, debug )
c *******************************************************************
c * *
c * rotate nodal vector values to crack x-y-z *
c * *
c *******************************************************************
cdispl(3,20)nodal displacements at crack x-y-zcvel(3,20)nodal velocities at crack x-y-zcaccel(3,20)nodal accelerations at crack x-y-z
- modify the nodal q-values to reflect linear interpolations for 20, 15, 12-node elements.
call dieliq: Domain Integral ELement Interpolate Q-values
call dieliq( qvals, debug, iout, nnode, etype, elemno, snodes,
& coord, front_elem_flag, num_front_nodes, front_nodes,
& front_coords, expanded_front_nodes, domain_type,
& domain_origin, front_order, qp_node, crack_curvature,
& max_exp_front)
c ******************************************************************
c * *
c * interpolate q-values at side nodes *
c * *
c * if any 1/4-point node is detected, special integration *
c * for crack-face traction on crack-front element faces *
c * is not used. *
c * *
c ******************************************************************
- compute coordinate jacobian at all gauss points and the center point of element. these are now in cracked coordinates.
6a. isoparametric coordinates of gauss point
call getgpts( etype, order, ptno, xsi, eta, zeta, weight )
c ****************************************************************
c * *
c * given the element type, order of integration and the *
c * integration point number, return the parametric *
c * coordinates for the point and the weight value for *
c * integration *
c * *
c ****************************************************************
call derivs( etype, xsi, eta, zeta, dsf(1,1,ptno),
& dsf(1,2,ptno), dsf(1,3,ptno) )
c ****************************************************************
c * *
c * given the element type, integration point coordinates *
c * return the parametric derivatives for each node *
c * shape function *
c * *
c ****************************************************************
6b. coordinate jacobian, it’s inverse, and determinant.
call dielcj: Domain Integral ELement Compute Jacobian
call dielcj( dsf(1,1,ptno), coord, nnode, jacob(1,1,ptno),
& jacobi(1,1,ptno), detvol(ptno), ierr, iout, debug )
c *******************************************************************
c * *
c * compute coordinate jacobian, its determinate, and inverse *
c * *
c *******************************************************************
jacob(3,3,28)coordinate jacobian matrix at integration pointsjacobi(3,3,28)coordinate inverse jacobian matrix at integration pointsdetvol(28)determinant of jacobian matrix at integration points
- rotate stresses and strains to crack-front coordinates at all gauss points. get the stress work density from warp results.
call dielrt: Domain Integral ELement Rotate Tensor
call dielrt( ptno, rotate, stresses(1,ptno), csig(1,ptno),
& 1, iout, debug )
call dielrt( ptno, rotate, e_strain(1,ptno), ceps_gp(1,ptno),
& 2, iout, debug )
c *******************************************************************
c * *
c * rotate stresses or strains from global->crack x-y-z *
c * *
c *******************************************************************
csig(10,27): stress in crack-front coordinates at gauss pointsceps_gp(9,27): strain in crack-front coordinates at gauss points
- rotate strains to crack-front coordinates at all nodes.
call dielrt( enode, rotate, strains(1,enode), ceps_n(1,enode),
& 3, iout, debug )
ceps_n(9,20): nodal strain in crack-front coordinates
- process temperature (inital) strains if required. rotate the thermal expansion coefficients from from global->crack coordinates in case they are anisotropic. elem_alpha(1:6) in global is replaced by elem_alpha(1:6) in crack. rotate values that are constant over the element (elem_alpha) and values at each node of the element (enode_alpha_ij). (nodal values enable computation of the x1 derivative of the alpha_ij term in J.) for isotropic CTEs this is some unnecessary work.
call dippie( rotate, iout, debug, elem_alpha )
c *******************************************************************
c * *
c * set up to process temperature effects for domain integral *
c * *
c *******************************************************************
- for single point integration of bricks, compute average stresses, energy density, and strain at center of element. set number of gauss points to one to control subsequent loop.
- loop over all gauss points and compute each term of the domain integral for j, and the interaction integral for i.
c jterm(1) : work density
c jterm(2) : traction - displacement gradient
c jterm(3) : kinetic energy denisty
c jterm(4) : accelerations
c jterm(5) : crack face loading (handled separately)
c jterm(6) : thermal strain (2 parts. set to zero for an fgm)
c jterm(7) : stress times partial of strain wrt x1 (for fgms)
c jterm(8) : partial of stress work density wrt x1 (for fgms)
c (jterms 7 & 8 are used to replace the explicit partial
c derivative of stress work density wrt x1 (for fgms))
c
c iterm(1) : stress * derivative of aux displacement
c iterm(2) : aux stress * derivative of displacement
c iterm(3) : mixed strain energy density (aux stress * strain)
c iterm(4) : first term of incompatibility (stress * 2nd deriv
c of aux displacement
c iterm(5) : second term of incompatibility (stress * deriv
c of aux strain)
c iterm(6) : deriv of constitutive tensor * strain * aux strain
c iterm(7) : (not yet verified)
c aux stress * deriv of cte * relative change in temp
c + aux stress * cte * deriv of relative change in temp
c iterm(8) : crack face traction * aux disp. derivative
- isoparametric coordinates of gauss point and weight.
call getgpts( etype, order, ptno, xsi, eta, zeta, weight )
- evaluate shape functions of all nodes at this point. used to find value of q function at integration point, the total velocity at the point, values of acceleration at the point, and alpha values at the point.
call shapef( etype, xsi, eta, zeta, sf )
call dielav: Domain Integral ELement gAuss point Values
call dielav( sf, evel, eaccel, point_velocity, point_accel_x,
& point_accel_y, point_accel_z, qvals, point_q,
& nnode, point_temp, e_node_temps, elem_alpha,
& enode_alpha_ij, linear_displ, fgm_alphas, elemno,
& ym_nodes, nu_nodes, point_ym, point_nu, fgm_e,
& fgm_nu, coord, point_x, point_y, point_z,
& seg_curves_flag, iout )
c *******************************************************************
c * *
c * get q, velocity, accelerations, temperature and alpha at *
c * gauss point. *
c * *
c *******************************************************************
point_q: q-valuse at gauss pointpoint_velocity: total velocity at gauss pointkin_energy: kinetic energy at gauss pointpoint_accel_x, point_accel_y, point_accel_z: accelerations at gauss pointpoint_temp: temperature at gauss pointelem_alpha(6): alpha at gauss pointpoint_ym: young’s modulus at gauss pointpoint_nu: poisson’s ratio at gauss pointpoint_x, point_y, point_z: coordinates of gauss point in the local crack-front system
- for this integration point, compute displacement derivatives, q-function derivatives and temperature derivative in crack coordinate system. strains will be used in the order eps11, eps22, eps33, eps12, eps23, eps13
dux, dvx, dwx: displacement derivatives with respect to x in crack coordinate systemdqx, dqy, dqz: q-function derivatives in crack coordinate systemdtx: temperature derivative with respect to x in crack coordinate systemdalpha_x1(6): alpha derivative with respect to x in crack coordinate systemdswd_x1: stress work density derivative with respect to x in crack coordinate systemdstrain_x1(9): nodal strain derivative with respect to x in crack coordinate systemcsig_dstrain_x1:de_x1: young’s modulus derivative with respect to x in crack coordinate systemdnu_x1: poisson’s ratio derivative with respect to x in crack coordinate system
- call routines to compute domain-integral terms for current integration point.
4a. compute j terms
call di_calc_j_terms(weight, evol, jterm, ptno, dqx,
& dqy, dqz, dux, dvx, dwx, dtx, csig,
& kin_energy, rho, point_q, point_accel_x,
& point_accel_y, point_accel_z, point_temp,
& process_temperatures, elem_alpha,
& dalpha_x1, csig_dstrain_x1, dstrain_x1,
& dswd_x1, omit_ct_elem, fgm_e, fgm_nu,
& elemno, myid, numprocs, seg_curves_flag,
& debug, iout)
4b. compute i terms
- perform edi evaluations for traction loaded faces for jterm(5) and then for iterm(8).
call di_calc_surface_integrals( elemno, etype, nnode, snodes,
& feload, cf_traction_flags,
& cf_tractions, rotate, dsf, jacobi,
& cdispl, qvals, coord, front_nodes,
& front_coords, domain_type,
& domain_origin, num_front_nodes,
& front_order, ym_front_node,
& nu_front_node, comput_j,
& comput_i, jterm, iterm,
& front_elem_flag, qp_node,
& crack_curvature, face_loading, iout,
& debug)
c***********************************************************************
c *
c subroutine to drive computation of surface-traction *
c integrals *
c***********************************************************************
- save edi values in result vectors
Compute i terms¶
call
di_calc_r_theta: Domain Integral CALCulate Radius and angle THETA- for straight element edges:
compute distance from integration point to the closest line that connects two adjacent front nodes.
compute angle between integration point, line connecting front nodes, and projection of integration point onto crack plane.
- for curved element edges:
compute distance from integration point to the curve fitted through adjacent front nodes.
compute angle between integration point, curve, and projection of integration point onto crack plane.
call di_calc_r_theta( 2, front_nodes, num_front_nodes,
& front_coords, domain_type, domain_origin,
& front_order, point_x, point_y, point_z,
& elemno, ptno, r, t, crack_curvature,
& debug, iout )
c ****************************************************************
c *
c subroutine to calculate radius "r", and angle "theta" of *
c a single point, measured in polar coordinates in the *
c local crack-front system. *
c *
c ****************************************************************
r: radius “r” in polar coordinates in the local crack-front systemt: angle “theta” in polar coordinates in the local crack-front system
- call
di_calc_constitutive: Domain Integral CALCulate CONSTITUTIVE tensor components
call di_calc_constitutive( dcijkl_x1, sijkl, dsijkl_x1,
& point_ym, point_nu, de_x1, dnu_x1,
& elemno, debug, iout )
c***************************************************************
c *
c subroutine to calculate constitutive tensors for *
c interaction integral *
c *
c***************************************************************
dcijkl_x1(3): constitutive tensor derivativesijkl(3): compliance tensordsijkl_x1(3): compliance tensor derivative
c the three non-zero components of cijkl are:
c (1) lambda + 2*mu = e(1-nu)/((1+nu)(1-2nu))
c d(1)/de = (1-nu)/((1+nu)(1-2nu))
c d(1)/dnu = 2e*nu(2-nu)/((1+nu)^2*(1-2nu)^2)
c (2) lambda = e*nu/((1+nu)(1-2nu))
c d(2)/de = nu/((1+nu)(1-2nu))
c d(2)/dnu = e(1+2nu^2)/((1+nu)^2*(1-2nu)^2)
c (3) 2*mu = e/(1+nu)
c d(3)/de = 1/(1+nu)
c d(3)/dnu = -e/(1+nu)^2
c
c the three non-zero components of sijkl are:
c (1) d(1) = 1/e
c d(1)/de = -1/e^2 d(1)/dnu = 0
c (2) d(2) = -nu/e
c d(2)/de = nu/e^2 d(2)/dnu = -1/e
c (3) d(3) = (1+nu)/e
c d(3)/de = -(1+nu)/e^2 d(3)/dnu = 1/e
- call
di_calc_aux_fields_k: Domain Integral CALCulate AUXiliary FIELDS for stress intensity factors K
call di_calc_aux_fields_k( elemno, ptno, r, t, ym_front_node,
& nu_front_node, dcijkl_x1, sijkl,
& dsijkl_x1, aux_stress,
& aux_strain, daux_strain_x1,
& du11_aux, du21_aux, du31_aux,
& du111_aux, du112_aux, du113_aux,
& du211_aux, du212_aux, du213_aux,
& du311_aux, du312_aux, du313_aux,
& iout )
c***************************************************************
c *
c subroutine to calculate auxiliary stresses and displacements *
c at element integration points using the expressions obtained *
c by Williams: On the stress distribution at the base of a *
c stationary crack. Journal of Applied Mechanics 24, 109-114, *
c 1957. *
c *
c***************************************************************
aux_stress(9,8): auxiliary stresses at integration pointaux_strain(9,8): auxiliary strains at integration pointdaux_strain_x1(9,8): auxiliary strains derivativesdu111_aux(8),du112_aux(8),du113_aux(8),du211_aux(8),du212_aux(8),du213_aux(8),du311_aux(8),du312_aux(8),du313_aux(8): second derivatives of displacement (uj,1i)
- call
di_calc_aux_fields_t: Domain Integral CALCulate AUXiliary FIELDS for T-stresses
call di_calc_aux_fields_t( elemno, ptno, r, t, ym_front_node,
& nu_front_node, dcijkl_x1, sijkl,
& dsijkl_x1, aux_stress,
& aux_strain, daux_strain_x1,
& du11_aux, du21_aux, du31_aux,
& du111_aux, du112_aux, du113_aux,
& du211_aux, du212_aux, du213_aux,
& du311_aux, du312_aux, du313_aux,
& iout )
c***************************************************************
c *
c subroutine to calculate auxiliary stresses and displacements *
c at integration points for t-stresses using the expressions *
c obtained by Michell for a point load on a straight crack: *
c (see Timoshenko and Goodier, Theory of Elasticity p. 110., *
c eq. 72 for stress sigma_r.) *
c *
c***************************************************************
- call
di_calc_i_terms: Domain Integral CALCulate Interaction integral terms for stress intensity factors
call di_calc_i_terms( ptno, dqx, dqy, dqz, dux, dvx, dwx,
& dtx, csig, aux_stress, ceps_gp,
& aux_strain, dstrain_x1,
& daux_strain_x1, dcijkl_x1,
& du11_aux, du21_aux, du31_aux,
& du111_aux, du211_aux, du311_aux,
& du112_aux, du212_aux, du312_aux,
& du113_aux, du213_aux, du313_aux,
& process_temperatures, elem_alpha,
& dalpha_x1, point_temp, point_q, weight,
& elemno, fgm_e, fgm_nu, iterm,
& iout, debug)
c***************************************************************
c *
c subroutine to calculate terms of i-integral *
c *
c***************************************************************
props(mxelpr,mxel)： Element Properties¶
elprp.f
global_data.propsDeclared as: REAL (mxelpr,mxel)global_data.ipropsDeclared as: INTEGER (mxelpr,mxel)global_data.lpropsDeclared as: LOGICAL (mxelpr,mxel)
\(\begin{array}{rll} \hline \textrm{Parameter} & \textrm{Type} & \textrm{Description} \\ \hline 1 & \textrm{INTEGER} & \textrm{element type} \\ 2 & \textrm{INTEGER} & \textrm{amount of nodes in an element} \\ 3 & \textrm{INTEGER} & = 24 \times \textrm{two16} + 6 \\ 4 & \textrm{INTEGER} & \textrm{number dof each element node} \\ 5 & \textrm{INTEGER} & \textrm{integration order} \\ 6 & \textrm{INTEGER} & \textrm{amount of integration points in an element} \\ 7 & & \textrm{young's modulus} \\ 8 & & \textrm{poisson's ratio} \\ 9 & & \textrm{thermal expansion coeff.} \\ & & \alpha \textrm{for isotropic,} \alpha_x \textrm{for anisotropic} \\ 10 & & \textrm{mass density} \\ 11 & \textrm{INTEGER} & = 6 \\ 12 & \textrm{INTEGER} & \textrm{output location} \\ 13 & & \textrm{thermal expansion coeff.} \alpha_y \textrm{for anisotropic} \\ 14 & & \textrm{kinematic hardening parameter(beta)} \\ 15 & & \textrm{h-prime for linear hardening} \\ 16 & \textrm{LOGICAL} & \textrm{output format for output} \\ & & \textrm{the default value is .true. for short.} \\ 17 & \textrm{LOGICAL} & \textrm{request for linear material formulation} \\ 18 & \textrm{LOGICAL} & \textrm{geometric nonlinearity flag} \\ & & \textrm{the default value is .true. for on} \\ 19 & \textrm{LOGICAL} & \textrm{bbar flag} \\ & & \textrm{the default value is .false. for off} \\ 20 & & \textrm{viscopls m-power} \\ 21 & & \textrm{power-law hardening n-power or} \\ & & \textrm{set number for segmental stress-strain curves} \\ 22 & & \textrm{viscoplastic ref strain rate} \\ & & \textrm{(the user inputs D which the inverse of eta,} \\ & & \textrm{here we change D to eta before passing the} \\ & & \textrm{value to props(22,elem) )} \\ 23 & & \textrm{uniaxial yield stress (sig-0)} \\ 24 & & \textrm{bit 1: request for debug of material computations} \\ & & \textrm{bit 2: allow material model to cutback adaptive} \\ & & \textrm{solution due to reversed plasticity} \\ & & \textrm{bit 3: set of segmental stress-strain} \\ & & \textrm{curves specified for element's material} \\ 25 & \textrm{INTEGER} & \textrm{material model type} \\ 26 & & \textrm{initial porosity for gurson model} \\ 27 & & q_1 \textrm{for gurson model} \\ 28 & & q_2 \textrm{for gurson model} \\ 29 & & q_3 \textrm{for gurson model} \\ 30 & \textrm{INTEGER} & \textrm{bit 1: logical nucleation flag for gurson model} \\ & & \textrm{bit 2: flag if element is killable} \\ 31 & & \textrm{s_n parameter for gurson model} \\ 32 & & \textrm{eps_n parameter for gurson model} \\ 33 & & \textrm{f_n parameter for gurson model} \\ 34 & & \textrm{thermal expansion coeff.} \alpha_z \textrm{for anisotropic} \\ 35 & & \textrm{thermal expansion coeff.} \alpha_{xy} \textrm{for anisotropic} \\ 36 & & \textrm{thermal expansion coeff.} \alpha_{yz} \textrm{for anisotropic} \\ 37 & & \textrm{thermal expansion coeff.} \alpha_{xz} \textrm{for anisotropic} \\ 38 & \textrm{INTEGER} & \textrm{material number during input} \\ \end{array}\)
Catalogue of Element Numbers (types) in WARP3D¶
\(\begin{array}{rll} \hline \textrm{Type} & \textrm{Name} & \textrm{Description} \\ \hline 1 & \textrm{q3disop} & \textrm{20 node hex} \\ 2 & \textrm{l3disop} & \textrm{8 nodes hex} \\ 3 & \textrm{ts12isop} & \textrm{12 node hex transition} \\ 4 & \textrm{ts15isop} & \textrm{15 node hex transition} \\ 5 & \textrm{ts9isop} & \textrm{9 node hex transition} \\ 6 & \textrm{tet10} & \textrm{10 node tetrahedron} \\ 7 & \textrm{wedge15} & \textrm{15 node wedge element} \\ 8 & \textrm{tri6} & \textrm{6 node 2-D triangle} \\ 9 & \textrm{quad8} & \textrm{8 node 2-D quadralateral} \\ 10 & \textrm{axiquad8} & \textrm{8 node axisymmetric quadralateral} \\ 11 & \textrm{axitri6} & \textrm{6 node triangular, axisymmetric} \\ 12 & \textrm{inter_8} & \textrm{8 node quadralateral interface} \\ 13 & \textrm{tet_4} & \textrm{4 node tetrahedron} \\ 14 & \textrm{trint6} & \textrm{6 node triangular interface} \\ 15 & \textrm{trint12} & \textrm{12 node triangular interface} \\ 16-17 & \textrm{***} & \textrm{reserved to use in obtaining integration} \\ & & \textrm{points, derivatives, shape functions} \\ 18 & \textrm{bar2} & \textrm{2 node bar} \\ 19 & \textrm{link2} & \textrm{2 node link element} \\ \hline \end{array}\)
matprp(300,500)： Material Properties¶
inmat.f
main_data.matprpDeclared as: REAL (300,500)main_data.imatprpDeclared as: INTEGER (300,500)main_data.lmtprpDeclared as: LOGICAL (300,500)main_data.smatprpDeclared as: CHARACTER (len=24)(300,500)
c assign default material values. the current
c ordering of material values is:
c
c (*) 1 -- young's modulus
c (*) 2 -- poisson's ratio
c 3 -- kinematic hardening ratio (beta)
c 4 -- tangent modulus for bilinear strain
c hardening (et)
c 5 -- inviscid yield stress (sigma-o)
c (*) 6 -- thermal expansion coefficient (isotropic)
c (*) 7 -- mass density
c 8 -- linear elastic material flag
c 9 -- material model type
c = 1 vectorized linear elastic
c and invisicid plasticity model.
c isotropic/kinematic hardening
c = 2 nonlinear elastic with linear +
c power law. small strains only.
c rate independent
c = 3 general gurson/mises model including
c nucleation, linear hardening,
c power-law hardening, matrix
c viscoplasticity
c = 4 Cohesive zone models: linear elastic,
c bilinear, ramp, exponential_1 and
c exponential_2, ppr, cavitation
c = 5 advanced cyclic plasticity model
c developed by kristine cochran
c = 6 creep model
c = 7 advanced mises model with hydrogen
c effects developed by yuiming liang
c = 8 general UMAT for warp3d.
c = 10 (matching back up with file #s)
c CP model by mark messner
c = 11 interface damage model
c 10 -- viscoplastic m power
c 11 -- power-law hardening n power
c 12 -- viscoplastic reference strain rate
c 13 -- debug material model computations
c 14 -- initial porosity (f sub 0)
c 15 -- gurson model parameter q1
c 16 -- gurson model parameter q2
c 17 -- gurson model parameter q3
c 18 -- nucleation flag (true or false)
c 19 -- nucleation parameter sn
c 20 -- nucleation parameter en
c 21 -- nucleation parameter fn
c 22 -- allow material model to cut step due to
c 23 -- flag to allow crack growth element killing
c excessive reversed plasticity
c 24 -- segmental stress-strain curve logical flag
c 25 -- flag to indicate anisotropic thermal
c coefficients are defined
c 26-31 anisotropic thermal expansion coefficients
c 32 -- interface stiffness in the longitudinal direction
c 33 -- interface stiffness in the transverse direction
c 34 -- interface stiffness in the normal direction
c 35 -- critical normal stress of the interface
c 36 -- critical shear stress of the interface
c 37 -- shape parameter for bilinear and ramp
c 38 -- second ( additional) shape parameter for ramp
c 39 -- critical separation distance in sliding
c 40 -- critical separation distance in opening
c 41 -- equivalent critical separation distance
c 42 -- a ratio to determine the equivalent separation
c under mixed mode loading ( =0 => mode I )
c 43 -- a flag for identifying the interface element
c 44 -- type of interface models
c 1 - linear elastic; 2- bilinear
c 3 - ramp; 4 - exponential_1; 5 - exponential_2
c 45 -- for segmental curve model, the segmental curve
c set number
c (*) 46 -- ductile material volume fracture for fgm
c 47 -- = 0 homogeneous cohesive material
c = 1 functionally graded cohesive material
c 48 -- critical separation distance ductile (fgm)
c 49 -- critical separation distance brittle (fgm)
c 50 -- critical stress ductile (fgm)
c 51 -- critical stress brittle (fgm)
c 52 -- beta_ductile (fgm)
c 53 -- beta_brittle (fgm)
c 54 -- compression stiffness multiplier for
c cohesive materials
c 55 -- start of props for cyclic plasticity model
c see inmat_cyclic
c 70 -- start of props for yuemin liang's adv.
c mises model that includes effects of
c staturated hydrogen
c 70 - yl_1
c 71 - yl_2
c 72 - yl_3
c 73 - yl_4
c 74 - yl_5
c 75 - yl_6
c 76 - yl_7
c 77 - yl_8
c 78 - yl_9
c 79 - yl_10
c
c 80-89 creep model
c
c 90-- PPR cohesive model
c 90-98 currently used. see comments
c in inmat_cohesive routine
c 100-- local tolerance for CP NR loop
c 101-- number of crystals at e. material point (or max n)
c 102-- angle convention (1=Kocks)
c 103-- angle type (1=degree, 2=radians)
c 104-- crystal input (1=single, 2=file)
c 105-- crystal number (for single)
c 106-- crystal offset (for list)
c 107-- orientation input (1=single, 2=file)
c 108-110-- psi, theta, phi (for single)
c 111-- orientation offset (for list)
c 112-- STRING crystal list (for offset/list)
c 113-- STRING orientation list (for offset/list)
c
c 115-- macroscale material model number
c 116-118-- s vector
c 119-121-- l vector
c 122-125-- t vector
c 126-- l_s
c 127-- l_l
c 128-- l_t
c 129-- alpha_dmg
c 130-- nstacks (temp, should calculate from element sz)
c 131-- nfail ("")
c 132-- macro_sz
c 133-- cp_sz
c
c 148-- link2 x-stiffness
c 149-- link2 y-stiffness
c 150-- link2 z-stiffness
c
c 151-200 -- Abaqus compatible UMAT
c 151 - um_1
c 151 - um_2
c ...
c ...
c 200 - um_50
c
c 201-230 cavity option for cohesive material
c
c the beta_fact is used to assist in construction
c of planar models containing one layer of elements.
c the stiffness and internal forces are mutiplied
c by the beta_fact to simulate the effect of a
c non-unit thickness.
c
c (*) some material values maybe specified as having the
c value 'fgm' (a string). In such cases the user
c supplied nodal values for the model are interpolated
c at the element gauss points during execution.
Reference: