Linear Elasticity
This tutorial illustrates the linear elasticity implementation in Modulus Sym. Modulus Sym offers the capability of solving the linear elasticity equations in the differential or variational form, allowing to solve a wide range of problems with a variety of boundary conditions. Three examples are presented in this chapter, namely the 3D bracket, the fuselage panel, and the plane displacement, to discuss the details of the linear elasticity in Modulus Sym. In this tutorial, you will learn:
How to solve linear elasticity equations using the differential and variational forms.
How to solve linear elasticity problems for 3D and thin 2D structures (plane stress).
How to nondimensionalize the elasticity equations.
This tutorial assumes that you have completed Introductory Example tutorial and have familiarized yourself with the basics of the Modulus Sym APIs. See Weak solution of PDEs using PINNs for more details on weak solutions of PDEs. Some of the boundary conditions are defined more elaborately in the tutorial Interface Problem by Variational Method .
The linear elasticity equations in Modulus Sym can be found in the
source code directory modulus/eq/pdes/linear_elasticity.py.
The Python package quadpy is required for this example.
Install using pip install quadpy (Also refer to install_modulus_docker).
Linear elasticity equations in the displacement form
The displacement form of the (steady state) linear elasticity equations, known as the Navier equations, is defined as
(149)\[(\lambda + \mu) u_{j,ji} + \mu u_{i,jj} + f_i = 0,\]
where \(u_i\) is the displacement vector, \(f_i\) is the body force per unit volume, and \(\lambda, \mu\) are the Lamé parameters defined as
(150)\[\lambda = \frac{E \nu}{(1+\nu)(1-2\nu)},\]
(151)\[\mu = \frac{E}{2(1+\nu)}.\]
Here, \(E\), \(v\) are, respectively, the Young’s modulus and Poisson’s ratio.
Linear elasticity equations in the mixed form
In addition to the displacement formulation, linear elasticity can also be described by the mixed-variable formulation. Based on internal experiments as well as the studies reported in 1, the mixed-variable formulation is easier for a neural network solver to learn, possibly due to the lower order differential terms. In the mixed-variable formulation, the equilibrium equations are defend as:
(152)\[ \sigma_{ji,j} + f_i = 0,\]
where \(\sigma_{ij}\) is the Cauchy stress tensor. the stress-displacement equations are also defined as
(153)\[\sigma_{ij} = \lambda \epsilon_{kk} \delta_{ij} + 2 \mu \epsilon_{ij},\]
where \(\delta_{ij}\) is the Kronecker delta function and \(\epsilon_{ij}\) is the strain tensor that takes the following form
(154)\[\epsilon_{ij} = \frac{1}{2}\left( u_{i,j} + u_{j,i} \right).\]
Nondimensionalized linear elasticity equations
It is often advantageous to work with the nondimensionalized and normalized variation of the elasticity equations to improve the training convergence and accuracy. The nondimensionalized variables can be defined as:
(155)\[\hat{x}_i = \frac{x_i}{L},\]
(156)\[\hat{u}_i = \frac{u_i}{U},\]
(157)\[\hat{\lambda} = \frac{\lambda}{\mu_c},\]
(158)\[\hat{\mu} = \frac{\mu}{\mu_c}\]
Here, \(L\) is the characteristic length, \(U\) is the characteristic displacement, and \(\mu_c\) is the nondimensionalizing shear modulus. The nondimensionalized Navier and equilibrium equations can be obtained by multiplying both sides of equations by \(L^2/\mu_c U\), as follows:
(159)\[(\hat{\lambda} + \hat{\mu}) \hat{u}_{j,ji} + \hat{\mu} \hat{u}_{i,jj} + \hat{f}_i = 0,\]
(160)\[\hat{\sigma}_{ji,j} + \hat{f_i} = 0,\]
where the nondimensionalized body force and stress tensor are
(161)\[\hat{f}_{i} = \frac{L^2}{\mu_c U} f_{i},\]
(162)\[\hat{\sigma}_{ij} = \frac{L}{\mu_c U}\sigma_{ij}.\]
Similarly, the nondimensionalized form of the stress-displacement equations are obtained by multiplying both sides of equations by \(L/\mu_c U\), as follows:
(163)\[ \hat{\sigma}_{ij} = \hat{\lambda} \hat{\epsilon}_{kk} \delta_{ij} + 2 \hat{\mu} \hat{\epsilon}_{ij},\]
(164)\[\hat{\epsilon}_{ij} = \frac{1}{2}\left( \hat{u}_{i,j} + \hat{u}_{j,i} \right).\]
Plane stress equations
In a plane stress setting for thin structures, it is assumed that
(165)\[\hat{\sigma}_{zz} = \hat{\sigma}_{xz} = \hat{\sigma}_{yz} = 0,\]
and therefore, the following relationship holds
(166)\[\hat{\sigma}_{zz} = \hat{\lambda} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} + \frac{\partial \hat{w}}{\partial \hat{z}} \right) + 2 \hat{\mu} \frac{\partial \hat{w}}{\partial \hat{z}} = 0 \Rightarrow \frac{\partial \hat{w}}{\partial \hat{z}} = \frac{-\hat{\lambda}}{(\hat{\lambda} +2\hat{\mu})} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} \right).\]
Accordingly, the equations for \(\hat{\sigma}_{xx}\) and \(\hat{\sigma}_{yy}\) can be updated as follows
(167)\[ \hat{\sigma}_{xx} = \hat{\lambda} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} + \frac{-\hat{\lambda}}{(\hat{\lambda} +2\hat{\mu})} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} \right) \right) + 2 \hat{\mu} \frac{\partial \hat{u}}{\partial \hat{x}}\]
(168)\[ \hat{\sigma}_{yy} = \hat{\lambda} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} + \frac{-\hat{\lambda}}{(\hat{\lambda} +2\hat{\mu})} \left(\frac{\partial \hat{u}}{\partial \hat{x}} + \frac{\partial \hat{v}}{\partial \hat{y}} \right) \right) + 2 \hat{\mu} \frac{\partial \hat{v}}{\partial \hat{y}}.\]
Problem 1: Deflection of a bracket
This linear elasticity example, shows a 3D bracket in Fig. 71. This example is partially adopted from the MATLAB PDE toolbox. The back face of this bracket is clamped, and a traction of \(4 \times 10^4 \, \text{Pa}\) is applied to the front face in the negative-\(z\) direction (this face is shown in red). The rest of the surface is considered as traction free boundaries. Other properties are set as, \((E, \nu) = (100 \,\text{GPa}, 0.3)\). You can nondimensionalize the linear elasticity equations by setting \(L=1 \,\text{m}, U=0.0001 \,\text{m}, \mu_c=0.01 \mu\).
Fig. 71 Geometry of the bracket. The back face of the bracket is clamped, and a shear stress is applied to the front face in the negative z-direction.
@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
# Specify parameters
nu = 0.3
E = 100e9
lambda_ = nu * E / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
mu_c = 0.01 * mu
lambda_ = lambda_ / mu_c
mu = mu / mu_c
characteristic_length = 1.0
characteristic_displacement = 1e-4
In general, the characteristic length can be chosen in such a way to bound the largest dimension of the geometry to \((-1,1)\). The characteristic displacement and \(\mu_c\) can also be chosen such that the maximum displacement and the applied traction are close to 1 in order. Although the maximum displacement is not known a priori, it is observed that the convergence is not sensitive to the choice of the characteristic displacement as long as it is reasonably selected based on an initial guess for the approximate order of displacement. More information on nondimensionalizing the PDEs can be found in Scaling of Differential Equations.
Case Setup and Results
The complete python script for this problem can be found at
examples/bracket/bracket.py. Two separate neural networks for displacement and stresses are used as follows
T = -4e4 * sigma_normalization
# make list of nodes to unroll graph on
le = LinearElasticity(lambda_=lambda_, mu=mu, dim=3)
disp_net = instantiate_arch(
input_keys=[Key("x"), Key("y"), Key("z")],
output_keys=[Key("u"), Key("v"), Key("w")],
cfg=cfg.arch.fully_connected,
)
stress_net = instantiate_arch(
input_keys=[Key("x"), Key("y"), Key("z")],
output_keys=[
Key("sigma_xx"),
Key("sigma_yy"),
Key("sigma_zz"),
Key("sigma_xy"),
Key("sigma_xz"),
Key("sigma_yz"),
],
cfg=cfg.arch.fully_connected,
)
nodes = (
le.make_nodes()
+ [disp_net.make_node(name="displacement_network")]
The mixed form of the linear elasticity equations is used here in this example, and therefore, the training constraints are defined as shown below:
)
# add constraints to solver
# make geometry
x, y, z = Symbol("x"), Symbol("y"), Symbol("z")
support_origin = (-1, -1, -1)
support_dim = (0.25, 2, 2)
bracket_origin = (-0.75, -1, -0.1)
bracket_dim = (1.75, 2, 0.2)
cylinder_radius = 0.1
cylinder_height = 2.0
aux_lower_origin = (-0.75, -1, -0.1 - cylinder_radius)
aux_lower_dim = (cylinder_radius, 2, cylinder_radius)
aux_upper_origin = (-0.75, -1, 0.1)
aux_upper_dim = (cylinder_radius, 2, cylinder_radius)
cylinder_lower_center = (-0.75 + cylinder_radius, 0, 0)
cylinder_upper_center = (-0.75 + cylinder_radius, 0, 0)
cylinder_hole_radius = 0.7
cylinder_hole_height = 0.5
cylinder_hole_center = (0.125, 0, 0)
support = Box(
support_origin,
(
support_origin[0] + support_dim[0],
support_origin[1] + support_dim[1],
support_origin[2] + support_dim[2],
),
)
bracket = Box(
bracket_origin,
(
bracket_origin[0] + bracket_dim[0],
bracket_origin[1] + bracket_dim[1],
bracket_origin[2] + bracket_dim[2],
),
)
aux_lower = Box(
aux_lower_origin,
(
aux_lower_origin[0] + aux_lower_dim[0],
aux_lower_origin[1] + aux_lower_dim[1],
aux_lower_origin[2] + aux_lower_dim[2],
),
)
aux_upper = Box(
aux_upper_origin,
(
aux_upper_origin[0] + aux_upper_dim[0],
aux_upper_origin[1] + aux_upper_dim[1],
aux_upper_origin[2] + aux_upper_dim[2],
),
)
cylinder_lower = Cylinder(cylinder_lower_center, cylinder_radius, cylinder_height)
cylinder_upper = Cylinder(cylinder_upper_center, cylinder_radius, cylinder_height)
cylinder_hole = Cylinder(
cylinder_hole_center, cylinder_hole_radius, cylinder_hole_height
)
cylinder_lower = cylinder_lower.rotate(np.pi / 2, "x")
cylinder_upper = cylinder_upper.rotate(np.pi / 2, "x")
cylinder_lower = cylinder_lower.translate([0, 0, -0.1 - cylinder_radius])
cylinder_upper = cylinder_upper.translate([0, 0, 0.1 + cylinder_radius])
curve_lower = aux_lower - cylinder_lower
curve_upper = aux_upper - cylinder_upper
geo = support + bracket + curve_lower + curve_upper - cylinder_hole
# Doamin bounds
bounds_x = (-1, 1)
bounds_y = (-1, 1)
bounds_z = (-1, 1)
bounds_support_x = (-1, -0.65)
bounds_support_y = (-1, 1)
bounds_support_z = (-1, 1)
bounds_bracket_x = (-0.65, 1)
bounds_bracket_y = (-1, 1)
bounds_bracket_z = (-0.1, 0.1)
# make domain
domain = Domain()
# back BC
backBC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"u": 0, "v": 0, "w": 0},
batch_size=cfg.batch_size.backBC,
lambda_weighting={"u": 10, "v": 10, "w": 10},
criteria=Eq(x, support_origin[0]),
)
domain.add_constraint(backBC, "backBC")
# front BC
frontBC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"traction_x": 0, "traction_y": 0, "traction_z": T},
batch_size=cfg.batch_size.frontBC,
criteria=Eq(x, bracket_origin[0] + bracket_dim[0]),
)
domain.add_constraint(frontBC, "frontBC")
# surface BC
surfaceBC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"traction_x": 0, "traction_y": 0, "traction_z": 0},
batch_size=cfg.batch_size.surfaceBC,
criteria=And((x > support_origin[0]), (x < bracket_origin[0] + bracket_dim[0])),
)
domain.add_constraint(surfaceBC, "surfaceBC")
# support interior
interior = PointwiseInteriorConstraint(
nodes=nodes,
geometry=geo,
outvar={
"equilibrium_x": 0.0,
"equilibrium_y": 0.0,
"equilibrium_z": 0.0,
"stress_disp_xx": 0.0,
"stress_disp_yy": 0.0,
"stress_disp_zz": 0.0,
"stress_disp_xy": 0.0,
"stress_disp_xz": 0.0,
"stress_disp_yz": 0.0,
},
batch_size=cfg.batch_size.interior_support,
bounds={x: bounds_support_x, y: bounds_support_y, z: bounds_support_z},
lambda_weighting={
"equilibrium_x": Symbol("sdf"),
"equilibrium_y": Symbol("sdf"),
"equilibrium_z": Symbol("sdf"),
"stress_disp_xx": Symbol("sdf"),
"stress_disp_yy": Symbol("sdf"),
"stress_disp_zz": Symbol("sdf"),
"stress_disp_xy": Symbol("sdf"),
"stress_disp_xz": Symbol("sdf"),
"stress_disp_yz": Symbol("sdf"),
},
)
domain.add_constraint(interior, "interior_support")
# bracket interior
interior = PointwiseInteriorConstraint(
nodes=nodes,
geometry=geo,
outvar={
"equilibrium_x": 0.0,
"equilibrium_y": 0.0,
"equilibrium_z": 0.0,
"stress_disp_xx": 0.0,
"stress_disp_yy": 0.0,
"stress_disp_zz": 0.0,
"stress_disp_xy": 0.0,
"stress_disp_xz": 0.0,
"stress_disp_yz": 0.0,
},
batch_size=cfg.batch_size.interior_bracket,
bounds={x: bounds_bracket_x, y: bounds_bracket_y, z: bounds_bracket_z},
lambda_weighting={
"equilibrium_x": Symbol("sdf"),
"equilibrium_y": Symbol("sdf"),
"equilibrium_z": Symbol("sdf"),
"stress_disp_xx": Symbol("sdf"),
"stress_disp_yy": Symbol("sdf"),
"stress_disp_zz": Symbol("sdf"),
"stress_disp_xy": Symbol("sdf"),
"stress_disp_xz": Symbol("sdf"),
"stress_disp_yz": Symbol("sdf"),
},
)
domain.add_constraint(interior, "interior_bracket")
# add validation data
mapping = {
"X Location (m)": "x",
"Y Location (m)": "y",
"Z Location (m)": "z",
"Directional Deformation (m)": "u",
}
mapping_v = {"Directional Deformation (m)": "v"}
mapping_w = {"Directional Deformation (m)": "w"}
mapping_sxx = {"Normal Stress (Pa)": "sigma_xx"}
mapping_syy = {"Normal Stress (Pa)": "sigma_yy"}
mapping_szz = {"Normal Stress (Pa)": "sigma_zz"}
mapping_sxy = {"Shear Stress (Pa)": "sigma_xy"}
mapping_sxz = {"Shear Stress (Pa)": "sigma_xz"}
mapping_syz = {"Shear Stress (Pa)": "sigma_yz"}
file_path = "commercial_solver"
if os.path.exists(to_absolute_path(file_path)):
commercial_solver_var = csv_to_dict(
to_absolute_path("commercial_solver/deformation_x.txt"),
mapping,
delimiter="\t",
)
commercial_solver_var_v = csv_to_dict(
to_absolute_path("commercial_solver/deformation_y.txt"),
mapping_v,
delimiter="\t",
)
commercial_solver_var_w = csv_to_dict(
to_absolute_path("commercial_solver/deformation_z.txt"),
mapping_w,
delimiter="\t",
)
commercial_solver_var_sxx = csv_to_dict(
to_absolute_path("commercial_solver/normal_x.txt"),
mapping_sxx,
delimiter="\t",
)
commercial_solver_var_syy = csv_to_dict(
to_absolute_path("commercial_solver/normal_y.txt"),
mapping_syy,
delimiter="\t",
)
commercial_solver_var_szz = csv_to_dict(
to_absolute_path("commercial_solver/normal_z.txt"),
mapping_szz,
delimiter="\t",
)
commercial_solver_var_sxy = csv_to_dict(
to_absolute_path("commercial_solver/shear_xy.txt"),
mapping_sxy,
delimiter="\t",
)
commercial_solver_var_sxz = csv_to_dict(
to_absolute_path("commercial_solver/shear_xz.txt"),
mapping_sxz,
delimiter="\t",
)
commercial_solver_var_syz = csv_to_dict(
to_absolute_path("commercial_solver/shear_yz.txt"),
mapping_syz,
delimiter="\t",
)
commercial_solver_var["x"] = commercial_solver_var["x"]
commercial_solver_var["y"] = commercial_solver_var["y"]
commercial_solver_var["z"] = commercial_solver_var["z"]
commercial_solver_var["u"] = (
commercial_solver_var["u"] / characteristic_displacement
)
commercial_solver_var["v"] = (
commercial_solver_var_v["v"] / characteristic_displacement
)
commercial_solver_var["w"] = (
commercial_solver_var_w["w"] / characteristic_displacement
)
commercial_solver_var["sigma_xx"] = (
commercial_solver_var_sxx["sigma_xx"] * sigma_normalization
)
commercial_solver_var["sigma_yy"] = (
commercial_solver_var_syy["sigma_yy"] * sigma_normalization
)
commercial_solver_var["sigma_zz"] = (
commercial_solver_var_szz["sigma_zz"] * sigma_normalization
)
commercial_solver_var["sigma_xy"] = (
commercial_solver_var_sxy["sigma_xy"] * sigma_normalization
)
commercial_solver_var["sigma_xz"] = (
commercial_solver_var_sxz["sigma_xz"] * sigma_normalization
)
commercial_solver_var["sigma_yz"] = (
commercial_solver_var_syz["sigma_yz"] * sigma_normalization
)
commercial_solver_invar = {
key: value
for key, value in commercial_solver_var.items()
if key in ["x", "y", "z"]
}
commercial_solver_outvar = {
key: value
for key, value in commercial_solver_var.items()
if key
in [
"u",
"v",
"w",
"sigma_xx",
"sigma_yy",
"sigma_zz",
"sigma_xy",
"sigma_xz",
"sigma_yz",
]
}
commercial_solver_validator = PointwiseValidator(
nodes=nodes,
invar=commercial_solver_invar,
true_outvar=commercial_solver_outvar,
batch_size=128,
)
domain.add_validator(commercial_solver_validator)
# add inferencer data
grid_inference = PointwiseInferencer(
nodes=nodes,
invar=commercial_solver_invar,
output_names=[
"u",
"v",
"w",
"sigma_xx",
"sigma_yy",
"sigma_zz",
"sigma_xy",
"sigma_xz",
"sigma_yz",
],
batch_size=128,
)
domain.add_inferencer(grid_inference, "inf_data")
else:
warnings.warn(
f"Directory{file_path}does not exist. Will skip adding validators. Please download the additional files from NGC https://catalog.ngc.nvidia.com/orgs/nvidia/teams/modulus/resources/modulus_sym_examples_supplemental_materials"
)
# make solver
slv = Solver(cfg, domain)
# start solver
slv.solve()
if __name__ == "__main__":
run()
The training constraints consists of two different sets of interior
points (i.e., interior_support and interior_bracket). This is
done only to generate the interior points more efficiently.
Fig. 72 shows the Modulus Sym results and also a comparison with a commercial solver results. The results of these two solvers show good agreement, with only a 8% difference in maximum bracket displacement.
Fig. 72 Modulus Sym linear elasticity results for the bracket deflection example
Problem 2: Stress analysis for aircraft fuselage panel
This example uses Modulus Sym for the analysis of stress concentration in an aircraft fuselage panel. Depending on the altitude of the flying plane, the fuselage panel is exposed to different values of hoop stress that can cause accumulated damage to the panel over the time. Therefore, in cumulative damage modeling of aircraft fuselage for the purpose of design and maintenance of aircrafts, it is required to perform several stress simulations for different hoop stress values. Consider a simplified aircraft fuselage panel as shown in Fig. 73. This example constructs a parameterized model that, once trained, can predict the stress and displacement in the panel for different values of hoop stress.
Fig. 73 Geometry and boundary conditions of the simplified aircraft fuselage panel.
Case Setup and Results
The panel material is Aluminium 2024-T3, with
\((E, \nu) = (73 \,\text{GPa}, 0.33)\). The objective is to train a parameterized
model with varying \(\sigma_{hoop} \in (46, 56.5)\). Since the
thickness of the panel is very small (i.e., \(2 \, \text{mm}\)),
plane stress equations can be used in this example. The python script
for this problem can be found at
examples/fuselage_panel/panel_solver.py. The plane stress form of
the linear elasticity equations in Modulus Sym can be called by using the
LinearElasticityPlaneStress class:
# make list of nodes to unroll graph on
le = LinearElasticityPlaneStress(lambda_=lambda_, mu=mu)
elasticity_net = instantiate_arch(
input_keys=[Key("x"), Key("y"), Key("sigma_hoop")],
output_keys=[
Key("u"),
Key("v"),
Key("sigma_xx"),
Key("sigma_yy"),
Key("sigma_xy"),
],
cfg=cfg.arch.fully_connected,
Fig. 74 shows the Modulus Sym results for panel displacements and stresses. For comparison, the commercial solver results are also presented in Fig. 75. The Modulus Sym and commercial solver results are in close agreement, with a difference of less than 5% in the maximum Von Mises stress.
Fig. 74 Modulus Sym linear elasticity results for the aircraft fuselage panel example with parameterized hoop stress. The results are for \(\sigma_{hoop}\) = 46
Fig. 75 Commercial solver linear elasticity results for the aircraft fuselage panel example with \(\sigma_{hoop}\) = 46
Linear elasticity equations in the variational form
In addition to the differential form of the linear elasticity equations, Modulus Sym also enables the use of variational form of these equations (Interface Problem by Variational Method). This section presents the linear elasticity equations in the variational form as implemented in Modulus Sym. Nondimensionalized variables will be used in this derivation. The inner product of Equation (152) and a vector test function \(v \in \mathbf{V}\), and integrating over the domain is given as follows
(169)\[\int_{\Omega} \hat{\sigma}_{ji,j} v_i d \mathbf{x}+ \int_{\Omega} \hat{f_i} v_i d \mathbf{x} = 0.\]
Using the integration by parts,
(170)\[ \int_{\partial \Omega} \hat{T_i} v_i d \mathbf{s} -\int_{\Omega} \hat{\sigma}_{ji} v_{j,i} d \mathbf{x} + \int_{\Omega} \hat{f_i} v_i d \mathbf{x} = 0,\]
where \(T_i\) is the traction. The first term is zero for the traction free boundaries. Equation (170) is the variational form of the linear elasticity equations that is adopted in Modulus Sym. The term \(\hat{\sigma}_{ji}\) in this equation is computed using the stress-displacement relation in Equation (163).
Problem 3: Plane displacement
This example solves the linear elasticity plane stress equations in the variational form. Consider a square plate that is clamped from one end and is under a displacement boundary condition on the other end, as shown in Fig. 76. The rest of the boundaries are traction free. The material properties are assumed to be \((E, \nu) = (10 \,\text{MPa}, 0.2)\). This example is adopted from 1.
Case Setup and Results
To solve this problem in the variational form, the
displacement boundary conditions are enforced in the differential form. Also
interior and boundary points are generated to be used in evaluation of the
integrals in Equation (170). The
python script for this problem can be found at
examples/plane_displacement/plane_displacement.py.
Fig. 76 Geometry and boundary conditions of the plane displacement example. This example is adopted from 1.
@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
# make list of nodes to unroll graph on
elasticity_net = instantiate_arch(
input_keys=[Key("x"), Key("y")],
output_keys=[Key("u"), Key("v")],
cfg=cfg.arch.fully_connected,
)
nodes = [elasticity_net.make_node(name="elasticity_net")]
# domain
square = Rectangle(
domain_origin,
(domain_origin[0] + domain_dim[0], domain_origin[1] + domain_dim[1]),
)
geo = square
# make domain
domain = Domain()
bottomBC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"u": 0.0, "v": 0.0},
batch_size=cfg.batch_size.bottom,
batch_per_epoch=5000,
lambda_weighting={"u": 10.0, "v": 10.0},
criteria=Eq(y, domain_origin[1]),
)
domain.add_constraint(bottomBC, "bottomBC_differential")
topBC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"u": 0.0, "v": 0.1},
batch_size=cfg.batch_size.top,
batch_per_epoch=5000,
lambda_weighting={"u": 10.0, "v": 10.0},
criteria=Eq(y, domain_origin[1] + domain_dim[1])
& (x <= domain_origin[0] + domain_dim[0] / 2.0),
)
domain.add_constraint(topBC, "topBC_differential")
# register variational data
batch_per_epoch = 1
variational_datasets = {}
batch_sizes = {}
# bottomBC, index : 0
invar = geo.sample_boundary(
batch_per_epoch * cfg.batch_size.bottom,
criteria=Eq(y, domain_origin[1]),
quasirandom=True,
)
invar["area"] *= batch_per_epoch
variational_datasets["bottom_bc"] = DictVariationalDataset(
invar=invar,
outvar_names=["u__x", "u__y", "v__x", "v__y"],
)
batch_sizes["bottom_bc"] = cfg.batch_size.bottom
# topBC, index : 1
invar = geo.sample_boundary(
batch_per_epoch * cfg.batch_size.top,
criteria=Eq(y, domain_origin[1] + domain_dim[1])
& (x <= domain_origin[0] + domain_dim[0] / 2.0),
quasirandom=True,
)
invar["area"] *= batch_per_epoch
variational_datasets["top_bc"] = DictVariationalDataset(
invar=invar,
outvar_names=["u__x", "u__y", "v__x", "v__y"],
)
batch_sizes["top_bc"] = cfg.batch_size.top
# Interior, index : 2
invar = geo.sample_interior(
batch_per_epoch * cfg.batch_size.interior,
bounds={x: bounds_x, y: bounds_y},
quasirandom=True,
)
invar["area"] *= batch_per_epoch
variational_datasets["interior"] = DictVariationalDataset(
invar=invar,
outvar_names=["u__x", "u__y", "v__x", "v__y"],
)
batch_sizes["interior"] = cfg.batch_size.interior
# make variational constraints
variational_constraint = VariationalConstraint(
datasets=variational_datasets,
batch_sizes=batch_sizes,
nodes=nodes,
num_workers=1,
loss=DGLoss(),
)
domain.add_constraint(variational_constraint, "variational")
# add inferencer data
inferencer = PointwiseInferencer(
nodes=nodes,
invar=geo.sample_interior(
2 * cfg.batch_size.interior,
bounds={x: bounds_x, y: bounds_y},
),
output_names=["u", "v"],
batch_size=2048,
plotter=InferencerPlotter(),
)
domain.add_inferencer(inferencer)
# make solver
slv = Solver(cfg, domain)
# start solver
slv.solve()
if __name__ == "__main__":
run()
The displacement boundary conditions have been included in
the regular PINN loss function. For the variational constraints, first, VariationalDataset for each bc is created by specifying the invar
and the required outvar_names. These output variables will be used while computing the variational loss. A VariationalConstraint is then constructed using the dictionary of variational_datasets and the nodes of the network. For loss, a DGLoss() is specified which is a custom loss function that includes the variational loss. The remainder
of this subsection, covers how to generate this variational
loss (DGLoss()). First, the neural network solution and the gradients at the interior and
boundary points are read:
class DGLoss(Loss):
def __init__(self):
super().__init__()
test_fn = Test_Function(
name_ord_dict={
Legendre_test: [k for k in range(10)],
Trig_test: [k for k in range(10)],
},
box=[
[domain_origin[0], domain_origin[1]],
[domain_origin[0] + domain_dim[0], domain_origin[1] + domain_dim[1]],
],
diff_list=["grad"],
)
self.v = Vector_Test(test_fn, test_fn, mix=0.02)
def forward(
self,
list_invar,
list_outvar,
step: int,
):
torch.cuda.nvtx.range_push("Make_DGLoss")
torch.cuda.nvtx.range_push("Make_DGLoss_Get_Data")
# self.v.sample_vector_test()
# get points on the interior
x_interior = list_invar[2]["x"]
y_interior = list_invar[2]["y"]
area_interior = list_invar[2]["area"]
# compute solution for the interior
u_x_interior = list_outvar[2]["u__x"]
u_y_interior = list_outvar[2]["u__y"]
v_x_interior = list_outvar[2]["v__x"]
v_y_interior = list_outvar[2]["v__y"]
# get points on the boundary
x_bottom_dir = list_invar[0]["x"]
y_bottom_dir = list_invar[0]["y"]
normal_x_bottom_dir = list_invar[0]["normal_x"]
normal_y_bottom_dir = list_invar[0]["normal_y"]
area_bottom_dir = list_invar[0]["area"]
x_top_dir = list_invar[1]["x"]
y_top_dir = list_invar[1]["y"]
normal_x_top_dir = list_invar[1]["normal_x"]
normal_y_top_dir = list_invar[1]["normal_y"]
area_top_dir = list_invar[1]["area"]
# compute solution for the boundary
u_x_bottom_dir = list_outvar[0]["u__x"]
u_y_bottom_dir = list_outvar[0]["u__y"]
v_x_bottom_dir = list_outvar[0]["v__x"]
v_y_bottom_dir = list_outvar[0]["v__y"]
u_x_top_dir = list_outvar[1]["u__x"]
u_y_top_dir = list_outvar[1]["u__y"]
In the next step, the test functions are defined, and the test functions and their required gradients are computed at the interior and boundary points:
torch.cuda.nvtx.range_pop()
torch.cuda.nvtx.range_push("Make_DGLoss_Test_Function")
# test functions
vx_x_interior, vy_x_interior = self.v.eval_test("vx", x_interior, y_interior)
vx_y_interior, vy_y_interior = self.v.eval_test("vy", x_interior, y_interior)
Here, a set of test functions consisting of Legendre polynomials and trigonometric functions is constructed, and 2% of these test functions are randomly sampled. Only the terms that appear in the variational loss in Equation (170) are computed here. For instance, it is not necessary to compute the derivative of the test functions with respect to input coordinates for the boundary points.
The next step is to compute the stress terms according to the plane stress equations in Equations (167), (168), and also the traction terms:
torch.cuda.nvtx.range_pop()
torch.cuda.nvtx.range_push("Make_DGLoss_Computation")
w_z_interior = -lambda_ / (lambda_ + 2 * mu) * (u_x_interior + v_y_interior)
sigma_xx_interior = (
lambda_ * (u_x_interior + v_y_interior + w_z_interior)
+ 2 * mu * u_x_interior
)
sigma_yy_interior = (
lambda_ * (u_x_interior + v_y_interior + w_z_interior)
+ 2 * mu * v_y_interior
)
sigma_xy_interior = mu * (u_y_interior + v_x_interior)
w_z_bottom_dir = (
-lambda_ / (lambda_ + 2 * mu) * (u_x_bottom_dir + v_y_bottom_dir)
)
sigma_xx_bottom_dir = (
lambda_ * (u_x_bottom_dir + v_y_bottom_dir + w_z_bottom_dir)
+ 2 * mu * u_x_bottom_dir
)
sigma_yy_bottom_dir = (
lambda_ * (u_x_bottom_dir + v_y_bottom_dir + w_z_bottom_dir)
+ 2 * mu * v_y_bottom_dir
)
sigma_xy_bottom_dir = mu * (u_y_bottom_dir + v_x_bottom_dir)
w_z_top_dir = -lambda_ / (lambda_ + 2 * mu) * (u_x_top_dir + v_y_top_dir)
sigma_xx_top_dir = (
lambda_ * (u_x_top_dir + v_y_top_dir + w_z_top_dir) + 2 * mu * u_x_top_dir
)
sigma_yy_top_dir = (
lambda_ * (u_x_top_dir + v_y_top_dir + w_z_top_dir) + 2 * mu * v_y_top_dir
)
sigma_xy_top_dir = mu * (u_y_top_dir + v_x_top_dir)
traction_x_bottom_dir = (
sigma_xx_bottom_dir * normal_x_bottom_dir
+ sigma_xy_bottom_dir * normal_y_bottom_dir
)
traction_y_bottom_dir = (
sigma_xy_bottom_dir * normal_x_bottom_dir
+ sigma_yy_bottom_dir * normal_y_bottom_dir
)
traction_x_top_dir = (
sigma_xx_top_dir * normal_x_top_dir + sigma_xy_top_dir * normal_y_top_dir
)
traction_y_top_dir = (
Finally, following the Equation (170), variational interior and boundary integral terms are defined and the variational loss is formulated and added to the overall loss as follows:
)
torch.cuda.nvtx.range_pop()
torch.cuda.nvtx.range_push("Make_DGLoss_Integral")
interior_loss = tensor_int(
area_interior,
sigma_xx_interior * vx_x_interior
+ sigma_yy_interior * vy_y_interior
+ sigma_xy_interior * (vx_y_interior + vy_x_interior),
)
boundary_loss1 = tensor_int(
area_bottom_dir,
traction_x_bottom_dir * vx_bottom_dir
+ traction_y_bottom_dir * vy_bottom_dir,
)
boundary_loss2 = tensor_int(
area_top_dir,
traction_x_top_dir * vx_top_dir + traction_y_top_dir * vy_top_dir,
)
torch.cuda.nvtx.range_pop()
torch.cuda.nvtx.range_push("Make_DGLoss_Register_Loss")
losses = {
"variational_plane": torch.abs(
interior_loss - boundary_loss1 - boundary_loss2
)
.pow(2)
.sum()
}
torch.cuda.nvtx.range_pop()
Fig. 77 shows the Modulus Sym results for this plane displacement example. The results are in good agreement with the FEM results reported in 1, verifying the accuracy of the Modulus Sym results in the variational form.
Fig. 77 Modulus Sym results for the plane displacement example.
References