ماتریس سختی المان خرپایی با جهت‌گیری دلخواه

ماتریس سختی ساده $\frac{E A}{L} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]$ (با فرض سطح مقطع ثابت و مدول یانگ) فقط برای یک عضو خرپایی که در راستای یک محور مختصات قرار دارد معتبر است. در سازه‌های واقعی، اعضای خرپا می‌توانند در هر زاویه‌ای جهت‌گیری شوند. برای مدیریت این موضوع، باید یک ماتریس سختی استخراج کنیم که نیروها و جابجایی‌ها را در یک سیستم مختصات جهانی مشترک (x, y) به هم مرتبط کند.

روش کار به شرح زیر است: ۱. با ماتریس سختی ساده در یک سیستم مختصات محلی (x’، y’) که با عضو هم‌راستاست شروع کنید. ۲. یک ماتریس تبدیل، T، که سیستم‌های مختصات محلی و جهانی را به هم مرتبط می‌کند، توسعه دهید. ۳. از این ماتریس تبدیل برای تبدیل ماتریس سختی محلی به ماتریس سختی جهانی مطلوب استفاده کنید.

ماتریس سختی محلی (K^محلی)

برای یک عضو خرپایی که با محور محلی x’ خود هم‌راستاست، فقط می‌تواند در برابر نیروها مقاومت کند و جابجایی‌ها را فقط در امتداد آن محور تجربه کند. هر نیرو یا جابجایی عمود بر عضو (در جهت y’) با سختی صفر همراه است زیرا باعث تغییر طول نمی‌شود.

بنابراین، می‌توانیم ماتریس سختی را برای عضو در سیستم مختصات دوبعدی محلی آن که دارای چهار درجه آزادی (u₁’، v₁’، u₂’، v₂’) است، به صورت زیر بنویسیم:

این ماتریس ۴×۴ ماتریس سختی در سیستم مختصات محلی، K^محلی است.

تبدیل مختصات

اجازه دهید سیستم مختصات جهانی (x، y) به اندازه یک زاویه θ نسبت به سیستم محلی (x’، y’) چرخیده باشد.

قبلاً آموختیم که هر بردار $𝐕$ را می‌توان با استفاده از یک ماتریس چرخش دوبعدی استاندارد از یک سیستم به سیستم دیگر تبدیل کرد. رابطه بین مؤلفه‌های محلی [V_x’, V_y’] و مؤلفه‌های جهانی [V_x, V_y] به صورت زیر است: با ترانهاده کردن دو طرف، به دست می‌آید: قرار دهید c = cos(θ) و s = sin(θ). ماتریس چرخش برابر است با: $𝐑 = [\begin{matrix} c & s \\ - s & c \end{matrix}]$ می‌توانیم این را به بردارهای جابجایی در هر گره اعمال کنیم. فرض کنید جابجایی‌های جهانی q و جابجایی‌های محلی u’ باشند:

در گره ۱:
در گره ۲:

می‌توانیم این‌ها را در یک تبدیل واحد برای کل عضو ترکیب کنیم. این، ماتریس تبدیل، T را تعریف می‌کند:

یک رابطه مشابه برای بردارهای نیرو وجود دارد: که در آن F بردار نیروهای گرهی در سیستم مختصات جهانی xy و F’ بردار نیروهای گرهی در سیستم مختصات محلی x’y’ است.

چون T یک ماتریس متعامد یکه است، معکوس آن برابر با ترانهاده آن است T^-1 = T^T. بنابراین، از معادله فوق نتیجه می‌شود که: که #### استخراج ماتریس سختی جهانی (K^جهانی)

اکنون می‌توانیم ماتریس سختی را در سیستم مختصات جهانی با ترکیب این روابط استخراج کنیم:

با رابطه بنیادی در سیستم محلی شروع کنید: $محلی$
تبدیل جابجایی u’ = Tq را جایگزین کنید: $محلی$
هر دو طرف را در T^T پیش‌ضرب کنید تا نیروهای جهانی F را به دست آورید: $محلی$
تشخیص دهید که سمت چپ، بردار نیروی جهانی، F = **T^TF’ است: $𝐅 = (𝐓^{𝖳} 𝐊^{محلی} 𝐓) 𝐪$

این معادله اکنون به شکل جهانی F = K^جهانیq است. بنابراین، معادله اصلی تبدیل برای ماتریس‌های سختی را اثبات کرده‌ایم:

