Composite Materials in caps2tacs (#266)

* add new titanium alloy

* add smeared stringer model and PCOMP cards into caps2tacs

* update sizing optimization to allow switching to stringers

* add missing BaseProperty class check

* working smeared stringer model

* correct the smeared stringer modulus factor

* test smeared stringer on simple wing..

* working composite materials with caps2tacs

* working composite material case

* rename composite material test to #9

examples/caps_wing/3_sizing_optimization.py

@@ -25,23 +25,31 @@
     tacs_model
 )
 
+use_stringers = False
 aluminum = caps2tacs.Isotropic.aluminum().register_to(tacs_model)
+aluminum_w_stringer = caps2tacs.Orthotropic.smeared_stringer(
+    aluminum, area_ratio=0.5
+).register_to(tacs_model)
+if use_stringers:
+    material = aluminum_w_stringer
+else:
+    material = aluminum
 
 # setup the thickness design variables + automatic shell properties
 nribs = int(tacs_model.get_config_parameter("nribs"))
 nspars = int(tacs_model.get_config_parameter("nspars"))
 nOML = nribs - 1
 for irib in range(1, nribs + 1):
     caps2tacs.ThicknessVariable(
-        caps_group=f"rib{irib}", value=0.03, material=aluminum
+        caps_group=f"rib{irib}", value=0.03, material=material
     ).register_to(tacs_model)
 for ispar in range(1, nspars + 1):
     caps2tacs.ThicknessVariable(
-        caps_group=f"spar{ispar}", value=0.03, material=aluminum
+        caps_group=f"spar{ispar}", value=0.03, material=material
     ).register_to(tacs_model)
 for iOML in range(1, nOML + 1):
     caps2tacs.ThicknessVariable(
-        caps_group=f"OML{iOML}", value=0.1, material=aluminum
+        caps_group=f"OML{iOML}", value=0.1, material=material
     ).register_to(tacs_model)
 
 # add constraints and loads

examples/caps_wing/9_composite_mat_test.py

@@ -0,0 +1,148 @@
+"""
+Sean Engelstad, Febuary 2023
+GT SMDO Lab, Dr. Graeme Kennedy
+Caps to TACS example
+"""
+
+from tacs import caps2tacs
+from mpi4py import MPI
+import numpy as np
+import matplotlib.pyplot as plt
+
+# run a steady elastic structural analysis in TACS using the tacsAIM wrapper caps2tacs submodule
+# -------------------------------------------------------------------------------------------------
+# 1: build the tacs aim, egads aim wrapper classes
+
+comm = MPI.COMM_WORLD
+tacs_model = caps2tacs.TacsModel.build(csm_file="simple_naca_wing.csm", comm=comm)
+tacs_model.mesh_aim.set_mesh(
+    edge_pt_min=15,
+    edge_pt_max=20,
+    global_mesh_size=0.25,
+    max_surf_offset=0.01,
+    max_dihedral_angle=15,
+).register_to(tacs_model)
+
+carbon_fiber = caps2tacs.Orthotropic.carbon_fiber().register_to(tacs_model)
+
+# setup the thickness design variables + automatic shell properties
+nribs = int(tacs_model.get_config_parameter("nribs"))
+nspars = int(tacs_model.get_config_parameter("nspars"))
+# makes unidirectional laminate composite properties automatically
+for irib in range(1, nribs + 1):
+    caps_group = f"rib{irib}"
+    thick = 0.05
+    caps2tacs.CompositeProperty(
+        caps_group=caps_group,
+        ply_materials=[carbon_fiber] * 4,
+        ply_thicknesses=[thick / 4] * 4,
+        ply_angles=[0, -45, 45, 90],
+    ).register_to(tacs_model)
+    # caps2tacs.ThicknessVariable(
+    #     caps_group=f"rib{irib}",
+    #     value=thick,
+    # ).register_to(tacs_model)
+for ispar in range(1, nspars + 1):
+    caps_group = f"spar{ispar}"
+    thick = 0.04
+    caps2tacs.CompositeProperty(
+        caps_group=caps_group,
+        ply_materials=[carbon_fiber] * 4,
+        ply_thicknesses=[thick / 4] * 4,
+        ply_angles=[0, -45, 45, 90],
+    ).register_to(tacs_model)
+    # caps2tacs.ThicknessVariable(
+    #     caps_group=caps_group,
+    #     value=thick,
+    # ).register_to(tacs_model)
+
+caps_group = "OML"
+thick = 0.03
+caps2tacs.CompositeProperty(
+    caps_group=caps_group,
+    ply_materials=[carbon_fiber] * 4,
+    ply_thicknesses=[thick / 4] * 4,
+    ply_angles=[0, -45, 45, 90],
+).register_to(tacs_model)
+# caps2tacs.ThicknessVariable(
+#     caps_group=caps_group,
+#     value=thick,
+# ).register_to(tacs_model)
+
+# add constraints and loads
+caps2tacs.PinConstraint("root").register_to(tacs_model)
+caps2tacs.GridForce("OML", direction=[0, 0, 1.0], magnitude=100).register_to(tacs_model)
+
+# add analysis functions to the model
+caps2tacs.AnalysisFunction.ksfailure(ksWeight=50.0, safetyFactor=1.5).register_to(
+    tacs_model
+)
+caps2tacs.AnalysisFunction.mass().register_to(tacs_model)
+
+caps2tacs.ShapeVariable("rib_a1", value=1.0).register_to(tacs_model)
+
+# run the pre analysis to build tacs input files
+# alternative is to call tacs_aim.setup_aim().pre_analysis() with tacs_aim = tacs_model.tacs_aim
+tacs_model.setup(include_aim=True)
+
+# ----------------------------------------------------------------------------------
+# 2. Run the TACS steady elastic structural analysis, forward + adjoint
+
+# choose method 1 or 2 to demonstrate the analysis : 1 - fewer lines, 2 - more lines
+method = 1
+
+# show both ways of performing the structural analysis in more or less detail
+if method == 1:  # less detail version
+    tacs_model.pre_analysis()
+    tacs_model.run_analysis()
+
+    # assembler =
+    # sx = np.array([0 for _ in range(8)])
+    # sx[0] = 1.0
+    # sy = np.array([0 for _ in range(8)])
+    # sy[1] = 1.0
+
+    tacs_model.post_analysis()
+
+    print("Tacs model static analysis outputs...\n")
+    for func in tacs_model.analysis_functions:
+        print(f"func {func.name} = {func.value}")
+
+        for var in tacs_model.variables:
+            derivative = func.get_derivative(var)
+            print(f"\td{func.name}/d{var.name} = {derivative}")
+        print("\n")
+
+elif method == 2:  # longer way that directly uses pyTACS & static problem routines
+    SPs = tacs_model.createTACSProbs(addFunctions=True)
+
+    # solve each structural analysis problem (in this case 1)
+    tacs_funcs = {}
+    tacs_sens = {}
+    function_names = tacs_model.function_names
+    for caseID in SPs:
+        SPs[caseID].solve()
+        SPs[caseID].evalFunctions(tacs_funcs, evalFuncs=function_names)
+        SPs[caseID].evalFunctionsSens(tacs_sens, evalFuncs=function_names)
+        SPs[caseID].writeSolution(
+            baseName="tacs_output", outputDir=tacs_model.analysis_dir
+        )
+
+    # print the output analysis functions and sensitivities
+    print("\nTacs Analysis Outputs...")
+    for func_name in function_names:
+        # find the associated tacs key (tacs key = (loadset) + (func_name))
+        for tacs_key in tacs_funcs:
+            if func_name in tacs_key:
+                break
+
+        func_value = tacs_funcs[tacs_key].real
+        print(f"\tfunctional {func_name} = {func_value}")
+
+        struct_sens = tacs_sens[tacs_key]["struct"]
+        print(f"\t{func_name} structDV gradient = {struct_sens}")
+
+        xpts_sens = tacs_sens[tacs_key]["Xpts"]
+        print(f"\t{func_name} coordinate derivatives = {xpts_sens}")
+
+        print("\n")

