Time-Domain Simulation of a Joined-Wing¶

In this example we use the same joined-wing model as used in the [previous example](@ref static-joined-wing), but with the following time varying loads applied at the wingtip:

A piecewise-linear load $F_L$ in the x and y-directions defined as follows:

$F_L(t) = \begin{cases} t10^6 \text{ N} & 0 \leq t \leq 0.01 \\ (0.02-t)10^6 & 0.01 \leq t \leq 0.02 \\ 0 & \text{otherwise} \end{cases}$

A sinusoidal load $F_S$ applied in the z-direction defined as follows:

$F_S(t) = \begin{cases} 0 & t \lt 0 \\ 5 \times 10^3 (1-\cos(\pi t /0.02)) \text{ N} & 0 \leq t \lt 0.02 \\ 10^4 \text{ N} & 0.02 \leq t \end{cases}$

We will also use the same compliance and mass matrix for all beams, in order to simplify the problem definition.

In [1]:

using GXBeam, LinearAlgebra

# Set endpoints of each beam
p1 = [0, 0, 0]
p2 = [-7.1726, -12, -3.21539]
p3 = [7.1726, -12,  3.21539]

Cab_1 = [
0.5         0.866025  0.0
0.836516    -0.482963  0.258819
0.224144     -0.12941   -0.965926
]

Cab_2 = [
0.5         0.866025  0.0
-0.836516    0.482963 0.258819
0.224144    -0.12941   0.965926
]

# beam 1
L_b1 = norm(p1-p2)
r_b1 = p2
nelem_b1 = 8
lengths_b1, xp_b1, xm_b1, Cab_b1 = discretize_beam(L_b1, r_b1, nelem_b1, frame=Cab_1)

# beam 2
L_b2 = norm(p3-p1)
r_b2 = p1
nelem_b2 = 8
lengths_b2, xp_b2, xm_b2, Cab_b2 = discretize_beam(L_b2, r_b2, nelem_b2, frame=Cab_2)

# combine elements and points into one array
nelem = nelem_b1 + nelem_b2
points = vcat(xp_b1, xp_b2[2:end])
start = 1:nelem
stop = 2:nelem + 1
lengths = vcat(lengths_b1, lengths_b2)
midpoints = vcat(xm_b1, xm_b2)
Cab = vcat(Cab_b1, Cab_b2)

# assign all beams the same compliance and mass matrix
compliance = fill(Diagonal([2.93944738387698e-10, 8.42991725049126e-10,
    3.38313996669689e-08, 4.69246721094557e-08, 6.79584100559513e-08,
    1.37068861370898e-09]), nelem)
mass = fill(Diagonal([4.86e-2, 4.86e-2, 4.86e-2, 1.0632465e-2, 2.10195e-4,
    1.042227e-2]), nelem)

# create assembly
assembly = Assembly(points, start, stop;
    compliance = compliance,
    mass = mass,
    frames = Cab,
    lengths = lengths,
    midpoints = midpoints)

F_L = (t) -> begin
    if 0.0 <= t < 0.01
        1e6*t
    elseif 0.01 <= t < 0.02
        -1e6*(t-0.02)
    else
        zero(t)
    end
end

F_S = (t) -> begin
    if t < 0.0
        zero(t)
    elseif 0.0 <= t < 0.02
        5e3*(1-cos(pi*t/0.02))
    else
        1e4
    end
end

# assign boundary conditions and point load
prescribed_conditions = (t) -> begin
    Dict(
        # fixed endpoint on beam 1
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        # force applied on point 4
        nelem_b1 + 1 => PrescribedConditions(Fx=F_L(t), Fy=F_L(t), Fz=F_S(t)),
        # fixed endpoint on last beam
        nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
    )
end

# time
t = range(0, 0.04, length=1001)

system, history, converged = time_domain_analysis(assembly, t;
    prescribed_conditions=prescribed_conditions,
    structural_damping=false)

nothing #hide

We can visualize tip displacements and the resultant forces accessing the post-processed results for each time step contained in the variable history. Note that the fore-root and rear-root resultant forces for this case are equal to the external forces/moments, but with opposite sign.

In [2]:

using Plots
pyplot()

point = vcat(fill(nelem_b1+1, 6), fill(1, 6))
field = [:u, :u, :u, :theta, :theta, :theta, :F, :F, :F, :M, :M, :M]
direction = [1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3]
ylabel = ["\$u_x\$ (\$m\$)", "\$u_y\$ (\$m\$)", "\$u_z\$ (\$m\$)",
    "Rodriguez Parameter \$\\theta_x\$", "Rodriguez Parameter \$\\theta_y\$",
    "Rodriguez Parameter \$\\theta_z\$", "\$F_x\$ at the forewing root (\$N\$)",
    "\$F_y\$ at the forewing root (\$N\$)", "\$F_z\$ at the forewing root (\$N\$)",
    "\$M_x\$ at the forewing root (\$Nm\$)", "\$M_y\$ at the forewing root (\$Nm\$)",
    "\$M_z\$ at the forewing root (\$N\$)"]