با انجام این ضرب ماتریسی، به شکل صریح ماتریس سختی جهانی برای یک عضو خرپایی با جهت‌گیری دلخواه می‌رسیم:

$𝐊^{global} = \frac{E A}{L} [\begin{matrix} c^{2} & c s & - c^{2} & - c s \\ c s & s^{2} & - c s & - s^{2} \\ - c^{2} & - c s & c^{2} & c s \\ - c s & - s^{2} & c s & s^{2} \end{matrix}]$

این ماتریس ۴×۴ اکنون می‌تواند در فرایند برهم‌نهی مستقیم (اسمبل) که در بخش اول توضیح داده شد، برای ساخت ماتریس سختی جهانی کل سازه خرپایی استفاده شود.

مثال زیر، فرایند کامل تبدیل مختصات و اسمبل سختی را پیاده‌سازی می‌کند.

مثال: تحلیل یک سازه خرپایی سه‌میله‌ای

خرپای زیر را در نظر بگیرید که یک نیروی افقی و یک نیروی قائم را تحمل می‌کند (شکل را ببینید). فرض کنید مدول یانگ و سطح مقطع برای همه اعضا یکسان است: E = 100 GPa، A = 200 mm².

سیستم مختصات جهانی و گره‌ها: یک سیستم مختصات جهانی $(x, y)$ با مبدأ در گره ۱ برقرار می‌کنیم.

گره ۱: (0, 0)
گره ۲: (0, 2)
گره ۳: (2, 2)

درجات آزادی جهانی (DOFs): سازه دارای ۳ گره با ۲ درجه آزادی در هر گره ( $u, v$ ) است، که در مجموع ۶ درجه آزادی جهانی را تشکیل می‌دهد. ابتدا درجات آزادی آزاد را شماره‌گذاری می‌کنیم و سپس درجات آزادی تکیه‌گاهی (محدود). این کار محاسبات بعدی را ساده‌تر می‌کند.

درجات آزادی آزاد (در گره ۳):
- $q_{1}$ : جابجایی قائم در گره ۱
- $q_{2}$ : جابجایی افقی در گره ۳
- $q_{3}$ : جابجایی قائم در گره ۳
درجات آزادی تکیه‌گاهی (در گره‌های ۱ و ۲):
- $q_{4}$ : جابجایی افقی در گره ۱
- $q_{5}$ : جابجایی افقی در گره ۲
- $q_{6}$ : جابجایی قائم در گره ۲

از آنجایی که گره ۱ روی یک غلتک و گره ۲ لولایی است، جابجایی‌های متناظر آن‌ها باید صفر باشند: $q_{4} = q_{5} = q_{6} = 0$

مرحله ۲: استخراج ماتریس سختی هر عضو در مختصات جهانی

اکنون، هر یک از سه عضو خرپایی را به‌طور جداگانه تحلیل می‌کنیم، ماتریس سختی ۴×۴ آن را در سیستم مختصات جهانی محاسبه کرده و سطرها و ستون‌ها را به‌طور صریح با شماره‌های درجه آزادی جهانی متناظر برچسب‌گذاری می‌کنیم.

عضو ① (گره‌های ۱ تا ۲):

طول: $L_{1} = 2 m$
زاویه: عضو قائم است، بنابراین $θ = 90^{\circ}$ .
کسینوس‌ها: $c = \cos (90^{\circ}) = 0$ ، $s = \sin (90^{\circ}) = 1$ .
نگاشت DOF: درجات آزادی محلی برای این عضو ( $u_{1}, v_{1}, u_{2}, v_{2}$ ) با درجات آزادی جهانی ( $q_{4}, q_{1}, q_{5}, q_{6}$ ) متناظر است.
ماتریس سختی: $\begin{array}{cc} \begin{array}{cccc} q_{4} & q_{1} & q_{5} & q_{6} \end{array} \\ 𝐊_{e}^{(1)} = \frac{E A}{L} & [\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & - 1 \\ 0 & 0 & 0 & 0 \\ 0 & - 1 & 0 & 1 \end{array}] \end{array} \begin{matrix} q_{4} \\ q_{1} \\ q_{5} \\ q_{6} \end{matrix}$

عضو ② (گره‌های ۲ تا ۳):