examples/caps_wing/large_naca_wing.csm

@@ -43,6 +43,10 @@ set margin1 5.0
 set margin2 3.0
 set margin3 10.0
 
+# composite material coordinate system
+dimension compositeCoord 9 1 0
+set compositeCoord "0;0;0;0.0;1.0;0.0;1.0;0.0;0.0;"
+
 
 ### 1. Wing Solid Body ###
 
@@ -200,6 +204,7 @@ endif
 
 # add AIM attribute to specify the analyses to use
 select body
+csystem   wing   compositeCoord
 attribute capsAIM $egadsTessAIM;tacsAIM
 
 end

examples/caps_wing/simple_naca_wing.csm

@@ -34,6 +34,13 @@ despmtr rib_a1 1.0
 despmtr rib_a2 0.0
 set rib_a3 1.0-rib_a1-rib_a2
 
+# composite material coordinate system
+dimension YZcsys 9 1 0
+set YZcsys "0;0;0;0.0;1.0;0.0;1.0;0.0;0.0;"
+
+dimension XZcsys 9 1 0
+set XZcsys "0;0;0;0.0;1.0;0.0;1.0;0.0;0.0;"
+
 
 # Depedendent parameters
 set sspan span/2.0
@@ -62,6 +69,10 @@ translate xf -sspan zf
 rule
 
 # attribute the solid wing
+#select face
+attribute capsGroup $OML
+attribute capsBound $OML
+csystem OML YZcsys
 select face
 attribute capsGroup $OML
 attribute _color $blue
@@ -93,16 +104,17 @@ patbeg ispar nspars
    # make the panel to intersect
    box xroot margin2 -margin1 0.0 -sspan-2*margin2 2*margin1
 
-   # rotate it to match the xtip location
-   translate -xroot 0 0
-   rotatez angle 0 0
-   translate xroot 0 0
-
    # add caps attributes
-   select face
    attribute capsGroup !$spar+ispar
+   attribute capsBound $spar
+   csystem spar YZcsys
    ATTRIBUTE AFLR_GBC $TRANSP_UG3_GBC
    attribute _color $green
+
+   # rotate it to match the xtip location
+   translate -xroot 0 0
+   rotatez angle 0 0
+   translate xroot 0 0
 patend
 
 # make the ribs except for end cap ribs
@@ -119,8 +131,9 @@ patbeg index ninnerRibs
    box -margin3 ypos -margin2 2*margin3 0 2*margin2
 
    # add caps attributes
-   select face
    attribute capsGroup !$rib+irib
+   attribute capsBound $rib
+   csystem rib XZcsys
    ATTRIBUTE AFLR_GBC $TRANSP_UG3_GBC
    attribute _color $green
 
@@ -169,6 +182,7 @@ udprim editAttr filename <<
 
 # add load attribute to OML
 select face $capsGroup $OML
+# set the composite coordinate system on the OML
 #ATTRIBUTE AFLR_GBC $TRANSP_UG3_GBC
 attribute capsLoad $OML