ph = Vector{Any}(undef, 12)

for i = 1:12
    ph[i] = plot(
        xlim = (0, 0.04),
        xticks = 0:0.01:0.04,
        xlabel = "Time (s)",
        ylabel = ylabel[i],
        grid = false,
        overwrite_figure=false
        )
    y = [getproperty(state.points[point[i]], field[i])[direction[i]] for state in history]

    if field[i] == :theta
        # convert to angle
        @. y = 4*atan(y/4)
    end

    if field[i] == :F || field[i] == :M
        y = -y
    end

    plot!(t, y, label="")
end

ph[1]

sys:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.

Out[2]:

In [3]:

ph[2]

Out[3]:

In [4]:

ph[3]

Out[4]:

In [5]:

ph[4]

Out[5]:

In [6]:

ph[5]

Out[6]:

In [7]:

ph[6]

Out[7]:

In [8]:

ph[7]

Out[8]:

In [9]:

ph[8]

Out[9]:

In [10]:

ph[9]

Out[10]:

In [11]:

ph[10]

Out[11]:

In [12]:

ph[11]

Out[12]:

In [13]:

ph[12]

Out[13]:

These graphs are identical to those presented in "GEBT: A general-purpose nonlinear analysis tool for composite beams" by Wenbin Yu and Maxwell Blair.

We can also visualize the time history of the system using ParaView. In order to view the small deflections we'll scale all the deflections up by a couple orders of magnitude. We'll also set the color gradient to match the magnitude of the deflections at each point.

In [14]:

airfoil  = [ #FX 60-100 airfoil
0000000 0.0000000;
0010700 0.0057400;
0042800 0.0114400;
0096100 0.0177500;
0170400 0.0236800;
0265300 0.0294800;
0380600 0.0352300;
0515600 0.0405600;
0669900 0.0460900;
0842700 0.0508600;
1033200 0.0556900;
1240800 0.0598900;
1464500 0.0640400;
1703300 0.0675400;
1956200 0.0708100;
2222100 0.0733900;
2500000 0.0756500;
2788600 0.0772000;
3086600 0.0783800;
3392800 0.0788800;
3705900 0.0789800;
4024500 0.0784500;
4347400 0.0775000;
4673000 0.0759600;
5000000 0.0740900;
5327000 0.0717400;
5652600 0.0691100;
5975500 0.0660800;
6294100 0.0627500;
6607200 0.0590500;
6913400 0.0551100;
7211400 0.0508900;
7500000 0.0465200;
7777900 0.0420000;
8043801 0.0374700;
8296700 0.0329800;
8535500 0.0286400;
8759201 0.0244700;
8966800 0.0205300;
9157300 0.0168100;
9330100 0.0134200;
9484400 0.0103500;
9619400 0.0076600;
9734700 0.0053400;
9829600 0.0034100;
9903900 0.0019300;
9957200 0.0008600;
9989300 0.0002300;
0000000 0.0000000;
9989300 0.0001500;
9957200 0.0007000;
9903900 0.0015100;
9829600 0.00251;
9734700 0.00377;
9619400 0.00515;
9484400 0.00659;
9330100 0.00802;
9157300 0.00941;
8966800 0.01072;
8759201 0.01186;
8535500 0.0128;
8296700 0.01347;
8043801 0.01381;
7777900 0.01373;
7500000 0.01329;
7211400 0.01241;
6913400 0.01118;
6607200 0.00951;
6294100 0.00748;
5975500 0.00496;
5652600 0.00217;
532700  -0.00092;
500000  -0.00405;
467300  -0.00731;
434740  -0.01045;
402450  -0.01357;
370590  -0.01637;
339280  -0.01895;
308660  -0.021;
278860  -0.02275;
250000  -0.02389;
222210  -0.02475;
195620  -0.025;
170330  -0.02503;
146450  -0.02447;
124080  -0.02377;
103320  -0.02246;
084270  -0.0211;
066990  -0.01913;
051560  -0.0173;
038060  -0.01481;
026530  -0.01247;
017040  -0.0097;
009610  -0.00691;
004280  -0.00436;
001070  -0.002;
0        0.0;
]

section = zeros(3, size(airfoil, 1))
for ic = 1:size(airfoil, 1)
    section[1,ic] = airfoil[ic,1] - 0.5
    section[2,ic] = 0
    section[3,ic] = airfoil[ic,2]
end

mkpath("dynamic-joined-wing-simulation")
write_vtk("dynamic-joined-wing-simulation/dynamic-joined-wing-simulation", assembly, history, t, scaling=1e2;
    sections = section)

This notebook was generated using Literate.jl.