طول: $L_{2} = 2 m$
زاویه: عضو افقی است، بنابراین $θ = 0^{\circ}$ .
کسینوس‌ها: $c = \cos (0^{\circ}) = 1$ ، $s = \sin (0^{\circ}) = 0$ .
نگاشت DOF: درجات آزادی محلی برای این عضو ( $u_{2}, v_{2}, u_{3}, v_{3}$ ) با درجات آزادی جهانی ( $q_{5}, q_{6}, q_{2}, q_{3}$ ) متناظر است.
ماتریس سختی: $\begin{array}{cc} \begin{array}{cccc} q_{5} & q_{6} & q_{2} & q_{3} \end{array} \\ 𝐊_{e}^{(2)} = \frac{E A}{L} & [\begin{array}{cccc} 1 & 0 & - 1 & 0 \\ 0 & 0 & 0 & 0 \\ - 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{array}] \end{array} \begin{matrix} q_{5} \\ q_{6} \\ q_{2} \\ q_{3} \end{matrix}$

عضو ③ (گره‌های ۱ تا ۳):

طول: $L_{3} = \sqrt{(2 - 0)^{2} + (2 - 0)^{2}} = 2 \sqrt{2} m$
زاویه: عضو در زاویه ۴۵ درجه قرار دارد، بنابراین $θ = 45^{\circ}$ .
کسینوس‌ها: $c = \cos (45^{\circ}) = 1 / \sqrt{2}$ ، $s = \sin (45^{\circ}) = 1 / \sqrt{2}$ .
نگاشت DOF: درجات آزادی محلی برای این عضو ( $u_{1}, v_{1}, u_{3}, v_{3}$ ) با درجات آزادی جهانی ( $q_{4}, q_{1}, q_{2}, q_{3}$ ) متناظر است.
ماتریس سختی: برای این عضو، $c^{2} = 1 / 2$ ، $s^{2} = 1 / 2$ ، $c s = 1 / 2$ . $\begin{array}{cc} \begin{array}{cccc} q_{4} & q_{1} & q_{2} & q_{3} \end{array} \\ 𝐊_{e}^{(3)} = \frac{E A}{2 \sqrt{2} L} & [\begin{array}{cccc} 1 & 1 & - 1 & - 1 \\ 1 & 1 & - 1 & - 1 \\ - 1 & - 1 & 1 & 1 \\ - 1 & - 1 & 1 & 1 \end{array}] \end{array} \begin{matrix} q_{4} \\ q_{1} \\ q_{2} \\ q_{3} \end{matrix}$

مرحله ۳: اسمبل ماتریس سختی جهانی ( $𝐊$ )

اکنون ماتریس سختی جهانی ۶×۶ را با افزودن مؤلفه‌های هر ماتریس عضو در محل‌های صحیح، با هدایت برچسب‌های DOF جهانی، اسمبل می‌کنیم. ماتریس نهایی $𝐊$ حاصل جمع سه ماتریس سختی عضو است که به فضای جهانی ۶×۶ گسترش یافته‌اند.

برای مثال، بیایید جمله جهانی $K_{11}$ را پیدا کنیم که متناظر با اندرکنش بین DOFهای جهانی $q_{1}$ و $q_{1}$ است. در هر ماتریس عضو به دنبال درایه (q₁, q₁) می‌گردیم:

عضو ①: درایه در سطر و ستون $\frac{E A}{L}$ است.
عضو ②: سطر یا ستون ندارد. سهم آن در $K_{11}$ صفر است.
عضو ③: درایه در سطر $q_{1}$ و ستون $q_{1}$ برابر $\frac{E A}{2 \sqrt{2} L}$ است.
مجموع: $K_{11} = \frac{E A}{L} + 0 + \frac{E A}{2 \sqrt{2} L} = E A (\frac{1}{2} + \frac{1}{4 \sqrt{2}})$ .

به عنوان مثالی دیگر، بیایید جمله خارج قطری $K_{16}$ (سطر ۱، ستون ۶) را پیدا کنیم:

عضو ①: درایه در سطر $q_{1}$ و ستون $q_{6}$ برابر $\frac{- E A}{L}$ است.
عضو ②: DOFهای $q_{6}$ و $q_{1}$ را ندارد. سهم آن صفر است.
عضو ③: DOFهای $q_{6}$ و $q_{1}$ را ندارد. سهم آن صفر است.
مجموع: $K_{16} = - \frac{E A}{L} + 0 + 0 = - \frac{E A}{L}$ .

این فرایند افزودن مقادیر بر اساس شاخص‌های سطر و ستون جهانی آن‌ها (برهم‌نهی مستقیم) برای تمام ۳۶ درایه ماتریس سختی جهانی تکرار می‌شود.

مرحله ۴: افرازبندی و حل دستگاه

دستگاه کامل $𝐏 = Kq$ است. این ماتریس ۶×۶ اسمبل‌شده $𝐊$ منفرد است زیرا هنوز شرایط تکیه‌گاهی را اعمال نکرده‌ایم. دستگاه را با افرازبندی آن به DOFهای آزاد و تکیه‌گاهی حل می‌کنیم.

${\begin{matrix} 𝐏_{f} \\ 𝐏_{s} \end{matrix}} = [\begin{matrix} 𝐊_{f f} & 𝐊_{f s} \\ 𝐊_{s f} & 𝐊_{s s} \end{matrix}] {\begin{matrix} 𝐪_{f} \\ 𝐪_{s} \end{matrix}}$

که در آن افرازها عبارت‌اند از: * جابجایی‌های آزاد: $𝐪_{f} = {\begin{matrix} q_{1} \\ q_{2} \\ q 3 \end{matrix}}$ (این‌ها مجهولات هستند). * جابجایی‌های تکیه‌گاهی: $𝐪_{s} = {\begin{matrix} q_{4} \\ q_{5} \\ q_{6} \end{matrix}} = 0$ . * نیروهای اعمالی: $𝐏_{f} = {\begin{matrix} P_{1} \\ P_{2} \\ P_{3} \end{matrix}} = {\begin{matrix} 0 \\ 20000 \\ - 30000 \end{matrix}}$ . * واکنش‌های تکیه‌گاهی: $𝐏_{s} = {\begin{matrix} P_{4} \\ P_{5} \\ P_{6} \end{matrix}}$ (این‌ها نیز مجهول هستند).

این معادله افرازشده منفرد به دو معادله ماتریسی جداگانه گسترش می‌یابد:

$𝐏_{f} = 𝐊_{f f} 𝐪_{f} + 𝐊_{f s} 𝐪_{s}$
$𝐏_{s} = 𝐊_{s f} 𝐪_{f} + 𝐊_{s s} 𝐪_{s}$

حل برای جابجایی‌های مجهول شرط مرزی شناخته‌شده $𝐪_{s} = 0$ را اعمال می‌کنیم. این، معادله اول را بسیار ساده می‌کند: $𝐏_{f} = 𝐊_{f f} 𝐪_{f} + 𝐊_{f s} (0) ⟹ 𝐏_{f} = 𝐊_{f f} 𝐪_{f}$ ماتریس $𝐊_{f f}$ افراز ۳×۳ بالا-چپ ماتریس جهانی $𝐊$ است که متناظر با DOFهای آزاد $q_{1}$ ، $q_{2}$ و $q_{3}$ می‌باشد. اکنون می‌توانیم جابجایی‌های مجهول $𝐪_{f}$ را با معکوس‌گیری از این ماتریس کوچک و غیرمنفرد ^۱ حل کنیم: $𝐪_{f} = 𝐊_{f f}^{- 1} 𝐏_{f}$ این محاسبه مقادیر عددی جابجایی‌های افقی و قائم در گره ۳ را به دست می‌دهد.

حل برای واکنش‌های تکیه‌گاهی هنگامی که بردار جابجایی $𝐪_{f}$ معلوم شد، می‌توانیم از معادله دوم (دوباره با $𝐪_{s} = 0$ ) برای یافتن نیروهای واکنش تکیه‌گاهی مجهول استفاده کنیم: $𝐏_{s} = 𝐊_{s f} 𝐪_{f} + 𝐊_{s s} (0) ⟹ 𝐏_{s} = 𝐊_{s f} 𝐪_{f}$ این گام نهایی یک ضرب ماتریس-بردار ساده است که چهار نیروی واکنش در تکیه‌گاه‌ها (گره‌های ۱ و ۲) را به دست می‌دهد. این، تحلیل خرپا را کامل می‌کند.

کد زیر تمام محاسبات را نشان می‌دهد. نشان می‌دهد که $𝐏_{s} = {\begin{matrix} 30 \\ - 50 \\ 30 \end{matrix}} k N .$

import numpy as np
# --------------------------------------
# Step 1: Input node coordinates
# --------------------------------------
# Example geometry: a triangle truss (3 nodes, 2 elements)
nodes = {
    1: [0.0, 0.0],
    2: [0.0, 2.0],
    3: [2.0, 2.0]
}
# --------------------------------------
# Step 2: Define DOF mapping (ID matrix)
# --------------------------------------
# 4 rows (u1, v1, u2, v2) per element
# Each column corresponds to one element
# Example: two elements (1–2) and (2–3)
ID = np.array([
    [4, 5, 4],   # u1 global DOF for elem1, elem2
    [1, 6, 1],   # v1
    [5, 2, 2],   # u2
    [6, 3, 3]    # v2
])
# Number of global DOF
num_dof = int(np.max(ID))
Ks = np.zeros((num_dof, num_dof))  # global stiffness matrix initialized to zero
# --------------------------------------
# Step 3: Define truss element stiffness
# --------------------------------------
def truss_element_stiffness(E, A, node1, node2, nodes):
    """
    Compute 4x4 stiffness matrix for a 2D truss element in global coordinates.
    Parameters
    ----------
    E : float
        Young's modulus
    A : float
        Cross-sectional area
    node1, node2 : int
        Node IDs of the element
    nodes : dict
        Dictionary of node coordinates, e.g. {1: [x1, y1], 2: [x2, y2], ...}
    """
    x1, y1 = nodes[node1]
    x2, y2 = nodes[node2]
    L = np.sqrt((x2 - x1)**2 + (y2 - y1)**2)
    c = (x2 - x1) / L
    s = (y2 - y1) / L
    k = (E * A / L) * np.array([
        [ c*c,  c*s, -c*c, -c*s],
        [ c*s,  s*s, -c*s, -s*s],
        [-c*c, -c*s,  c*c,  c*s],
        [-c*s, -s*s,  c*s,  s*s]
    ])
    return k
# --------------------------------------
# Step 4: Assembly function (Python version of ADDK)
# --------------------------------------
def ADDK(Ks, Ke, ID, elem_id):
    """
    Assemble element stiffness matrix Ke into global Ks using mapping ID.
    - ID: 4 x N array (u1,v1,u2,v2 per element)
    - elem_id: which element (column index)
    """
    n = 4  # local DOF per truss element
    for i in range(n):
        for j in range(n):
            ig = ID[i, elem_id]
            jg = ID[j, elem_id]
            if ig != 0 and jg != 0:
                val = Ke[i, j]
                Ks[ig - 1, jg - 1] += val
    return Ks
    
# --------------------------------------
# Step 5: Example assembly
# --------------------------------------
E = 100e9  # Young's modulus (Pa)
A = 2e-4   # cross-sectional area (m^2)
# Element 1:
Ke1 = truss_element_stiffness(E, A, 1, 2, nodes)
Ks = ADDK(Ks, Ke1, ID, elem_id=0)
# Element 2: Node 2 → Node 3
Ke2 = truss_element_stiffness(E, A, 2, 3, nodes)
Ks = ADDK(Ks, Ke2, ID, elem_id=1)
# Element 3: Node 2 → Node 3
Ke2 = truss_element_stiffness(E, A, 1, 3, nodes)
Ks = ADDK(Ks, Ke2, ID, elem_id=2)
# --------------------------------------
# Step 6: Print the result
# --------------------------------------
np.set_printoptions(precision=3, suppress=True)
print("Global stiffness matrix Ks: (1e7)")
print(Ks/1e7)
Kff=Ks[0:3,0:3]
Kfs=Ks[0:3, 3:6]
Ksf=Ks[3:6,0:3]
Kss=Ks[3:6,3:6]
Pf=np.array([[0],[20e3],[-30e3]])
qf = np.linalg.inv(Kff) @ Pf
print(qf)
Ps =  Ksf @ qf
Ps

array([[ 30000.],
       [-50000.],
       [ 30000.]])

روش‌های بهتر و کم‌هزینه‌تری از نظر محاسباتی برای حل $𝐏_{f}$ نسبت به معکوس کردن $𝐊_{f f}$ وجود دارد.↩︎

ماتریس سختی محلی (Kمحلی)